Path-integral Monte Carlo simulation of time-reversal noninvariant bulk systems with a case study of rotating Yukawa gases
PPath-integral Monte Carlo simulation of time-reversal noninvariant bulk systems witha case study of rotating Yukawa gases
Tam´as Haidekker Galambos , , and Csaba T˝oke , BME-MTA Exotic Quantum Phases “Lend¨ulet” Research Group,Budapest University of Technology and Economics,Institute of Physics, Budafoki ´ut 8, H-1111 Budapest, Hungary Department of Theoretical Physics,Budapest University of Technology and Economics,Institute of Physics, Budafoki ´ut 8, H-1111 Budapest, Hungary and Department of Physics, University of Basel, Klingelbergstrasse 82, CH-4056 Basel, Switzerland (Dated: October 9, 2018)We elaborate on the methodology to simulate bulk systems in the absence of time-reversal symme-try by the phase-fixed path-integral Monte Carlo method under (possibly twisted) periodic boundaryconditions. Such systems include two-dimensional electrons in the quantum Hall regime and rotat-ing ultracold Bose and Fermi gases; time-reversal symmetry is broken by an external magnetic fieldand the Coriolis force, respectively. We provide closed-form expressions in terms of Jacobi ellipticfunctions for the thermal density matrix (or the Euclidean propagator) of a single particle on a flattorus under very general conditions. We then modify the multi-slice sampling method in order tosample paths by the magnitude of the complex-valued thermal density matrix. Finally, we demon-strate that these inventions let us study the vortex melting process of a two-dimensional Yukawagas in terms of the de Boer interaction strength parameter, temperature, and rotation (Coriolisforce). The bosonic case is relevant to ultracold Fermi-Fermi mixtures of widely different massesunder rotation.
PACS numbers:
I. INTRODUCTION
The path-integral Monte Carlo (PIMC) method [1] letsus simulate many-body systems at finite temperature in acontrolled manner. Equilibrium properties are obtainedfrom the many-body density matrix ρ ( R, R (cid:48) ; β ) = (cid:88) n e − β(cid:15) n Ψ n ( R )Ψ ∗ n ( R (cid:48) ) , (1)where R ≡ ( r , r , . . . , r N ) collects dN particle coordi-nates, d is the dimensionality of the system, N is thenumber of particles, { Ψ n } is a complete set of many-body eigenstates, and { (cid:15) n } are the corresponding ener-gies. The convolution identity of the density matrix, ρ ( R, R (cid:48) ; β + β ) = (cid:90) dR (cid:48)(cid:48) ρ ( R, R (cid:48)(cid:48) ; β ) ρ ( R (cid:48)(cid:48) , R (cid:48) ; β ) , (2)is applied iteratively to yield the imaginary-time path-integral representation ρ ( R, R (cid:48) ; β ) = (cid:90) dR · · · (cid:90) dR M − ρ ( R, R ; τ ) × ρ ( R , R ; τ ) . . . ρ ( R M − , R (cid:48) ; τ ) . (3)Here, the time-step τ ≡ β/M corresponds to a muchhigher temperature than the system temperature. Thehigh-temperature density matrix that connects adjacentslices, ρ ( R m − , R m ; τ ), can be approximated by severalplausible schemes [1]. Estimators of physical quanti-ties are defined by integrals that involve ρ ( R, R (cid:48) ; β ); inmost cases the diagonal element ρ ( R, R ; β ) is sufficient.The Metropolis-Hastings Monte Carlo method [2, 3] is applicable to path integration if the product of high-temperature density matrices in Eq. (3) can be inter-preted as a probability density function.For time-reversal invariant bosonic systems this alwaysholds, and PIMC is an unbiased and essentially exactmethod in this case. For fermions, however, the noto-rious sign problem arises, because the contribution of aparticular path can have either sign due to the presenceof nondiagonal factors ρ ( R m − , R m ; τ ) in the integrand ofestimators. The generic means to overcome this problem,the use of restricted or constrained paths that avoid thenodal surfaces of a preconceived trial many-body densitymatrix [4, 5], makes PIMC variational in character.On the other hand, if time-reversal is not a sym-metry of the system, either because charged particlesare exposed to an external magnetic field or the sys-tem is rotated, the density matrices are complex-valuedin general and hence the prescription of the nodal sur-faces is insufficient. A consequent method would be tosample paths by the probability density function (PDF) (cid:81) Mm =1 | ρ ( R m − , R m ; τ ) | (here we assume integration withthe diagonal density matrix as the kernel of the estima-tor, and we define R = R (cid:48) ≡ R = R M ), and sum themup with the complex phase factor (cid:81) Mm =1 ρ ( R m − ,R m ; τ ) | ρ ( R m − ,R m ; τ ) | .This procedure would result in a more severe form of thesign problem: contributions with different phase factorswould cancel almost completely. This issue is equally se-vere for bosons and fermions, and it arises even in thenonphysical case of distinguishable particles (“bolzman-nons”). In analogy to the phase-fixing extension [6, 7]of zero-temperature methods such as diffusion quantumMonte Carlo [8], phase fixing is an obvious route to a r X i v : . [ c ond - m a t . s t a t - m ec h ] F e b adapt PIMC to such problems. Unlike the case of zero-temperature methods, the function whose phase needs tobe fixed is the many-body density matrix in Eq. (1), nota wave function. While the fixed-phase extension of thePIMC method is often mentioned in the literature [9], itis hardly ever applied, in contrast to the similar extensionof zero-temperature methods [6, 7, 10, 11].We address several issues related to the use of PIMCin time-reversal non-invariant bulk systems. (Finite sys-tems such as quantum dots are not our primary interesthere.) First, if we want to simulate bulk systems con-sequently, we have to use periodic boundary conditions,possibly with twist angles that let us reduce finite-size ef-fects such as shell effects in finite-size representations ofFermi liquids [12], which have analogs in strongly corre-lated electron systems in magnetic fields [13]. One shouldbase any PIMC simulation on the single-particle thermaldensity matrix (equivalently, kinetic action) that is exactunder the chosen boundary conditions. We show thatthe free propagation of a charged particle (equivalently,the thermal density matrix) on a flat torus subjected toa perpendicular magnetic field already exhibits a ratherrich structure, although these patterns lose their signif-icance for small imaginary times or large system sizes.This result lets us define the kinetic action in a way thatis compatible with the torus.The PIMC method is applicable beyond toy modelsonly because the sampling of paths could be made ef-ficient by the introduction of multi-slice moves. Thesereplace entire segments of the path [14] according tothe PDF (cid:81) Mm =1 ρ ( R m − , R m ; τ ). If, however, the densitymatrix is complex-valued and the probability density ofpaths is determined by its magnitude, the familiar bisec-tion method [1] that relies on the L´evy construction of aBrownian bridge, runs into difficulties because the convo-lution property in Eq. (2) is not applicable to magnitudes.We elaborate on a modification of the multi-slice movealgorithm that takes the external magnetic field and theperiodicity of the torus into account.Finally, we demonstrate the use of phase-fixed PIMCfor bulk systems in a case study of rotating two-dimensional Yukawa gases. Yukawa bosons arise eitherin type-II superconductors, where the Abrikosov vortexlines interact by a repulsive modified-Bessel-function po-tential ∝ K ( r ) [15–17], or in strongly interacting Fermi-Fermi mixtures of ultracold atoms, if the mass ratio ofthe two species, M/m , is very far from unity and the mo-tion of both species is confined to two dimensions [18].A flux density can be introduced to cold atomic systemsby rotating the gas, a technique that has been appliedfrequently in the past two decades [19–22]. In the modelwe consider particles that interact via a modified-Bessel-function potential ∝ K ( r ). This is a good approxima-tion also to the inter-atomic interaction in a Fermi-Fermimixture at sufficiently long range [18]. We do not claim,however, to represent either problem faithfully: we do notinclude the nonuniversal short-range repulsion betweenFermi-Fermi bound states, and the inclusion of additional flux density would be difficult to justify for Abrikosovvortices. We have deliberately chosen this system forcomputational convenience in order to demonstrate theadequacy of our methodology. On the one hand, K ( r )is mildly divergent at short range, thus even the simplestapproximation to the high-temperature density matrix,the primitive action, is a reasonable starting point. Onthe other, as K ( r ) decays exponentially at large range,the intricacies of Ewald summation can be avoided.As a first approach, we use the density matrix of thefree Bose and Fermi gases to fix the phase of the many-body density matrix. We are encouraged in this by thefact that in the case of the node fixing problem, whicharises analogously for time-reversal invariant fermionicsystems, significant progress was possible both for He[23] and the hydrogen plasma [24, 25] using the nodalsurfaces of either the noninteracting system or some well-tested variational ground state wave function. (The twoapproaches are somewhat complementary.) Simple as itis, we demonstrate that phase-fixed PIMC captures thecrystallization of rotating Yukawa bosons and fermionsas a function of interaction strength, flux density, andtemperature. We emphasize that unlike for the diffusionMonte Carlo or Green’s Function Monte Carlo methods,no trial wave function of the proper symmetry serves asinput to such a calculation; but we do choose the aspectratio of the unit cell so that it can accommodate a finitepiece of a triangular lattice.The paper is structured as follows. In Sec. II we presentthe density matrix for a single particle in a magneticfield on the torus, with some mathematical details of thederivation delegated to Appendix A, and the consider-ations of its efficient computation to Appendix B. Theadaptation of the multi-slice sampling algorithm is dis-cussed in Sec. III, with a detour to periodic, but time-reversal-invariant systems. Sec. IV presents a case study,where the phase-fixed path-integral Monte Carlo methodis applied to rotating systems of two-dimensional Yukawagases under periodic boundary conditions. In Sec. V wesummarize our results and discuss further research direc-tions. Appendix C presents the technical details of thephase-fixing methodology for PIMC.
II. THE THERMAL DENSITY MATRIX
We consider a flat torus pierced by a perpendicu-lar magnetic field. Consider the parallelogram spannedby two nonparallel vectors L = ( L ,
0) and L =( L cos θ, L sin θ ). A torus is obtained by identifying theopposite sides of this unit cell; cf. Fig. 1(a). We will referto a similar parallelogram that has the origin as its centeras the principal domain.We use the Landau gauge A = − By ˆx throughout thisarticle. Electrons are characterized by complex coordi-nates z = x + iy , and we define τ = L L e iθ , (4) (0,0) L L xy L (1+ τ )/2L (1- τ )/2-L (1+ τ )/2-L (1- τ )/2 θ (a) (0,0) L (1+ τ )L (1- τ )-L (1+ τ ) -L (1- τ ) L τ -L τ -L L θ (b) FIG. 1: (a) The principal domain of the torus. We alsodepict L , L and θ as defined in the text; we identify theplane with the complex plane, and indicate the corners of theprincipal domain using the complex parameter τ defined inEq. (4). (b) The quadruple domain used for finding the zerosof the density matrix in the low-temperature limit. so that L and L τ span the parallelogram on the com-plex plane. In the presence of a perpendicular magneticfield, magnetic translations [26] are useful: t ( L ) = exp (cid:18) i (cid:126) L · p − i ˆz · ( L × r ) (cid:96) (cid:19) , (5)where p = (cid:126) i ∇ − e A . In the current gauge, these act as t ( L ) ψ ( r ) = exp( ix ˆy · L (cid:96) ) ψ ( r + L ). We will require eachstate and the implied density matrix to obey twistedboundary conditions with twist angles φ , , t ( L , ) ψ ( r ) = e iφ , ψ ( r ) . (6)The two conditions are mutually compatible only if theparallelogram is pierced by an integral number of fluxquanta, N φ = | L × L | π(cid:96) = L L sin θ π(cid:96) . (7)Then the principal domain is also a magnetic unit cell.If N φ (cid:60) τ = k is an integer, i.e., L cos θ = kL N φ , (8)straightforward but tedious algebra yields the single-particle density matrix ρ PBC ( r , r (cid:48) ; β ) = 1 N φ ρ open ( r , r (cid:48) ; β ) × N φ − (cid:88) m =0 (cid:26) ϑ (cid:20) a m (cid:21) (cid:16) z (cid:12)(cid:12)(cid:12) τ (cid:17) ϑ (cid:20) b (cid:48) m (cid:21) ( z | τ )++( − k ϑ (cid:20) a m + (cid:21) (cid:16) z (cid:12)(cid:12)(cid:12) τ (cid:17) ϑ (cid:20) b (cid:48) m (cid:21) ( z | τ ) (cid:27) , (9)where we have factored out ρ open , the density matrix foropen boundary conditions: ρ open ( r , r (cid:48) ; β ) = 12 π(cid:96) √ u − u × exp (cid:32) − u − u | r − r (cid:48) | (cid:96) + i ( x (cid:48) − x )( y + y (cid:48) )2 (cid:96) (cid:33) , (10) where (cid:96) = (cid:113) (cid:126) eB is the magnetic length, u = e − β (cid:126) ω c , and ω c = eBm is the cyclotron frequency [27]. Above, we haveused Jacobi elliptic functions with characteristics [28, 29] ϑ (cid:20) ab (cid:21) ( z | τ ) = (cid:88) n e iπτ ( n + a ) +2 i ( n + a )( z + bπ ) . (11)The arguments in Eq. (9) are defined as τ = iπ (cid:18) L (cid:96)N φ (cid:19) u − u ,z = L (cid:96) N φ (cid:18) y + y (cid:48) + i ( x (cid:48) − x ) 1 + u − u (cid:19) ,τ = iπ (cid:18) (cid:96)N φ L (cid:19) u − u ,z = N φ πL (cid:18) x + x (cid:48) + i ( y − y (cid:48) ) 1 + u − u (cid:19) ; (12)and the constants related to boundary conditions are a m = φ πN φ + mN φ ,b m = − φ π − N φ (cid:60) τ ,b (cid:48) m = b m + N φ a m (cid:60) τ. (13)The derivation of Eq. (9) is delegated to Appendix A. -6-4-2 0 2 4 6-6 -4 -2 0 2 4 6(a) y / l x/l 1e-050.00010.0010.01 -6-4-2 0 2 4 6-6 -4 -2 0 2 4 6(b) y / l x/l 1e-050.00010.0010.01 -6-4-2 0 2 4 6-6 -4 -2 0 2 4 6(c) y / l x/l 1e-050.00010.0010.01 -6-4-2 0 2 4 6-6 -4 -2 0 2 4 6(d) y / l x/l 1e-050.00010.0010.01 FIG. 2: The dependence of | ρ PBC ( r , r (cid:48) ; β ) | on imaginary time β . There are N φ = 6 flux quanta in the principal domain, L /L = 1 . θ ≈ ◦ , φ = φ = 0, and we have fixed r (cid:48) at the origin. The panels correspond to β (cid:126) ω c = 0 .
3, 0.7, 1.1,and 5, respectively.
The behavior of the density matrix is shown in Fig. 2for the most general case, an oblique unit cell. For smallimaginary time (high temperature) | ρ PBC ( r , r (cid:48) ; β ) | has asmall Gaussian peak around r (cid:48) , which is fixed at the ori-gin in the figure. This peak spreads out by diffusion as β is increased, and eventually the Gaussians from neigh-boring unit cells start to overlap appreciably. However,the density matrix also has a phase due to the externalmagnetic field, which gives rise to an interference patternin this time range. There is destructive interference atcertain points, which effectively arrests the diffusion. [Wewill analyze the zeros of ρ PBC ( r , r (cid:48) ; β ) below.] Beyond acertain value of β , the picture is essentially stationary. We note that | ρ PBC ( r , r (cid:48) ; β ) | is not invariant for a si-multaneous displacement of both r and r (cid:48) by the samevector d , which corresponds to choosing a shifted mag-netic unit cell on the plane for compactification by pe-riodic boundary conditions, except for special choices of d . This is understood easily by noting that the secondcharacteristic b m appears in Eq. (11) as a simple additiveconstant to the variable z , letting us rewrite Eq. (9) as ρ PBC ( r , r (cid:48) ; β ) = 1 N φ ρ open ( r , r (cid:48) ; β ) × N φ − (cid:88) m =0 (cid:26) ϑ (cid:20) a m + L π(cid:96) N φ ( y + y (cid:48) ) (cid:21) (cid:18) πN φ τ (cid:48) L ( x (cid:48) − x ) (cid:12)(cid:12)(cid:12) τ (cid:48) (cid:19) ϑ (cid:20) b (cid:48) m + N φ L ( x + x (cid:48) ) (cid:21) (cid:18) π ( y − y (cid:48) ) τ L sin θ (cid:12)(cid:12)(cid:12) τ (cid:19) ++( − k ϑ (cid:20) a m + + L π(cid:96) N φ ( y + y (cid:48) ) (cid:21) (cid:18) πN φ τ (cid:48) L ( x (cid:48) − x ) (cid:12)(cid:12)(cid:12) τ (cid:48) (cid:19) ϑ (cid:20) b (cid:48) m + N φ L ( x + x (cid:48) ) (cid:21) (cid:18) π ( y − y (cid:48) ) τ L sin θ (cid:12)(cid:12)(cid:12) τ (cid:19)(cid:27) . (14)Then it is clear that the arguments of the ϑ functionsdepend on the coordinate differences only, and the dis-placement of the center of mass can be incorporated inthe characteristics as b m → b m + N φ L d x , a m → a m + L π(cid:96) N φ d y . (15)These in turn correspond to fluxes [30, 31], and the shiftof the center of mass corresponds to a change in the twistangles according to Eq. (13): φ → φ − πN φ L d x , φ → φ + L (cid:96) d y . (16)Thus the twisted boundary conditions in Eq. (6), and,consequently, | ρ PBC ( r , r (cid:48) ; β ) | , are invariant only if d = (cid:18) L N φ n , π(cid:96) L n (cid:19) (17)for integral n and n .In the β → ρ ( r , r (cid:48) ; β ) → δ ( r − r (cid:48) ), and this holds for the den-sity matrix appropriate for open boundary conditions inEq. (10). Using Eq. (9) and the identities of the tradi-tionally defined Jacobi elliptic functions [29] ϑ , ( z | τ ) = (cid:114) iτ ∞ (cid:88) n = −∞ ( ± n exp (cid:18) − iπτ (cid:16) n + zπ (cid:17) (cid:19) one can check that ρ PBC ( r , r (cid:48) ; β →
0) = (cid:88) k ,k e ik φ + ik φ − ixk L θ(cid:96) × δ ( x − x (cid:48) − k L − k L cos θ ) × δ ( y − y (cid:48) − k L sin θ ) , (18) which complies with the discrete magnetic translationsymmetries t r ( n L + m L ) ρ PBC ( r , r (cid:48) ; β ) = e i ( nφ + mφ ) ρ PBC ( r , r (cid:48) ; β ) ,t ∗ r (cid:48) ( n L + m L ) ρ PBC ( r , r (cid:48) ; β ) = e − i ( nφ + mφ ) ρ PBC ( r , r (cid:48) ; β ) , (19)which hold for any β .In the low-temperature limit, β → ∞ ( u → ρ PBC ( r , r (cid:48) ; β ) simplifies significantly.Notice that both for open and periodic boundary con-ditions, the value of the density matrix goes to zero atany fixed coordinates r and r (cid:48) . This is an artifact of thezero-point energy (cid:126) ω c , and it does not appear in averagesas they involve normalization by the partition function Z ( β ) = (cid:80) ∞ n =0 u n +1 / = √ u − u . We study the analyticstructure in the low-temperature limit by factoring outthe nonzero factor ρ open ( r , r (cid:48) ; β ) for convenience:lim β →∞ ρ PBC ( r , r (cid:48) ; β ) ρ open ( r , r (cid:48) ; β ) = f ∞ ( z, z (cid:48) ) , (20)where f ∞ ( z, z (cid:48) ) = 1 N φ N φ − (cid:88) m =0 (cid:26) ϑ (cid:20) a m (cid:21) (cid:18) iL (cid:96) N φ (cid:0) z (cid:48)∗ − z (cid:1) (cid:12)(cid:12)(cid:12) ˜ τ (cid:19) × ϑ (cid:20) b (cid:48) m (cid:21) (cid:18) N φ πL (cid:0) z + z (cid:48)∗ (cid:1) | ˜ τ (cid:19) ++( − k ϑ (cid:20) a m + (cid:21) (cid:18) iL (cid:96) N φ (cid:0) z (cid:48)∗ − z (cid:1) (cid:12)(cid:12)(cid:12) ˜ τ (cid:19) × ϑ (cid:20) b (cid:48) m (cid:21) (cid:18) N φ πL (cid:0) z + z (cid:48)∗ (cid:1) | ˜ τ (cid:19)(cid:27) , (21)where ˜ τ = iπ (cid:16) L (cid:96)N φ (cid:17) and ˜ τ = iπ (cid:16) (cid:96)N φ L (cid:17) . f ∞ ( z, z (cid:48) )is holomorphic in z , and antiholomorphic in z (cid:48) , on theentire complex plane. Fixing z (cid:48) , the zeros of f ∞ ( z, z (cid:48) )can be counted by the argument principle of complexcalculus. Consider the quadruple domain Q with corners z (cid:48) + L ( ± ± τ ); cf. Fig. 1(b). We have (cid:73) ∂Q ddz ln ( f ∞ ( z, z (cid:48) )) dz = − πiN φ , (22)which, exploiting the periodicities in Eq. (19) and thefact that ρ open ( r , r (cid:48) ; β ) is nonzero, implies that the ther-mal propagator ρ PBC ( r , r (cid:48) ; β → ∞ ) has N φ zeros in theprincipal domain in Fig. 1(a). At nonzero temperature,the analytic structure of ρ PBC ( r , r (cid:48) ; β ) is not simple. Nev-ertheless, we have found numerically that the number ofzeros in the principal domain is the same at any finite β , and the zeros very quickly reach their final location.See Fig. 3 for illustration. If N φ is odd, there are zerosthat do not move at all. For φ = φ = 0, in particular,one of them is located in the corners of the principal do-main (which are identical by periodicity). Fig. 4 shows -6-4-2 0 2 4 6 -6 -4 -2 0 2 4 6 y / l x/l 0.0001 0.001 0.01 -2.4-2.2-2-1.8-1.6-1.4-1.2-1 -3 -2.8 -2.6 -2.4 -2.2 -2(a) y / l x/l β -8-6-4-2 0 2 4 6 8 -8-6-4-2 0 2 4 6 8 y / l x/l 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 -3-2.5-2-1.5-1-0.5-5.5 -5 -4.5 -4 -3.5 -3 -2.5(b) y / l x/l β FIG. 3: The low-temperature limit β → ∞ of the thermaldensity matrix. As no change is discernible beyond β = 100,the density plot has been generated using this value. (a) N φ =4 particles, θ ≈ ◦ , | L | / | L | = 1 . φ = φ = 0. In thezoomed area we show how the zeros move to their asymptoticposition as a function of inverse temperature β . (b) The samefor N φ = 5, θ ≈ ◦ , | L | / | L | = 1 . φ = φ = 0. Notethat one of the zeros is fixed at the corner of the principalregion, which is the generic behavior when N φ is odd. the structure of zeros for different geometries. Multiplezeros occur in regular cases, as for the square unit cell inpanel (b).In Fig. 5 we show the motion of the zeros of the ther-mal density matrix as we tune the twist angles. Qual-itatively, the motion of the zeros shows an interesting -6-4-2 0 2 4 6 -4 -2 0 2 4(a) y / l x/l 1e-050.00010.0010.01 -6-4-2 0 2 4 6 -6 -4 -2 0 2 4 6(b) y / l x/l 1e-050.00010.0010.01 -8-6-4-2 0 2 4 6 8 -6 -4 -2 0 2 4 6(c) y / l x/l 1e-050.00010.0010.01 -8-6-4-2 0 2 4 6 8 -8 -6 -4 -2 0 2 4 6 8(d) y / l x/l 1e-050.00010.0010.01 FIG. 4: The structure of zeros of the thermal density matrixfor β (cid:126) ω c = 200, where the picture is stationary for differentgeometries and flux quanta. We show | ρ PBC ( r , r (cid:48) ; β ) | , thezeros are the darkest spots. We set φ = φ = 0 and fix r (cid:48) atthe origin. (a) Generic torus with N φ = 6, L /L = 1 .
13 and θ ≈ ◦ ; (b) Square principal domain ( θ = 90 ◦ , L /L = 1)with N φ = 7; (c) Generic torus with N φ = 11, L /L =1 .
19 and θ ≈ ◦ ; (d) Hexagonal principal domain ( θ = 60 ◦ , L /L = 1) with N φ = 12. analogy with the Hall current: tuning φ moves them inthe L direction–the direction of the electromotive forceon a charged particle induced by the change of flux–,and conversely. As a deeper explanation of the motionof the zeros is not crucial to the present work, we leavethe analysis of this issue as an open problem. III. MULTI-SLICE SAMPLING
For noninteracting particles and open boundary condi-tions, the familiar construction of multi-slice moves [14]by the bisection method [1] builds a Brownian bridge R L +1 , R L +2 , . . . , R R − between the fixed configurations R L and R R at possibly distant slices L and R = L + 2 l (mod M ) on the path. The deviation of the high-temperature density matrix used in the simulation fromthe ideal gas case can be taken into account either at eachor just the last level of this recursive procedure. At eachlevel of this recursive construction we need to know thePDF of configuration R i , which is to be inserted between R i − s and R i + s at time distances ± sτ on the path. If theideal gas density matrix ρ ( R, R (cid:48) ; β ) is real , this is simply p ( R i ) = ρ ( R i − s , R i ; sτ ) ρ ( R i , R i + s ; sτ ) ρ ( R i − s , R i + s ; 2 sτ ) ; (23) -4-3-2-1 0 1 2 3 4 -4 -3 -2 -1 0 1 2 3 4(a) y / l x/l -4-3-2-1 0 1 2 3 4 -4 -3 -2 -1 0 1 2 3 4(b) y / l x/l FIG. 5: The trajectories of the zeros of the thermal density matrix as we tune (a) the twist angle φ and (b) the twist angle φ between 0 and 2 π . We set β (cid:126) ω c = 200, r (cid:48) = 0, N φ = 2, L /L = 1 .
19 and θ ≈ ◦ . The marks at specific points on thetrajectories correspond to multiples of π/
5. The speed of the zeros is not uniform, as visible from the distance between adjacentlabeled points. the convolution property in Eq. (2) ensures that this isa normalized PDF. If we can sample p ( R i ) directly, weimplement the heat-bath rule for noninteracting particles.(In fact, with open boundary conditions and zero externalmagnetic field, p ( R i ) is a Gaussian.) On the other hand,if the free density matrix ρ ( R, R (cid:48) , τ ) is complex , pathsmust be sampled from the PDF (cid:81) Mm =1 | ρ ( R m − , R m ; τ ) | .As | ρ ( R, R (cid:48) , τ ) | does not satisfy a convolution propertyanalogous to Eq. (2), (cid:101) p ( R i ) = | ρ ( R i − s , R i ; sτ ) || ρ ( R i , R i + s ; sτ ) || ρ ( R i − s , R i + s ; 2 sτ ) | (24)is not a normalized PDF. This is not a problem for single-slice moves, but it plagues the bisection method.First consider how one could adapt multi-slice movesto periodic boundary conditions in the absence of a mag-netic field in one dimension. The single-particle densitymatrix is [1] ρ PBC0 ( x, x (cid:48) ; β ) = 1 L ϑ (cid:18) πL ( x − x (cid:48) ) (cid:12)(cid:12)(cid:12) πiλβL (cid:19) == 1 √ πλβ ∞ (cid:88) n = −∞ exp (cid:18) − ( x − x (cid:48) + nL ) λβ (cid:19) , (25)where L is the period. (The second equality involves amodular transformation of the function ϑ ( z | τ )). Opti-mal sampling could be achieved by the heat-bath rule onslice mT ∗ ( x (cid:48) m | x m − , x m +1 ) = ρ PBC0 ( x m − , x m ; τ ) ρ PBC0 ( x m , x m +1 ; τ ) ρ PBC0 ( x m − , x m +1 ; 2 τ ) . Sampling x (cid:48) m from this PDF results in moves thatare always accepted for noninteracting particles. Withstraightforward algebra, T ∗ ( x (cid:48) m | x m − , x m +1 ) == α ∞ (cid:88) k = −∞ exp (cid:18) − (( x m +1 + x m − ) / − x (cid:48) m + kL ) λτ (cid:19) ++ α ∞ (cid:88) k = −∞ exp (cid:18) − (( x m +1 + x m − + L ) / − x (cid:48) m + kL ) λτ (cid:19) , (26)where α i = 1 √ πλτ (cid:80) k (cid:48) exp (cid:16) − (( x m +1 − x m − + iL ) / k (cid:48) L ) λτ (cid:17)(cid:80) k (cid:48) exp (cid:16) − ( x m +1 − x m − + k (cid:48) L ) λτ (cid:17) .T ∗ ( x (cid:48) m | x m − , x m +1 ) has a very simple structure: the firstterm is a collection of the periodic copies of the Gaussianpeak centered at ( x m +1 + x m − ) /
2, the second term col-lects peaks at periodic copies of ( x m +1 + x m − + L ) / p = α / ( α + α ) we sample a Gaussian of variance λτ at( x m +1 + x m − ) /
2, with probability 1 − p we sample a sim-ilar Gaussian at ( x m +1 + x m − + L ) /
2. (With no loss ofgenerality we can choose any of the equivalent peaks, andmap x (cid:48) m back to the interval ( − L/ , L/ T ∗ in Eq. (26) can be applied on any level of the bisectionmethod to construct a free-particle trajectory betweentwo slices separated by imaginary time 2 l τ . With inter-actions present, the deviation of the high-temperaturedensity matrix that defines the PDF of paths from ρ PBC0 could be taken into account by a rejection step on thelast level of recursion. (For alternative approaches to pe-riodicity in zero magnetic field, see Ref. 32.)In the presence of an external magnetic field, the den-sity matrix in Eq. (9) is complex-valued. We samplepaths by the product of the magnitudes of the densitymatrices that connect subsequent slices. If we considermoving a bead z m on slice m with all other beads fixed. T ∗ ( z (cid:48) m | z m − , z m +1 ) = | ρ PBC ( z m − , z m ; τ ) || ρ PBC ( z m , z m +1 ; τ ) || ρ PBC ( z m − , z m +1 ; 2 τ ) | is not a normalized PDF, but this would not impairthe Metropolis algorithm. As in the β → | ρ PBC ( z, z (cid:48) ; τ ) | with fixed z (cid:48) tends to a system of Gaus-sian peaks centered at z (cid:48) + nL + mL τ , just like in thenonmagnetic case, we try the following. We choose the a priori sampling PDF T ( z (cid:48) m | z m − , z m +1 ) as a collectionof four Gaussian peaks centered at Z z m − ,z m +1 = ( z m − + z m +1 ) / ,Z z m − ,z m +1 = ( z m − + z m +1 + L ) / ,Z z m − ,z m +1 = ( z m − + z m +1 + L τ ) / ,Z z m − ,z m +1 = ( z m − + z m +1 + L (1 + τ )) / . (27)The height of these peaks is proportional to α i = | ρ PBC ( z m − , Z i ; τ ) || ρ PBC ( Z i , z m +1 ; τ ) || ρ PBC ( z m − , z m +1 ; 2 τ ) | , (28)for 0 ≤ i ≤
3. We choose peak i with probability p i = α i / ( (cid:80) j =0 α j ). We take into account the fact thatthe diffusive motion described by both | ρ open ( R, R (cid:48) , τ ) | and | ρ PBC ( R, R (cid:48) , τ ) | is different from the diffusion inthe absence of magnetic field. Thus the sampled Gaus-sian has variance − u u (cid:96) with u = e − (cid:126) ω c τ . Notice that − u u (cid:96) < λτ .As the heat-bath rule is not obeyed, the acceptanceprobability is less than unity even for noninteracting par-ticles in single-slice moves: A ( z m → z (cid:48) m ) = | ρ PBC ( z m − , z (cid:48) m ; τ ) || ρ PBC ( z (cid:48) m , z m +1 ; τ ) | ρ PBC ( z m − , z m ; τ ) || ρ PBC ( z m , z m +1 ; τ ) × T ( z m | z m − , z m +1 ) T ( z (cid:48) m | z m − , z m +1 ) . (29)For multi-slice moves, we proceed as follows.(i) A trial path is constructed recursively between slices L and R = L + 2 l . Midway between slices L and R ,we choose z (cid:48) ( L + R ) / from one of four Gaussian peaksat Z z L ,z R i of variance − u u (cid:96) , where u = e − (cid:126) ω c τ and τ = 2 l − τ . Then we sample z (cid:48) L +2 l − from one of fourGaussian peaks at Z z L ,z (cid:48) ( L + R ) / i and z (cid:48) R − l − from one offour Gaussian peaks at Z z (cid:48) ( L + R ) / ,z R i , all having variance − u u (cid:96) , where u = e − (cid:126) ω c τ and τ = 2 l − τ . We continueon subsequent levels, until the trial path z (cid:48) L +1 , . . . z (cid:48) R − is complete. During this construction, the ratio of the apriori sampling PDFs P = T ( z L +1 , . . . z R − | z L , z R ) T ( z (cid:48) L +1 , . . . z (cid:48) R − | z L , z R ) (30)is stored.(ii) Once the trial path is available, the ratio of thePDF of the new and the old paths is calculated, P = (cid:81) Rm = L +1 | ρ ( z (cid:48) m − , z (cid:48) m ; τ ) | (cid:81) Rm = L +1 | ρ ( z m − , z m ; τ ) | . (31)The constructed trial path is then accepted with proba-bility A ( z → z (cid:48) ) = P P . FIG. 6: Acceptance ratios for sampling the motion of a singleparticle on a rectangular torus pierced by N φ = 2 flux quanta.The inverse temperature of the system is β (cid:126) ω c = 2, and thenumber of slices ranged between M = 8 and 256. Level l means that 2 l − For testing the efficiency of the above algorithm, inFig. 6 we show the acceptance ratio for the simplest pos-sible case, namely the simulation of a single free particleon the torus. The phase was fixed to the density matrixin Eq. (9); we set β (cid:126) ω c = 2, and there are N φ = 2 fluxquanta through a rectangular torus. (For the computa-tional advantage of choosing N φ even, see Appendix B.)For N particles, the acceptance ratio is roughly raised tothe N -th power; this is the baseline that interactions areexpected to reduce further. We have checked systemat-ically that the acceptance ratio depends only weakly onthe aspect ratio or the twist angles. IV. APPLICATION: ROTATING YUKAWAGASES
We consider particles that interact by a repulsivemodified-Bessel-function interaction. The system rotatesabout the z -axis with angular velocity Ω. In the co-rotating frame it is described by the Hamiltonian H = − (cid:126) m N (cid:88) i =1 (cid:18) ∇ i − im (cid:126) Ω × r (cid:19) + (cid:15) (cid:88) i 3 2 1 3 2 1 -3 -2 -1 0 1 2 3 δ x-2-1 0 1 2 δ y (b) 3 2 1 3 2 1 -3 -2 -1 0 1 2 3 δ x-2-1 0 1 2 δ y (c) -3 -2 -1 0 1 2 3 δ x-2-1 0 1 2 δ y (d) -3 -2 -1 0 1 2 3 δ x-2-1 0 1 2 δ y -1-0.8-0.6-0.4-0.2 0 0.2 0.4 0.6 (e) -3 -2 -1 0 1 2 3 δ x-2-1 0 1 2 δ y -1.4-1.2-1-0.8-0.6-0.4-0.2 0 0.2 0.4 0.6 FIG. 7: (a-c) The pair-correlation function for N = 12 bosons at density ρa = 0 . 02 ( κ = 0 . ν = 2(i.e., N φ = 6 flux quanta piercing the torus) and interaction strength Λ = 0 . M = 32 sliceswere used, the imaginary time-step is τ = 0 . β = 114, 99,and 88 in panels (a) to (c). Panels (d) and (e) show the differences of the pair-correlation functions, g Λ=0 . − g Λ=0 . and g Λ=0 . − g Λ=0 . , respectively, as Λ is changed for systems shown in the top row. The triangular lattice of dark spots shows thedecreasing crystalline correlation as Λ is increased. The small deviations from perfect C symmetry in panels (b) and (c) canbe attributed to imperfect thermalization, and could be reduced by longer Monte Carlo runs. Taking the differences betweenpair-correlation functions in panels (d) and (e) amplifies these small errors. For Ω = 0, Yukawa bosons are known to exhibit non-monotonic behavior as a function of density: at fixed in-teraction strength Λ the system first crystallizes with in-creasing density, then at sufficiently high density it meltsagain. Due to computational limitations, we have onlybeen able to verify this for the Fermi system. Fig. 9 showsthe evolution of the first peak of the pair-correlation func-tion as the density changes at fixed β ∗ and Λ values forfermions. Apparently, crystalline order prevails only forintermediate densities, just like for bosons at zero tem-perature in the absence of rotation [16]. Determining thephase boundary will require more extensive simulations.In the β ∗ > β ∗ = 0 . β ∗ (cid:46) . β ∗ ≈ 1. Clearly, more comprehensive cal-culations in the large- β region are necessary to ascertainthat this tendency is robust. If so, it indicates the com-petition of the homogeneous integer quantum Hall liquidstate (the ground state candidate for this particular den- sity) and the density-wave ordering, which requires ther-mal excitations above the cyclotron gap that the interac-tion can organize in a crystalline order. This competitionis, of course, not expected for bosons or bolzmannons; forthe latter we have checked the monotonic evolution upto β ∗ = 1 . ρ ∗ isheld fixed. Again, we could study this only for fermionsand bolzmannons; some of the results are shown inFig. 11. (Notice that while β ∗ is kept constant, the sys-tem becomes colder on the interaction energy scale as˜ β = β ∗ ν/ (4 π Λ ρ ∗ ) with ν = N/N φ ; the ratio of the in-teraction and the magnetic length scale also changes as κ = (cid:112) πρ ∗ /ν .) We see that the system becomes morecrystalline as the number of flux quanta is decreased,which is only possible in very crude steps with N = 16,the largest system we simulated routinely. The tendencyis qualitatively the same for fermions and bolzmannons,but it is stronger for fermions. Note that the flux den-sity would localize particles on the scale of the magneticlength, which is greater than the lattice constant for κ < 1. On the other hand, it is more difficult to ob-0 (a) -3 -2 -1 0 1 2 3 δ x-3-2-1 0 1 2 3 δ y -0.6-0.4-0.2 0 0.2 0.4 0.6 0.8 (b) 3 2 1 3 2 1 -3 -2 -1 0 1 2 3 δ x-3-2-1 0 1 2 3 δ y (e) -3 -2 -1 0 1 2 3 δ x-3-2-1 0 1 2 3 δ y -0.4-0.2 0 0.2 0.4 0.6 0.8 (d) 3 2 1 3 2 1 -3 -2 -1 0 1 2 3 δ x-3-2-1 0 1 2 3 δ y (e) 3 2 1 3 2 1 -3 -2 -1 0 1 2 3 δ x-3-2-1 0 1 2 3 δ y (f) 4 3 2 1 4 3 2 1 -3 -2 -1 0 1 2 3 δ x-3-2-1 0 1 2 3 δ y (g) -3 -2 -1 0 1 2 3 δ x-3-2-1 0 1 2 3 δ y -2-1.5-1-0.5 0 0.5 1 (h) 4 3 2 1 4 3 2 1 -3 -2 -1 0 1 2 3 δ x-3-2-1 0 1 2 3 δ y (i) 0 -0.5 0 -0.5 -3 -2 -1 0 1 2 3 δ x-3-2-1 0 1 2 3 δ y -0.8-0.6-0.4-0.2 0 0.2 0.4 FIG. 8: Second row (d)-(f): the pair-correlation function for N = 16 fermions at density ρa = 0 . 02 ( κ = 0 . ν = 2 (i.e., N φ = 8 flux quanta piercing the torus) and interaction strength Λ = 0 . 035 at inverse temperature β ∗ = 0 . M = 16 slices were used, τ = 0 . 025 and ˜ β = 132; in (e) M = 16, τ = 0 . β = 114;in (f) M = 24, τ = 0 . 025 and ˜ β = 99. Panels (a) and (c) show the differences of the pair-correlation functions g β =0 . − g β =0 . and g β =0 . − g β =0 . , respectively, between colder and warmer systems shown in consecutive panels in the second row. Thetriangular lattice of bright spots shows the increasing crystalline correlation as the temperature is decreased. Second column(b), (e), and (h): the pair-correlation function as the temperature is held fixed at β ∗ = 0 . 5, but the de Boer parameter is tunedfrom Λ = 0 . 03 in panel (h) to Λ = 0 . 04 in panel (b). Panels (g) and (i) show the differences of the pair-correlation functions g Λ=0 . − g Λ=0 . and g Λ=0 . − g Λ=0 . , respectively, as Λ is tuned for systems shown in the second column. The triangularlattice of dark spots shows the decreasing crystalline correlation as Λ is increased. tain converged results for smaller flux densities, which isno doubt related to the shortening of the length scale onwhich the change of the phase of the many-body wavefunction can be considered smooth for the phase-fixingprocedure; in the limit of vanishing magnetic field, weapproach the sudden sign changes that are treated bynode fixing in time-reversal-symmetric simulations.We note that the PIMC calculations for N = 12 bosonsin Fig. 7 required about one day of thermalization andtwo days of data collection on a single Intel Xeon X5660 CPU core at 2.8 GHz, while the calculations for N = 16fermions in Fig. 8 were about half that long. With in-creasing inverse temperature the number of slices also hasto be increased; the most expensive calculation we per-formed was for β = 1 . FIG. 9: The height of the first peak of the pair-correlationfunction for N = 12 fermions at β ∗ = 0 . N φ = 6, forvarious Λ de Boer parameter values as the function of density ρa . The nonmonotonic evolution indicates that crystallineorder exists only for a limited range of densities.FIG. 10: The evolution of the first peak of the pair-correlation function for N = 16 fermions at flux N φ = 8at density ρa = 0 . 02 ( κ = 0 . cheaper. As the computing requirement of PIMC scalesas a moderate power, typically N , of the system size,and no attempt has yet been made to parallelize the code,we expect we can routinely simulate dozens of particlesusing the method we elaborated. (a) -0.4 -0.2 0 0.2 0.4 δ x-0.4-0.2 0 0.2 0.4 δ y -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 (b) -0.4 -0.2 0 0.2 0.4 δ x-0.4-0.2 0 0.2 0.4 δ y -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 FIG. 11: The difference g N φ =2 ( r ) − g N φ =6 ( r ) between pair-correlation functions at different flux densities (6 and 2 fluxquanta through the torus) at β ∗ = 0 . N = 16 and ρa =0 . 02 for fermions (a) and bolzmannons (b). The total area ofthe simulation call is scaled to unity, thus the peak locationsmay coincide. (The filling factor corresponding to N φ = 6 , ν = , 8, respectively.) The triangular lattice of bright peakscorrespond to stronger crystalline correlations at smaller fluxdensity. Bolzmannons in panel (b) are still liquid-like; a smallrotation of the hexagonally distorted rings from directionswhere crystalline structure will emerge can be attributed toimperfect thermalization. V. CONCLUSION AND OUTLOOK We have explored the feasibility of the path-integralMonte Carlo simulation of systems that do not obeytime-reversal symmetry under periodic boundary condi-tions. Technically, this requires the use of the single-particle thermal density matrix that is appropriate forthe boundary conditions in the presence of a magneticfield. We have derived several equivalent closed-form ex-2pressions for this purpose. The multi-slice sampling al-gorithm was modified for the case in which the weightof a path is determined by the magnitude of the den-sity matrix, which does not obey a convolution property.We have illustrated the use of these techniques in thesimulation of two-dimensional Yukawa systems, wheretime-reversal symmetry is broken by the Coriolis-force,as commonly done in experiments on cold atomic sys-tems. We have shown that in spite of the crudeness ofthe phase-fixing we used, the interaction-driven transi-tion between a crystalline phase and a correlated liquidcan be captured qualitatively by a PIMC simulation. Acomprehensive quantitative study of this system is dele-gated to future work. Eventually, fermions that interactby the Coulomb potential are of more fundamental inter-est. For such systems the primitive approximation to theaction is clearly not an adequate starting point. Moresophisticated approximations exist, but in their currentform they rely upon the consequences of time-reversalinvariance. The development of suitable approximationsfor the non-time-reversal-invariant case is underway andis delegated to future publications. Acknowledgments This research was supported by the National ResearchDevelopment and Innovation Office of Hungary withinthe Quantum Technology National Excellence Program(Project No. 2017-1.2.1-NKP-2017-00001), and by theHungarian Scientific Research Funds No. K105149. Weare grateful to the HPC facility at the Budapest Univer-sity of Technology and Economics. We thank P´eter L´evayand Bal´azs Het´enyi for useful discussions. C. T. was sup-ported by the Hungarian Academy of Sciences. H. G. T.acknowledges support from the “Quantum Computingand Quantum Technologies” PhD School of the Univer-sity of Basel. Appendix A: The derivation of the single-particledensity matrix1. Single-particle states on the torus In the gauge A = − By ˆx the states in the lowest Lan-dau level assume the form [36] ψ ( z ) = f ( z ) e − y (cid:96) , (A1)where f ( z ) is a holomorphic function. We seek the holo-morphic part f ( z ) of the lowest Landau level eigenstatesin terms of Jacobi elliptic functions, see Eq. (11). Thetwisted boundary conditions we impose in Eq. (6) yield N φ distinct states [37, 38], ψ m ( z ) = 1 (cid:112) (cid:96)L √ π ϑ (cid:20) a m b m (cid:21) (cid:18) πN φ zL (cid:12)(cid:12)(cid:12) N φ τ (cid:19) e − y (cid:96) , (A2) m = 0 , , . . . , ( N φ − a m , b m defined in Eq. (13).We note that ψ m ( z ), together with its higher Landaulevel descendants that follow later, is normalized for themagnetic unit cell, (cid:90) L sin θ dy (cid:90) y cot θ + L y cot θ dxψ ∗ n (cid:48) m (cid:48) ( x + iy ) ψ nm ( x + iy ) == δ nn (cid:48) δ mm (cid:48) . (A3)This particular basis corresponds to a string arrangement[36] of zeros of the holomorphic function f ( z ) in the prin-cipal domain.The orbitals in higher Landau levels are obtained bythe application of the Landau level ladder operators, ψ nm ( z ) = ( a † ) n √ n ! ψ m ( z ) , (A4)where ˆ a † = i(cid:96) √ ∂ z − iA z ) , (A5)with ∂ z = ( ∂ x − i∂ y ) and A z = ( A x − iA y ). In ourparticular gauge, ˆ a † = (cid:96) √ (cid:0) i∂ x + ∂ y − y(cid:96) (cid:1) . The degener-acy of each Landau level is N φ . Straightforward algebrayields ψ nm ( z ) = ( − n (cid:112) n n ! (cid:96)L √ π ∞ (cid:88) p = −∞ H n (cid:18) y + C p,m (cid:96) (cid:19) × exp (cid:0) iπτ N φ ( p + a m ) + 2 πi ( p + a m ) b m (cid:1) × exp (cid:32) C p,m (cid:96) + iC p,m x(cid:96) − ( y + C p,m ) (cid:96) (cid:33) , (A6)where C p,m = πN φ (cid:96) L ( p + a m ) = L ( p + a m ) (cid:61) τ . 2. The thermal density matrix If we substitute Eq. (A6) in the definition of the densitymatrix, Eq. (1), the summation over n can be performedby Mehler’s formula, and we get ρ PBC ( r , r (cid:48) ; β ) = √ u(cid:96)L √ π √ − u N φ − (cid:88) m =0 ∞ (cid:88) p,p (cid:48) = −∞ × exp (cid:0) iπτ N φ ( p + a m ) − iπτ ∗ N φ ( p (cid:48) + a m ) ++2 πi ( p + a m ) b m − πi ( p (cid:48) + a m ) b m ++ C p,m (cid:96) + iC p,m x(cid:96) + C p (cid:48) ,m (cid:96) − iC p (cid:48) ,m x (cid:48) (cid:96) − (cid:96) u − u (cid:0) ( y + C p,m ) + ( y (cid:48) + C p (cid:48) ,m ) (cid:1) ++ 2 u ( y + C p,m )( y (cid:48) + C p (cid:48) ,m )(1 − u ) (cid:96) (cid:19) . (A7)3Introducing new summation variables n = p + p (cid:48) and n = p − p (cid:48) , double-counting is avoided if n , n are ei-ther both even or both odd. This decouples the sum-mation variables in all terms except for a factor ofexp( iπN φ n n (cid:60) τ ). This can be omitted if Eq. (8) holds.As L /N φ is the separation of the guiding centers of or-bitals in the L direction, this condition simply meansthat a translation by L should be compatible with theseguiding center positions. By simple algebra and the ap-plication of ϑ functions in Eq. (11) we obtain ρ PBC ( r , r (cid:48) ; β ) = √ u(cid:96)L √ π √ − u × exp (cid:18) − (cid:96) u − u ( y + y (cid:48) ) + 2 u − u yy (cid:48) (cid:96) (cid:19) × N φ − (cid:88) m =0 (cid:26) ϑ (cid:20) a m (cid:21) ( z (cid:48) | τ (cid:48) ) ϑ (cid:20) b (cid:48) m (cid:21) ( z | τ )++( − k ϑ (cid:20) a m + (cid:21) ( z (cid:48) | τ (cid:48) ) ϑ (cid:20) b (cid:48) m (cid:21) ( z | τ ) (cid:27) , (A8)where we have used the definitions in Eq. (12), and τ (cid:48) = iπ (cid:18) (cid:96)N φ L (cid:19) − u u ,z (cid:48) = N φ πL (cid:18) x − x (cid:48) + i ( y + y (cid:48) ) 1 − u u (cid:19) . (A9)The density matrix in Eq. (A8) can be cast in a dif-ferent form by the application of a modular transforma- tion τ (cid:48) → τ = − τ (cid:48) , z (cid:48) → z = z (cid:48) τ (cid:48) in the correspond-ing ϑ functions. The result is Eq. (9). The structureof Eq. (9) is more transparent perhaps because the x -and y -components of the difference vector r − r (cid:48) and thecenter-of-mass vector r + r (cid:48) appear on the same footing inthe ϑ functions. Appendix B: Computational considerations While our first formula for the thermal density matrix,Eq. (A8), and the one we obtain by a modular transfor-mation, Eq. (9), are mathematically equivalent, they dodiffer from a computational point of view. As each ϑ function is computed as a sum of Gaussians with subse-quently shifted arguments, it is essential that those Gaus-sians should be narrow. This is ensured if the parameters( τ , τ (cid:48) , τ ) of those ϑ ’s have a large magnitude. Noticethat τ and τ are pure imaginary, andlim β →∞ | τ (cid:48) | = lim β →∞ | τ | = 2 N φ L sin θL , lim β →∞ | τ | = L N φ L sin θ , lim β → | τ | = lim β → | τ | = ∞ , lim β → | τ (cid:48) | = 0 . (B1)Hence it is advantageous to use Eq. (A8) for large β andEq. (9) for small β . Spelling out the summations implicitin the Jacobi ϑ functions, ρ PBC ( r , r (cid:48) ; β ) = 1 (cid:96)L √ π (cid:114) u − u N φ − (cid:88) m =0 (cid:40) ∞ (cid:88) n = −∞ A ( (cid:48) )0 mn ∞ (cid:88) n = −∞ B ( (cid:48) )0 mn + ( − k ∞ (cid:88) n = −∞ A ( (cid:48) ) mn ∞ (cid:88) n = −∞ B ( (cid:48) ) mn (cid:41) , (B2)where A dmn = exp (cid:40) iπτ (cid:48) (cid:18) n + a m + d + y + y (cid:48) L sin θ (cid:19) + 2 πiN φ ( n + a m + d ) x − x (cid:48) L (cid:41) (B3) A (cid:48) dmn = (cid:115) iτ (cid:48) exp (cid:40) i ( x (cid:48) − x )( y + y (cid:48) )2 (cid:96) + πiτ (cid:48) (cid:18) n + N φ x (cid:48) − xL (cid:19) + 2 πin (cid:18) y + y (cid:48) L sin θ + a m + d (cid:19)(cid:41) , (B4)and B dmn = exp (cid:40) iπτ (cid:18) n + d + y − y (cid:48) L sin θ (cid:19) + 2 πi ( n + d ) (cid:18) N φ x − x (cid:48) L + 2 b (cid:48) m (cid:19)(cid:41) , (B5) B (cid:48) dmn = (cid:114) iτ exp (cid:40) i ( y (cid:48) − y )( x + x (cid:48) )2 (cid:96) + 2 πib (cid:48) m ( y (cid:48) − y ) L sin θ + πiτ (cid:18) n − N φ x + x (cid:48) L − b (cid:48) m (cid:19) + 2 πin (cid:18) y − y (cid:48) L sin θ + d (cid:19)(cid:41) , (B6)Here, the A (cid:48) , B (cid:48) terms come from Eq. (9) and the un-primed ones are from Eq. (A8). Notice that A (cid:48) dmn (cid:54) = A dmn and B (cid:48) dmn (cid:54) = B dmn , the primed and unprimed4expressions are interchangeable only within the summa-tion over n and n , respectively. We have found it con-venient to use Eq. (B3) in the low-temperature rangetanh (cid:16) β (cid:126) ω c (cid:17) > L N φ L sin θ , and Eq. (B4) otherwise (hightemperature). For the other term, B dmn in Eq. (B5) isalmost always preferable to B (cid:48) dmn in Eq. (B6), except if N φ and θ are small and β large. Using Eq. (13), B dmn is independent of m iff (cid:60) τ is an integer, i.e., k (cid:48) = kN φ (B7)is an integer. Notice that this condition is stricter thanEq. (8). (Both conditions hold trivially for a rectangulartorus.) Then, using A (cid:48) mn d in Eq. (B4) and B mn d inEq. (B5), the summation over m can be performed. If,furthermore, N φ is even, an extremely compact formulais obtained: ρ PBC ( r , r (cid:48) ; β ) = 12 π(cid:96) √ u − u exp (cid:18) i ( x (cid:48) − x )( y + y (cid:48) )2 (cid:96) (cid:19) × ∞ (cid:88) n = −∞ exp (cid:18) − u − u (cid:96) ( x − x (cid:48) − n L ) ++ iπn (cid:18) N φ y + y (cid:48) L sin θ + φ π (cid:19)(cid:19) × ∞ (cid:88) n = −∞ exp (cid:18) − u − u (cid:96) ( y − y (cid:48) + n L sin θ ) ++ iπn (cid:18) N φ x + x (cid:48) L − φ − k (cid:48) φ π (cid:19)(cid:19) . (B8)Notice that Eq. (B8) amounts to obtaining the densitymatrix for twisted periodic boundary conditions from thecorresponding object for the infinite plain [Eq. (10)] asthe sum ∞ (cid:88) n ,n = −∞ e − in φ − in φ t r ( n L + n L ) ρ open ( r , r (cid:48) ; β ) . (B9)However, the two infinite summations in this formula donot decouple unless the condition in Eq. (B7) holds and N φ is even. Appendix C: Phase fixing As phase fixing for PIMC has already been describedin the literature [9], we just review the relevant formulasfor completeness. The thermal density matrix satisfiesBloch’s equation ∂∂β ρ ( R, R (cid:48) ; β ) = H ρ ( R, R (cid:48) ; β ) , (C1)where H = N (cid:88) i =1 λ (cid:16) ∇ i − e (cid:126) A ( r i ) (cid:17) + V ( R ) (C2) is the Hamiltonian that acts on the unprimed coordi-nates, and λ = (cid:126) m . We let ∇ ≡ ( ∇ , . . . , ∇ N ) and A ( R ) ≡ ( A ( r ) , . . . , A ( r N )). Separating the magnitudeand the phase of the density matrix as ρ ( R, R (cid:48) ; β ) = | ρ ( R, R (cid:48) ; β ) | e iϕ ( R,R (cid:48) ; β ) , (C3)Eq. (C1) maps to two coupled partial differential equa-tions ∂ | ρ | ∂β = λ ∇ | ρ | − (cid:20) V + λ (cid:16) ∇ ϕ − e (cid:126) A (cid:17) (cid:21) | ρ | , (C4) ∂ϕ∂β = λ (cid:18) ∇ ϕ + 2 ∇| ρ | · ∇ ϕ | ρ | − e (cid:126) A · ∇| ρ || ρ | − e (cid:126) ∇ · A (cid:19) , where we have suppressed the arguments ( R, R (cid:48) ; β ) for ρ and ϕ , and ( R ) for V and A , respectively. Consider somevariational many-body density matrix ρ T ( R, R (cid:48) ; β ) = | ρ T ( R, R (cid:48) ; β ) | e iϕ T ( R,R (cid:48) ; β ) . We seek the density matrix ρ ( R, R (cid:48) ; β ) under the assumption that ϕ ( R, R (cid:48) ; β ) = ϕ T ( R, R (cid:48) ; β ), i.e., with its phase fixed. Then Eq. (C4) isformally equivalent to a Bloch equation for | ρ ( R, R (cid:48) ; β ) | with effective potential ( R (cid:48) is fixed) V eff ( R ) = V ( R ) + λ (cid:16) ∇ ϕ T ( R, R (cid:48) ; β ) − e (cid:126) A ( R ) (cid:17) . (C5)Thus PIMC with phase fixing samples paths with realand nonnegative weight, using a fixed-phase dependenteffective interaction.If we know ϕ T ( R m , R m − ; β ) and its gradient ∇ R m ϕ T ( R m , R m − ; β ), we can apply the following ap-proximation. The gradient is decomposed into compo-nents parallel and perpendicular to the semiclassical pathbetween ( R m − , 0) and ( R m , τ ): G (cid:107) ( R ) = ∇ ϕ T ( R ) · R m − R m − | R m − R m − | ,G ⊥ ( R ) = (cid:113) |∇ ϕ T ( R ) | − ( G (cid:107) ( R )) . (C6)The perpendicular component is taken into account bythe primitive action. On the other hand, the evolutionof the phase is approximated by a cubic polynomial onthe semiclassical trajectory, and the contribution of theparallel component of the gradient of ϕ T is integrated onthis trajectory as in the semiclassical approximation tothe action. Technically, we assume the following quanti-ties are known: ϕ = lim τ ∗ → lim R → R m − ϕ T ( R, R m − ; τ ∗ ) = 0 ,g = lim τ ∗ → lim R → R m − ∇ R ϕ T ( R, R m − ; τ ∗ ) ,ϕ = ϕ T ( R m , R m − ; τ ) ,g = ∇ R m ϕ T ( R m , R m − ; τ ) , (C7)and g ⊥ , g (cid:107) , g ⊥ , g (cid:107) are magnitudes of the perpendicularand parallel components of g and g , respectively, in thesense of Eq. (C6). (If the phase is fixed to a single-particle5density matrix, g = − y (cid:48) ˆx /(cid:96) both for open and periodicboundary conditions. If the phase of the free Fermi orBose gas is used, cf. Eqs. (36-37), g = − (cid:80) i y (cid:48) i ˆx i /(cid:96) .)The perpendicular component is taken into account bythe primitive action: U FP,0 ( R m , R m − ; τ ) = λτ (cid:0) ( g ⊥ ) + ( g ⊥ ) (cid:1) . (C8)The next contribution is the line integral of ( G (cid:107) ) on thestraight path between R m − and R m , if ϕ T is approxi-mated by a cubic polynomial on this route. U FP,1 ( R m , R m − ; τ ) = λτ (cid:104) (cid:16) ( g (cid:107) ) + ( g (cid:107) ) ) (cid:17) − g (cid:107) g (cid:107) −− g (cid:107) + g (cid:107) )( ϕ − ϕ ) δR + 18 ( ϕ − ϕ ) δR (cid:35) . (C9)We proceed in the same way for the dot product of thephase gradient and the vector potential. A · G ⊥ con-tributes at the end points: U FP,2 ( R m , R m − ; τ ) = λ(cid:96) (cid:88) j =1 g ⊥ j × (cid:118)(cid:117)(cid:117)(cid:116) N (cid:88) i =1 y m − j,i − (cid:32) (cid:80) Ni =1 y m − j,i ( x m,i − x m − ,i ) δR (cid:33) , (C10) and for A · G (cid:107) we again use the semiclassical action withthe cubic approximation for ϕ T : U FP,3 ( R m , R m − ; τ ) = 2 λ(cid:96) N (cid:88) i =1 ( x m, − x m − ,i ) × (cid:34)(cid:32) c δR c g (cid:107) δR (cid:33) y m − ,i ++ (cid:32) c δR c g (cid:107) δR (cid:33) y m,i (cid:35) , (C11)where c = ϕ − ϕ ) δR − g (cid:107) + g (cid:107) δR and c = g (cid:107) + g (cid:107) δR − ϕ − ϕ ) δR .Finally, the semiclassical contribution of the A term is U FP,4 ( R m , R m − ; τ ) = λ (cid:96) × N (cid:88) i =1 (cid:0) y m − ,i + y m,i + y m − ,i y m,i (cid:1) . (C12)The total contribution is the sum of Eqs. (C8) to (C12). [1] D. M. Ceperley, Rev. Mod. Phys. , 279 (1995),URL http://link.aps.org/doi/10.1103/RevModPhys.67.279 .[2] N. Metropolis, A. W. Rosenbluth, M. N. Rosen-bluth, A. H. Teller, and E. Teller, The Jour-nal of Chemical Physics , 1087 (1953),http://dx.doi.org/10.1063/1.1699114, URL http://dx.doi.org/10.1063/1.1699114 .[3] W. K. Hastings, Biometrika , 97 (1970), URL http://dx.doi.org/10.1093/biomet/57.1.97 .[4] D. M. Ceperley, Journal of Statistical Physics , 1237(1991), ISSN 1572-9613, URL https://doi.org/10.1007/BF01030009 .[5] D. M. Ceperley, in Monte Carlo and molecular dynamicsof condensed matter systems , edited by K. Binder andG. Ciccotti (Italian Physical Society, 1996).[6] G. Ortiz, D. M. Ceperley, and R. M. Martin, Phys. Rev.Lett. , 2777 (1993), URL http://link.aps.org/doi/10.1103/PhysRevLett.71.2777 .[7] F. Bolton, Phys. Rev. B , 4780 (1996), URL https://link.aps.org/doi/10.1103/PhysRevB.54.4780 .[8] W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Ra-jagopal, Rev. Mod. Phys. , 33 (2001), URL http://link.aps.org/doi/10.1103/RevModPhys.73.33 .[9] V. Akkineni, dissertation, University of Illinois atUrbana-Champaign (2008), URL http://hdl.handle.net/2142/80567 . [10] V. Melik-Alaverdian and N. E. Bonesteel, Phys. Rev. B , R17032 (1995), URL http://link.aps.org/doi/10.1103/PhysRevB.52.R17032 .[11] V. Melik-Alaverdian, N. E. Bonesteel, and G. Ortiz,Phys. Rev. Lett. , 5286 (1997), URL http://link.aps.org/doi/10.1103/PhysRevLett.79.5286 .[12] C. Lin, F. H. Zong, and D. M. Ceperley, Phys. Rev. E , 016702 (2001), URL http://link.aps.org/doi/10.1103/PhysRevE.64.016702 .[13] J. Shao, E.-A. Kim, F. D. M. Haldane, and E. H. Rezayi,Phys. Rev. Lett. , 206402 (2015), URL http://link.aps.org/doi/10.1103/PhysRevLett.114.206402 .[14] E. L. Pollock and D. M. Ceperley, Phys. Rev. B , 2555 (1984), URL https://link.aps.org/doi/10.1103/PhysRevB.30.2555 .[15] D. R. Nelson and H. S. Seung, Phys. Rev. B , 9153 (1989), URL https://link.aps.org/doi/10.1103/PhysRevB.39.9153 .[16] W. R. Magro and D. M. Ceperley, Phys. Rev. B ,411 (1993), URL https://link.aps.org/doi/10.1103/PhysRevB.48.411 .[17] H. Nordborg and G. Blatter, Phys. Rev. Lett. , 1925 (1997), URL https://link.aps.org/doi/10.1103/PhysRevLett.79.1925 .[18] D. S. Petrov, G. E. Astrakharchik, D. J. Papoular,C. Salomon, and G. V. Shlyapnikov, Phys. Rev. Lett. , 130407 (2007), URL https://link.aps.org/doi/ .[19] K. W. Madison, F. Chevy, W. Wohlleben, and J. Dal-ibard, Phys. Rev. Lett. , 806 (2000).[20] J. R. Abo-Shaeer, C. Raman, J. M. Vogels, and W. Ket-terle, Science , 476 (2001).[21] E. Hodby, G. Hechenblaikner, S. A. Hopkins, O. M.Marag`o, and C. J. Foot, Phys. Rev. Lett. ,010405 (2001), URL https://link.aps.org/doi/10.1103/PhysRevLett.88.010405 .[22] P. C. Haljan, I. Coddington, P. Engels, and E. A. Cornell,Phys. Rev. Lett. , 210403 (2001).[23] D. M. Ceperley, Phys. Rev. Lett. , 331 (1992), URL https://link.aps.org/doi/10.1103/PhysRevLett.69.331 .[24] C. Pierleoni, D. M. Ceperley, B. Bernu, and W. R. Magro,Phys. Rev. Lett. , 2145 (1994), URL https://link.aps.org/doi/10.1103/PhysRevLett.73.2145 .[25] W. R. Magro, D. M. Ceperley, C. Pierleoni, and B. Bernu,Phys. Rev. Lett. , 1240 (1996), URL https://link.aps.org/doi/10.1103/PhysRevLett.76.1240 .[26] J. Zak, Phys. Rev. , A1602 (1964), URL http://link.aps.org/doi/10.1103/PhysRev.134.A1602 .[27] Just like its time-reversal symmetric counterpart, ρ PBC ( r , r (cid:48) ; β ) describes diffusive motion in imaginarytime (inverse temperature), but it also has a gauge-dependent phase factor. Compared to the case withoutthe external magnetic field, the diffusion process de-scribed by Eq. (10) is slower; for large imaginary time β the width of the Gaussian | ρ open ( r , r (cid:48) ; β ) | tends to afinite value determined by the magnetic length. In theother limit, β (cid:126) ω c (cid:28) Tata Lectures on Theta I , Mod-ern Birkh¨auser Classics (Springer, 1987), ISBN9780817645779, URL .[29] Traditionally defined Jacobi elliptic functions correspondto ϑ (cid:20) (cid:21) ( z | τ ) = − ϑ ( z | τ ), ϑ (cid:20) (cid:21) ( z | τ ) = ϑ ( z | τ ), ϑ (cid:20) (cid:21) ( z | τ ) = ϑ ( z | τ ), and ϑ (cid:20) (cid:21) ( z | τ ) = ϑ ( z | τ ). Noticethat sometimes a factor π is factored out of the argument,i.e., ϑ (cid:20) ab (cid:21) (cid:0) zπ (cid:12)(cid:12) τ (cid:1) is defined.[30] Y. Aharonov and D. Bohm, Phys. Rev. ,485 (1959), URL http://link.aps.org/doi/10.1103/PhysRev.115.485 .[31] N. Byers and C. N. Yang, Phys. Rev. Lett. ,46 (1961), URL http://link.aps.org/doi/10.1103/PhysRevLett.7.46 .[32] J. Cao, Phys. Rev. E , 882 (1994), URL https://link.aps.org/doi/10.1103/PhysRevE.49.882 .[33] N. K. Wilkin, J. M. F. Gunn, and R. A. Smith, Phys.Rev. Lett. , 2265 (1998), URL http://link.aps.org/doi/10.1103/PhysRevLett.80.2265 .[34] N. R. Cooper, Advances in Physics , 539 (2008).[35] Due to the computational complexity of the permanent,phase-fixing by the free Bose gas is actually more difficultthan by the free Fermi gas. In this case study we willrestrict our work to relatively small systems.[36] F. D. M. Haldane and E. H. Rezayi, Phys. Rev. B ,2529 (1985), URL http://link.aps.org/doi/10.1103/PhysRevB.31.2529 .[37] P. L´evay, Journal of Mathematical Physics , 2792(1995), URL http://scitation.aip.org/content/aip/journal/jmp/36/6/10.1063/1.531066 .[38] N. Read and E. Rezayi, Phys. Rev. B , 16864 (1996),URL http://link.aps.org/doi/10.1103/PhysRevB.54.16864http://link.aps.org/doi/10.1103/PhysRevB.54.16864