Efficient linear scaling mapping for permutation symmetric Fock spaces
EEfficient linear scaling mapping for permutation symmetric Fock spaces
M. Ahsan Zeb ∗ Department of Physics, Quaid-i-Azam University, Islamabad 45320, Pakistan
Abstract
Numerically solving a second quantised many-body model in the permutation symmetric Fock space canbe challenging for two reasons: ( i ) an increased complication in the calculations of the matrix elements ofvarious operators, and ( ii ) a poor scaling of the cost of these calculations with the Fock space size. Wepresent a method that solves both these problems. We find a mapping that can be used to simplify thecalculations of the matrix elements. The mapping is directly generated so its computational cost scales onlylinearly with the space size and is negligible even for large enough sizes that approach the thermodynamiclimit. Keywords:
Order N, second quantised, Fock space, permutation symmetric.
PROGRAM SUMMARY
Program Title: FockMapLicensing provisions: GPLv3Programming language: FORTRANNature of problem: Solving second quantised many-bodymodels in permutation symmetric Fock spaceSolution method: A mapping between the Fock statesexists that can be used to calculate the matrix elements ofvarious operators. The pattern in the mapping is foundand it is directly generated.
Solving a model in a Fock space on computer re-quires computation of matrix elements of the Hamil-tonian and other operators. A particular indexingsystem for the basis states can be more efficient thanothers [1]. However, since the size of the Hilbert spaceof a many-body system scales exponentially with thesystem size, only the relevant subspace of the fullHilbert space is often used. For example, to calculateground state properties or absorption spectrum of a ∗ Corresponding author.
E-mail address: [email protected] system of coupled excitons and phonons, the permu-tation symmetric subspace does the job. The mainadvantage is that the size of the permutation sym-metric subspace scales only polynomially with thesystem size. However, the computation of matrix el-ements becomes non-trivial and a brute force calcu-lation of these scales as the square of the subspacesize, which in return limits the size of the systemsthat can be dealt with this method. We find thatthis problem boils down to the calculation of a certainmapping. We present an efficient method to calculatethis mapping. Our method is based on the recogni-tion of the pattern that this mapping acquires whenthe indexing of the basis states follows a certain or-der. This method can be adopted to solve any secondquantised model in the permutation symmetric Fockspace whenever the large system sizes are to be stud-ied. For example, in Ref. [2], we apply this methodon a complex many-body problem with three typesof excitations.Here we illustrate the method for a generic many-body system with N bosonic modes. In Sec. 1, we in-troduce the permutation symmetric subspace of theFock space and a specific indexing scheme for its Preprint submitted to Computer Physics Communications February 2, 2021 a r X i v : . [ phy s i c s . c o m p - ph ] J a n tates. Section 2 defines the mapping discussed aboveand Sec. 3 gives two examples of how this mappingcan be used to calculate the matrix elements of oper-ators. We describe the pattern this mapping acquiresand a possible algorithm to generate it in Sec. 4.At the end, in Sec. 5., we describe the input andoutput of “FockMap” code that efficiently calculatesthis mapping and some other properties of the basisstates, which can be employed to solve any many-body problem in the permutation symmetric Fockspace. The code is open source so it can also beintegrated into other complex programs.
1. Permutation symmetric Fock space
Consider the Fock space of N identical bosonmodes, e.g., N identical harmonic oscillators. Wecan make subsets of the Fock states such that eachsubset contains the states that are related by permu-tation of their occupation numbers. If we denote { ν } as the set of the occupation numbers for such a sub-set, then the permutation symmetric superposition ofthe states in it S N ( { ν } ), is given by, |S N { ν }(cid:105) ≡ (cid:80) P P[ |{ ν }(cid:105) ] (cid:112) P N ( { ν } ) , (1)where the sum over P indicates a sum over the permu-tations, and P N ( { ν } ) counts the number of distinctpermutations or the size of the subset, which will de-pend on the pattern of occupations in { ν } . The set ofall possible permutation symmetric states {S N ( { ν } ) } spans the permutation symmetric subspace of the Fockspace. If we label the frequency f ν i as the number oftimes each value ν i appears in the set { ν } , then thenumber of permutations is the multinomial coefficient P N ( { ν } ) = N ! / ( (cid:89) i f ν i !) . (2)For example, for the set of occupations { } , thefrequencies are 1 , , P ( { } ) = 12, andthe permutation symmetric state is: |S { }(cid:105) ≡ (cid:0) | (cid:105) + | (cid:105) + | (cid:105) + | (cid:105) + | (cid:105) + | (cid:105) + | (cid:105) + | (cid:105) + | (cid:105) + | (cid:105) + | (cid:105) + | (cid:105) (cid:1) √ . To perform numerical calculations, the occupationnumbers need to be truncated. That is, we need tointroduce a cutoff M , such that ν i ∈ [0 , M ]. In anymodel, M needs to be sufficiently large for the resultsto be converged and hence reliable. The total num-ber of distinct permutation symmetric states for N modes is M + N C M compared to a total of ( M + 1) N states, which increases only polynomially with N ,much slower than the exponential scaling of the fullHilbert space. The counting comes from the numberof ways to pick N numbers in the range [0 , M ] ignor-ing order. This far better scaling makes it possibleto calculate the lowest energy eigenstate and someother properties for large values of N, M to see thebehaviour of the model under study in the thermo-dynamic limit.
To use the the permutation symmetric basis states {S N ( { ν } ) } on computer, a suitable indexing is re-quired. That is, an integer I N ( { ν } ) for every S N ( { ν } ). There is no unique way to do it, however,there are important advantages if we choose to lex-icographically order the occupations in the set { ν } and index the basis states in order of increasing theoccupation from left to right, as shown in Table 1. Table 1: Indexing of the basis states with lexicographicallyordered sets of occupation numbers, illustrated for N = 5 , M =2. { ν } I N ( { ν } ) { } { } { } { } { } { } { } { }
7. . . . . . { } . Mapping for the permutation symmetricFock states The downside of using the permutation symmetricspace is that the calculation of the matrix elementsof the Hamiltonian and other operators becomes non-trivial. To this end, as we show in the following sec-tion (3), we find that we can always depend on a(many-to-one) mapping M N − (cid:55)→ N from the permu-tation symmetric states of N − N modes. The mapping M N − (cid:55)→ N takes us from theindex I N − ( { µ } ) and the occupation of a single mode n (which can itself be treated as an index) to I N ( { ν } ) if adding n to the set { µ } makes the set { ν } . For theindexing described in sec. 1.1, not only I N − ( { µ } )but also M N − (cid:55)→ N turns out to work for all possi-ble cases we could be interested in with permuta-tion symmetric states of N, N − , N − , ..., M N − (cid:55)→ N can also be used for M p − (cid:55)→ p with p = 1 , , ..., N .For N modes, a naive calculation of M N − (cid:55)→ N should scale as the square of the size N N of the per-mutation symmetric Fock space, i.e., as O ( N N ), asthere are N N × M integers (in the space of N modes)to compare to N N − × M × M integers. (Somerestrictions can be imposed to improve the scalingthough.) But, thanks to our indexing, M N − (cid:55)→ N in-herits a pattern that can be easily recognised andgenerated directly thus completely avoiding this bot-tleneck. This direct generation of M N − (cid:55)→ N is notonly an O ( N ) process, it has a very small prefactorsuch that the computational cost even for N ∼ is negligible.
3. Using the Mapping: Matrix elements of op-erators
Let’s calculate the matrix elements of a few opera-tors to illustrate how the mapping M N − (cid:55)→ N can beused. In this section, we discuss how the mapping M N − (cid:55)→ N can help us determine the reduced den- sity matrix. We can write eigenstate as follows: | Ψ (cid:105) = (cid:88) k N ≡I N ( { ν } ) ψ k N |S N { ν }(cid:105) , (3)where { ψ k N } is the array of coefficients that is com-puted. A crucial step to calculating observables is todefine the reduced density matrix ρ r that describesthe (mixed) state of a single mode. It requires takinga trace over the states of all modes other than theone in question. This can be written as: ρ rm,m (cid:48) ≡ (cid:104) Ψ | m (cid:48) (cid:105) (cid:104) m | Ψ (cid:105) (4)Here i is an arbitrary mode (since states are per-mutation symmetric), and m, m (cid:48) denote occupationnumber states on the molecule in question.To find the element ρ rm,m (cid:48) , we need to trace outthe states of the N − N excited modes which arereduced to the same N − m, m (cid:48) are taken out. If we denote k N = I N ( { ν } ) and k (cid:48) N = I N ( { ν (cid:48) } ) as the indices of a pair of states { ν } and { ν (cid:48) } of N modes that reduce to the same state { ν (cid:48)(cid:48) } of N − j N − = I N − ( { ν (cid:48)(cid:48) } ), wecan write k N = M N − (cid:55)→ N ( m, j N − ) ,k (cid:48) N = M N − (cid:55)→ N ( m (cid:48) , j N − ) . (5)With these maps, we can then trace over j N − , de-scribing the state of the other modes. Assuming theabove relations between j N − , k N , and k (cid:48) N , the re-duced density matrix takes the form: ρ rm,m (cid:48) = N N − (cid:88) j N − =1 ψ k N ψ ∗ k (cid:48) N P N − ( j N − ) (cid:112) P N ( k N ) P N ( k (cid:48) N ) . (6)Here N N − is the total number of the permutationalsymmetric Fock states involving N − P N − ( j N − ) counts the numberof matching terms in the permutation symmetric su-perposition of the Fock states k N , k (cid:48) N — so give unitoverlap — after taking out the states of our subject3olecule. Using Eq. 2 and keeping in mind the rela-tionship between states j N − , k N , k (cid:48) N , we can simplifyEq. 6 to ρ rm,m (cid:48) = N N − (cid:88) j N − =1 ψ k N ψ ∗ k (cid:48) N N (cid:112) f m ( { ν } ) f m (cid:48) ( { ν (cid:48) } ) . (7)Algorithm 1 (with f m,k N ≡ f m ( { ν } ), etc.) sum-marises this computation. Algorithm 1
Matrix elements of ρ r function getrho for j N − in range(0, N N − − ) dofor m in range(0,M) dofor m (cid:48) in range(0,M) do k N = M N − (cid:55)→ N ( m, j N − ) k (cid:48) N = M N − (cid:55)→ N ( m (cid:48) , j N − ); x = ψ kN ψ ∗ k (cid:48) N N (cid:113) f m,kN f m (cid:48) ,k (cid:48) N ρ rm,m (cid:48) = ρ rm,m (cid:48) + x end forend forend forreturn ρ r end function (cid:80) Ni ˆ b † i Consider the operator (cid:80) Ni ˆ b † i . The matrix elementcan be written explicitly as a sum over permutations: (cid:104){ ν (cid:48) } N | N (cid:88) j =1 ˆ b † j |{ ν } N (cid:105) = (cid:88) P (cid:48) P [ (cid:104) ν (cid:48) ν (cid:48) . . . ν (cid:48) N | ] P N ( { ν (cid:48) } ) N (cid:88) j =1 ˆ b † j (cid:88) P P [ | ν ν . . . ν N (cid:105) ] (cid:112) P N ( { ν } ) . (8)Consider (cid:80) Nj =1 ˆ b † j (cid:80) P P [ | ν ν . . . ν N (cid:105) ]. Each permu-tation gives terms such as √ ν + 1 times the statewith ν → ν + 1. In general this leads to overlaps ofthe form: √ ν i + 1 × P [ (cid:104) ν (cid:48) ν (cid:48) . . . ν (cid:48) N | ] P [ | ν i → ν i + 1 in { ν } N (cid:105) ] , which are non-zero only if { ν (cid:48) } is the same as { ν } ex-cept ν i → ν i + 1. In other words, the only differencebetween these two states is that their frequencies of ν i and ν i + 1 are different but still related by f ν i ( { ν } ) = f ν i ( { ν (cid:48) } ) + 1 , (9a) f ν i +1 ( { ν } ) = f ν i +1 ( { ν (cid:48) } ) − . (9b)If so, every ket in the permutations finds its dual.Since, there are P p ( { ν } ) permutations of { ν } , we willget P p ( { ν } ) ×√ ν i + 1 for one such term. The element ν i may occur multiple times in the set { ν } ; we denotethe frequency with which it occurs as f ν i ( { ν } ). Thematrix element then becomes (cid:104){ ν (cid:48) } N | N (cid:88) j =1 ˆ b † j |{ ν } N (cid:105) = (cid:115) P N ( { ν } ) P N ( { ν (cid:48) } ) f ν i ( { ν } ) √ ν i + 1= (cid:112) ( ν i + 1) f ν i ( { ν } )( f ν i +1 ( { ν } ) + 1) , (10)where we have used the definition of P N ( { ν } ), Eq. 2.To use the above expression in a computer pro-gram, we need to know the indexes I N ( { ν } ) and I N ( { ν (cid:48) } ) of the two states |S N { ν }(cid:105) and |S N { ν (cid:48) }(cid:105) in-volved. To determine that we can use the conditionsin Eq. 9. Suppose we have already indexed statesas in Sec. 1.1 and calculated the frequencies f ν i ( { ν } )(and P N ( { ν } ), etc.). Now, here comes the fun part.Eq. 9 can be satisfied if we start with a given permu-tation symmetric state |S N − { ν (cid:48)(cid:48) }(cid:105) of N − { ν (cid:48)(cid:48) } being the common subset of the two sets { ν } , { ν (cid:48) } , which can in fact be any set! The mapping M N − (cid:55)→ N takes us from the index I N − ( { ν (cid:48)(cid:48) } ) and ν i to I N ( { ν } ). Similarly, it gives us I N ( { ν (cid:48) } ) from I N − ( { ν (cid:48)(cid:48) } ) and ν i + 1. In summary, Algorithm 2 isquite an efficient way to calculate these matrix ele-ments.
4. Calculating the mapping M N − (cid:55)→ N In the follwing, we first describe the pattern thatthe mapping M N − (cid:55)→ N acquires if we use the index-ing scheme presented in sec. 1.1. A simple method togenerate M N − (cid:55)→ N will be presented afterwards.4 lgorithm 2 Matrix elements of ˆ H = (cid:80) i ˆ b i function getH for j N − in range(0, N N − − ) dofor m in range(0,M) do k N = M N − (cid:55)→ N ( m, j N − ) k (cid:48) N = M N − (cid:55)→ N ( m + 1 , j N − ); x = (cid:113) ( m + 1) f m,k N ( f m +1 ,k (cid:48) N + 1) H k N ,k (cid:48) N = H k N ,k (cid:48) N + x end forend forreturn H end function M N − (cid:55)→ N If we represent the mapping M N − (cid:55)→ N for a fixedM, say M = 5, and a few values of N (say N =2 , , , d arrays with the indices of N − N th modeas the row and column indices, and the indices of N mode states as array elements, we observe a clearrecursive pattern with some symmetries and well de-fined features. The map M N − (cid:55)→ N can be dividedin blocks corresponding to addition of each site, andeach block consists of a series of square transpose-symmetric sub-blocks with sizes following a recursivepattern.To illustrate this, the mapping for N = 3 and M = 5 is shown in Fig.1. Here, rows 0 − N = 2. and there is a series of smallersub-blocks in the second block, i.e., rows 6 −
20. Togenerate this map, we can divide it into these twoblocks. The first block is simple to generate. If wesee the red shaded triangular region, we find that theindex simply starts from 0 on the top left, increasesone by one as we move to the right, and leaves i column for the ith row due to the condition that theoccupation for the second site should not be less thanthat of the first one. The pattern gets a little com-plicated in the second block, rows 6 −
20. First, seethe green shaded regions. These are a set of trian-gles with smaller and smaller sizes, again due to therestriction on the occupation of the third site. Theseshaded regions are the part of the map that can berecorded when the basis are formed. Our main task I m Figure 1: Map for N=3, M=5. The map can be dividedinto two blocks, rows 0 − −
20. The latter can begenerated with a recursive function of depth 5. Fig. 2 containsthe additional blocks if we increase N by 1, i.e., N = 4. is to calculate the unshaded part of the map. Thepattern that the unshaded regions follow is also easyto see, however. In the first block, it’s simply thesymmetric image of the upper triangular part. Inthe second block, the lower triangular parts of thesquare blocks around the recursion of triangles fol-low the same rule. The one extra complication is thecolumns on the left side not included in the squareregions, but, they also follow a regular pattern. Thefirst column starting from the row 6 continues thecount of the first column from the first block until theend of the second block. The other columns startingfrom row 11 , , , ,
20, follow the same rule.The depth of the recursion for the triangles in thefigure for N = 3 is d = 5. For N = 4, there arerecursions of depth 5 , , , , N , there are recur-5 Figure 2: Additional blocks in the mapping for N=4 (andM=5), with recursive depths d = 4 , , , sions of order d, d − , d − , ..., d in N − Generating the N = 2 map is trivial. We will firstdescribe here how the N > N = 3 case shown in Fig.1, a recursivefunction can generate the second part of the map,i.e., rows 6 −
21, taking M = 5, the depth of therecursion d = 5, initial index i = 21, and the start-ing index of the column on the left i c = 6 as inputs.For N >
3, we have to generate recursions of lowerorder/depths that would require starting values for multiple columns on the left side, M − d values fororder d recursion, to be specific. In summary, wecan divide the map into blocks corresponding to re-cursions of various orders and a function can generateeach block given the required starting indices and therecursion depth.These blocks can be generated sequentially or inparallel, and combined with the block for two sites.The sequential implementation is simpler, as all thearguments of the recursive function — starting indexfor the first triangle, the depth of the recursion, andthe set of starting indices for the leftover columns —are available at each call. However, we use this onlyfor small system sizes. For larger systems, the codeuses a parallel implementation. There are two pointsto consider for parallelisation. First, the argumentsof the recursive functions need to be calculated be-forehand, and second, the workload needs to be dis-tributed evenly between all processes. The argumentsare calculated sequentially, starting from the set forthe first function call and, using the information onhow much the indices are going to advance in thatcall, calculating the arguments for the next functioncall. To distribute the load evenly, we calculate theblock sizes (heights, ( d ( d + 1) / d ) that every recursion would produce. Us-ing this list, the total size among all processes canbe evenly divided. This way, different processes canhave different number of function calls but the to-tal workload determined by the total size of the mapthey produce is approximately the same. Since, weonly fill in a grid with indices that are obtained by ei-ther adding an integer to another or just plain rangesof integers, this method does not just scale linearlywith the basis size, the prefactor or the slope of thescaling is also very small, which makes it computa-tionally extraordinarily efficient.
5. Input and output of FockMap program
The number of modes N and the cutoff on the oc-cupations M are given as the input. The programcalculates the mapping M N − (cid:55)→ N as well as someother states related properties that could be requiredto compute the matrix elements of various operators.It writes the output as unformatted binary files that6an be read by the user’s programs. For the outputdata structure, see the routines writebasis() and writebasisf() in the source file basis.f. For testingpurposes, the same routines can be modified by theuser to write the output in a human readable format.Some additional properties of the states thatthe code calculates follow. It gives the totalnumber of excitations in a state, i.e., n ( { ν } ) ≡ (cid:80) i ν i , the number of permutations P N ( k N ), andthe frequencies f ν i ( { ν } ) for any number of modesup to given N . It also calculates the ratios r ν i ( { ν (cid:48) } ) ≡ P N − ( j N − ) / P N ( k N ) for ν i , j N − ≡I N − ( { ν (cid:48) } ) , k N ≡ I N ( { ν } ) satisfying k N = M N − (cid:55)→ N ( ν i , j N − ), which can be used instead ofbare numbers P N ( k N ) , P N − ( j N − ) to increase theefficiency in some cases.
6. How to use FockMap program?
FockMap uses a “Makefile” for the compilation soif you have gfortran or some other fortran compiler(we have tested gfortran only) installed, FockMap canbe compiled simply by running make from the sourcedirectory on the command line on Unix or terminal on Mac.
On the command line, FockMap can be run like anyother serial program, i.e., by executing the command fockmap . The program asks for the input
N, M (onthe standard output), you give desired integer valuesto these and hit
Enter . The program writes the out-put in binary files basis.dat and map.dat as well asa human readable file basis-f.dat just for inspec-tion.
FockMap is open source. The source code ofFockMap can be easily integrated into other pro-grams, especially if they are written in Fortran. Thesource code is also fairly understandable so the algo-rithms used can also be copied to write the routinesin other languages.
7. Conclusion
A highly efficient method to solve many-body mod-els in the second quantised form in permutation sym-metric Fock space is presented. To facilitate the useof the method, a fortran code FockMap is providedthat can facilitate efficient calculations of any suchmodel.
References [1] A. I. Streltsov, O. E. Alon, L. S. Cederbaum,General mapping for bosonic and fermionicoperators in fock space, Phys. Rev. A 81 (2010)022124. doi:10.1103/PhysRevA.81.022124 .URL https://link.aps.org/doi/10.1103/PhysRevA.81.022124 [2] M. A. Zeb, P. G. Kirton, J. Keeling, Incoherentcharge transport in an organic polariton conden-sate (2020). arXiv:2004.09790arXiv:2004.09790