Path-integral Monte Carlo simulation of time-reversal noninvariant bulk systems with a case study of rotating Yukawa gases

TL;DR

This paper develops a path-integral Monte Carlo method for simulating bulk systems without time-reversal symmetry, applying it to study vortex melting in rotating Yukawa gases relevant to ultracold atomic systems.

Contribution

It introduces a phase-fixed Monte Carlo approach with closed-form density matrix expressions and modified sampling techniques for non-time-reversal-invariant systems.

Findings

Successfully simulated vortex melting in rotating Yukawa gases.

Provided analytical expressions for thermal density matrices on a torus.

Demonstrated applicability to ultracold Fermi-Fermi mixtures under rotation.

Abstract

We elaborate on the methodology to simulate bulk systems in the absence of time-reversal symmetry by the phase-fixed path-integral Monte Carlo method under (possibly twisted) periodic boundary conditions. Such systems include two-dimensional electrons in the quantum Hall regime and rotating ultracold Bose and Fermi gases; time-reversal symmetry is broken by an external magnetic field and the Coriolis force, respectively. We provide closed-form expressions in terms of Jacobi elliptic functions for the thermal density matrix (or the Euclidean propagator) of a single particle on a flat torus under very general conditions. We then modify the multi-slice sampling method in order to sample paths by the magnitude of the complex-valued thermal density matrix. Finally, we demonstrate that these inventions let us study the vortex melting process of a two-dimensional Yukawa gas in terms of the de Boer interaction strength parameter, temperature, and rotation (Coriolis force). The bosonic case is relevant to ultracold Fermi-Fermi mixtures of widely different masses under rotation.

Equations

$$\rho (R,R^{\prime };\beta )=\sum _n{e}^{-\beta {ϵ}_n}{\Psi }_n(R){\Psi }_n^{\ast }(R^{\prime }),$$

$$\rho (R,R^{\prime };{\beta }_1+{\beta }_2)=\int d{R}^{\prime \prime }\rho (R,R^{\prime \prime };{\beta }_1)\rho (R^{\prime \prime },R^{\prime };{\beta }_2),$$

$$\begin{gathered} \rho (R,R^{\prime };\beta )=\int d{R}_1\cdots \int d{R}_{M-1}\rho (R,R_1;\tau ) \\ \times \rho (R_1,R_2;\tau )\dots \rho (R_{M-1},R^{\prime };\tau ). \end{gathered}$$

$$\tau =\frac{L_2}{L_1}{e}^{i\theta },$$

$$t(\mathbf{L})=\exp \left(\frac{i}{\hslash }\mathbf{L}\cdot \mathbf{p}-i\frac{\hat{\mathbf{z}}\cdot (\mathbf{L}\times \mathbf{r})}{{\ell }^2}\right),$$

$$t({\mathbf{L}}_{1,2})\psi (\mathbf{r})=e^{i{\varphi }_{1,2}}\psi (\mathbf{r}).$$

$$N_{\varphi }=\frac{\left|{\mathbf{L}}_1\times {\mathbf{L}}_2\right|}{2\pi {\ell }^2}=\frac{L_1{L}_2\sin \theta }{2\pi {\ell }^2}.$$

$$L_2\cos \theta =\frac{k{L}_1}{N_{\varphi }},$$

$$\rho^{\text{PBC}}(\mathbf{r},\mathbf{r}^{\prime};\beta)=\frac{1}{N_{\phi}}\rho^{\text{open}}(\mathbf{r},\mathbf{r}^{\prime};\beta)\\ \times\sum_{m=0}^{N_{\phi}-1}\left\{\vartheta\begin{bmatrix}0\\ a_{m}\end{bmatrix}\left(z_{1}\Big{|}\tau_{1}\right)\vartheta\begin{bmatrix}0\\ 2b_{m}^{\prime}\end{bmatrix}(z_{2}|\tau_{2})+\right.\\ \left.+(-1)^{k}\vartheta\begin{bmatrix}0\\ a_{m}+\frac{1}{2}\end{bmatrix}\left(z_{1}\Big{|}\tau_{1}\right)\vartheta\begin{bmatrix}\frac{1}{2}\\ 2b_{m}^{\prime}\end{bmatrix}(z_{2}|\tau_{2})\right\},$$

$$\begin{gathered} {\rho }^{\text{open}}(\mathbf{r},{\mathbf{r}}^{\prime };\beta )=\frac{1}{2\pi {\ell }^2}\frac{\sqrt{u}}{1-u} \\ \times \exp \left(-\frac{1+u}{1-u}\frac{{\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|}^2}{4{\ell }^2}+\frac{i(x^{\prime }-x)(y+y^{\prime })}{2{\ell }^2}\right), \end{gathered}$$

$$\vartheta \left[\begin{array}{c}a\\ b\end{array}\right](z|\tau )=\sum _n{e}^{i\pi \tau (n+a{)}^2+2i(n+a)(z+b\pi )}.$$

$$\begin{array}{rl}{\tau }_1& =\frac{i}{\pi }{\left(\frac{L_1}{2\ell N_{\varphi }}\right)}^2\frac{1+u}{1-u},\\ z_1& =\frac{L_1}{4{\ell }^2{N}_{\varphi }}\left(y+y^{\prime }+i(x^{\prime }-x)\frac{1+u}{1-u}\right),\\ {\tau }_2& =i\pi {\left(\frac{2\ell N_{\varphi }}{L_1}\right)}^2\frac{1+u}{1-u},\\ z_2& =\frac{N_{\varphi }\pi }{L_1}\left(x+x^{\prime }+i(y-y^{\prime })\frac{1+u}{1-u}\right);\end{array}$$

$$\begin{array}{rl}{a}_m& =\frac{{\varphi }_1}{2\pi N_{\varphi }}+\frac{m}{N_{\varphi }},\\ b_m& =-\frac{{\varphi }_2}{2\pi }-\frac{N_{\varphi }\Re \tau }{2},\\ b_m^{\prime }& =b_m+N_{\varphi }{a}_m\Re \tau .\end{array}$$

$$\rho^{\text{PBC}}(\mathbf{r},\mathbf{r}^{\prime};\beta)=\frac{1}{N_{\phi}}\rho^{\text{open}}(\mathbf{r},\mathbf{r}^{\prime};\beta)\\ \times\sum_{m=0}^{N_{\phi}-1}\left\{\vartheta\begin{bmatrix}0\\ a_{m}+\frac{L_{1}}{4\pi\ell^{2}N_{\phi}}(y+y^{\prime})\end{bmatrix}\left(\frac{\pi N_{\phi}\tau_{1}^{\prime}}{L_{1}}(x^{\prime}-x)\Big{|}\tau_{1}^{\prime}\right)\vartheta\begin{bmatrix}0\\ 2b_{m}^{\prime}+\frac{N_{\phi}}{L_{1}}(x+x^{\prime})\end{bmatrix}\left(\frac{\pi(y-y^{\prime})\tau_{2}}{2L_{2}\sin\theta}\Big{|}\tau_{2}\right)+\right.\\ \left.+(-1)^{k}\vartheta\begin{bmatrix}0\\ a_{m}+\frac{1}{2}+\frac{L_{1}}{4\pi\ell^{2}N_{\phi}}(y+y^{\prime})\end{bmatrix}\left(\frac{\pi N_{\phi}\tau_{1}^{\prime}}{L_{1}}(x^{\prime}-x)\Big{|}\tau_{1}^{\prime}\right)\vartheta\begin{bmatrix}\frac{1}{2}\\ 2b_{m}^{\prime}+\frac{N_{\phi}}{L_{1}}(x+x^{\prime})\end{bmatrix}\left(\frac{\pi(y-y^{\prime})\tau_{2}}{2L_{2}\sin\theta}\Big{|}\tau_{2}\right)\right\}.$$

$$b_m\to b_m+\frac{N_{\varphi }}{L_1}{d}_x,a_m\to a_m+\frac{L_1}{2\pi {\ell }^2{N}_{\varphi }}{d}_y.$$

$${\varphi }_2\to {\varphi }_2-\frac{2\pi N_{\varphi }}{L_1}{d}_x,{\varphi }_1\to {\varphi }_1+\frac{L_1}{{\ell }^2}{d}_y.$$

$$\mathbf{d}=\left(\frac{L_1}{N_{\varphi }}{n}_1,\frac{2\pi {\ell }^2}{L_1}{n}_2\right)$$

$${\vartheta }_{3,2}(z|\tau )=\sqrt{\frac{i}{\tau }}\sum _{n=-\infty }^{\infty }(\pm 1{)}^n\exp \left(-\frac{i\pi }{\tau }{\left(n+\frac{z}{\pi }\right)}^2\right)$$

$$\begin{gathered} {\rho }^{\text{PBC}}(\mathbf{r},{\mathbf{r}}^{\prime };\beta \to 0)=\sum _{k_1,k_2}{e}^{i{k}_1{\varphi }_1+i{k}_2{\varphi }_2-\frac{ix{k}_2{L}_2\sin \theta }{{\ell }^2}} \\ \times \delta \left(x-x^{\prime }-k_1{L}_1-k_2{L}_2\cos \theta \right) \\ \times \delta \left(y-y^{\prime }-k_2{L}_2\sin \theta \right), \end{gathered}$$

$$\begin{array}{rl}{t}_{\mathbf{r}}(n{\mathbf{L}}_1+m{\mathbf{L}}_2){\rho }^{\text{PBC}}(\mathbf{r},{\mathbf{r}}^{\prime };\beta )& =e^{i(n{\varphi }_1+m{\varphi }_2)}{\rho }^{\text{PBC}}(\mathbf{r},{\mathbf{r}}^{\prime };\beta ),\\ t_{{\mathbf{r}}^{\prime }}^{\ast }(n{\mathbf{L}}_1+m{\mathbf{L}}_2){\rho }^{\text{PBC}}(\mathbf{r},{\mathbf{r}}^{\prime };\beta )& =e^{-i(n{\varphi }_1+m{\varphi }_2)}{\rho }^{\text{PBC}}(\mathbf{r},{\mathbf{r}}^{\prime };\beta ),\end{array}$$

$$\underset{\beta \to \infty }{\lim }\frac{{\rho }^{\text{PBC}}(\mathbf{r},{\mathbf{r}}^{\prime };\beta )}{{\rho }^{\text{open}}(\mathbf{r},{\mathbf{r}}^{\prime };\beta )}=f_{\infty }(z,z^{\prime }),$$

$$f_{\infty}(z,z^{\prime})=\frac{1}{N_{\phi}}\sum_{m=0}^{N_{\phi}-1}\left\{\vartheta\begin{bmatrix}0\\ a_{m}\end{bmatrix}\left(\frac{iL_{1}}{4\ell^{2}N_{\phi}}\left({z^{\prime}}^{\ast}-z\right)\Big{|}\tilde{\tau}_{1}\right)\right.\\ \left.\times\vartheta\begin{bmatrix}0\\ 2b_{m}^{\prime}\end{bmatrix}\left(\frac{N_{\phi}\pi}{L_{1}}\left(z+{z^{\prime}}^{\ast}\right)|\tilde{\tau}_{2}\right)+\right.\\ \left.+(-1)^{k}\vartheta\begin{bmatrix}0\\ a_{m}+\frac{1}{2}\end{bmatrix}\left(\frac{iL_{1}}{4\ell^{2}N_{\phi}}\left({z^{\prime}}^{\ast}-z\right)\Big{|}\tilde{\tau}_{1}\right)\right.\\ \left.\times\vartheta\begin{bmatrix}\frac{1}{2}\\ 2b_{m}^{\prime}\end{bmatrix}\left(\frac{N_{\phi}\pi}{L_{1}}\left(z+{z^{\prime}}^{\ast}\right)|\tilde{\tau}_{2}\right)\right\},$$

$${\oint }_{\partial Q}\frac{d}{dz}\ln \left(f_{\infty }(z,z^{\prime })\right)dz=-8\pi i{N}_{\varphi },$$

$$p(R_i)=\frac{{\rho }_0(R_{i-s},R_i;s\tau ){\rho }_0(R_i,R_{i+s};s\tau )}{{\rho }_0(R_{i-s},R_{i+s};2s\tau )};$$

$$\tilde{p}(R_i)=\frac{|{\rho }_0(R_{i-s},R_i;s\tau )||{\rho }_0(R_i,R_{i+s};s\tau )|}{|{\rho }_0(R_{i-s},R_{i+s};2s\tau )|}$$

$$\rho^{\text{PBC}}_{0}(x,x^{\prime};\beta)=\frac{1}{L}\vartheta_{3}\left(\frac{\pi}{L}(x-x^{\prime})\Big{|}\frac{4\pi i\lambda\beta}{L^{2}}\right)=\\ =\frac{1}{\sqrt{4\pi\lambda\beta}}\sum_{n=-\infty}^{\infty}\exp\left(-\frac{(x-x^{\prime}+nL)^{2}}{4\lambda\beta}\right),$$

$$T^{\ast }(x_m^{\prime }|x_{m-1},x_{m+1})=\frac{{\rho }_0^{\text{PBC}}(x_{m-1},x_m;\tau ){\rho }_0^{\text{PBC}}(x_m,x_{m+1};\tau )}{{\rho }_0^{\text{PBC}}(x_{m-1},x_{m+1};2\tau )}.$$

$$\begin{gathered} T^{\ast }(x_m^{\prime }|x_{m-1},x_{m+1})= \\ ={\alpha }_0\sum _{k=-\infty }^{\infty }\exp \left(-\frac{((x_{m+1}+x_{m-1})/2-x_m^{\prime }+kL{)}^2}{2\lambda \tau }\right)+ \\ +{\alpha }_1\sum _{k=-\infty }^{\infty }\exp \left(-\frac{((x_{m+1}+x_{m-1}+L)/2-x_m^{\prime }+kL{)}^2}{2\lambda \tau }\right), \end{gathered}$$

$${\alpha }_i=\frac{1}{\sqrt{2\pi \lambda \tau }}\frac{{\sum }_{k^{\prime }}\exp \left(-\frac{((x_{m+1}-x_{m-1}+iL)/2+k^{\prime }L{)}^2}{2\lambda \tau }\right)}{{\sum }_{k^{\prime }}\exp \left(-\frac{(x_{m+1}-x_{m-1}+k^{\prime }L{)}^2}{8\lambda \tau }\right)}.$$

$$T^{\ast }(z_m^{\prime }|z_{m-1},z_{m+1})=\frac{|{\rho }^{\text{PBC}}(z_{m-1},z_m;\tau )||{\rho }^{\text{PBC}}(z_m,z_{m+1};\tau )|}{|{\rho }^{\text{PBC}}(z_{m-1},z_{m+1};2\tau )|}$$

Full text

I Introduction

The path-integral Monte Carlo (PIMC) method Ceperley (1995) lets us simulate many-body systems at finite temperature in a controlled manner. Equilibrium properties are obtained from the many-body density matrix

[TABLE]

where $R\equiv(\mathbf{r}_{1},\mathbf{r}_{2},\dots,\mathbf{r}_{N})$ collects $dN$ particle coordinates, $d$ is the dimensionality of the system, $N$ is the number of particles, $\{\Psi_{n}\}$ is a complete set of many-body eigenstates, and $\{\epsilon_{n}\}$ are the corresponding energies. The convolution identity of the density matrix,

[TABLE]

is applied iteratively to yield the imaginary-time path-integral representation

[TABLE]

Here, the time-step $\tau\equiv\beta/M$ corresponds to a much higher temperature than the system temperature. The high-temperature density matrix that connects adjacent slices, $\rho(R_{m-1},R_{m};\tau)$, can be approximated by several plausible schemes Ceperley (1995). Estimators of physical quantities are defined by integrals that involve $\rho(R,R^{\prime};\beta)$; in most cases the diagonal element $\rho(R,R;\beta)$ is sufficient. The Metropolis-Hastings Monte Carlo method Metropolis et al. (1953); Hastings (1970) is applicable to path integration if the product of high-temperature density matrices in Eq. (3) can be interpreted as a probability density function.

For time-reversal invariant bosonic systems this always holds, and PIMC is an unbiased and essentially exact method in this case. For fermions, however, the notorious sign problem arises, because the contribution of a particular path can have either sign due to the presence of nondiagonal factors $\rho(R_{m-1},R_{m};\tau)$ in the integrand of estimators. The generic means to overcome this problem, the use of restricted or constrained paths that avoid the nodal surfaces of a preconceived trial many-body density matrix Ceperley (1991, 1996), makes PIMC variational in character.

On the other hand, if time-reversal is not a symmetry of the system, either because charged particles are exposed to an external magnetic field or the system is rotated, the density matrices are complex-valued in general and hence the prescription of the nodal surfaces is insufficient. A consequent method would be to sample paths by the probability density function (PDF) $\prod_{m=1}^{M}|\rho(R_{m-1},R_{m};\tau)|$ (here we assume integration with the diagonal density matrix as the kernel of the estimator, and we define $R=R^{\prime}\equiv R_{0}=R_{M}$), and sum them up with the complex phase factor $\prod_{m=1}^{M}\frac{\rho(R_{m-1},R_{m};\tau)}{|\rho(R_{m-1},R_{m};\tau)|}$. This procedure would result in a more severe form of the sign problem: contributions with different phase factors would cancel almost completely. This issue is equally severe for bosons and fermions, and it arises even in the nonphysical case of distinguishable particles (“bolzmannons”). In analogy to the phase-fixing extension Ortiz et al. (1993); Bolton (1996) of zero-temperature methods such as diffusion quantum Monte Carlo Foulkes et al. (2001), phase fixing is an obvious route to adapt PIMC to such problems. Unlike the case of zero-temperature methods, the function whose phase needs to be fixed is the many-body density matrix in Eq. (1), not a wave function. While the fixed-phase extension of the PIMC method is often mentioned in the literature Akkineni (2008), it is hardly ever applied, in contrast to the similar extension of zero-temperature methods Ortiz et al. (1993); Melik-Alaverdian and Bonesteel (1995); Bolton (1996); Melik-Alaverdian et al. (1997).

We address several issues related to the use of PIMC in time-reversal non-invariant bulk systems. (Finite systems such as quantum dots are not our primary interest here.) First, if we want to simulate bulk systems consequently, we have to use periodic boundary conditions, possibly with twist angles that let us reduce finite-size effects such as shell effects in finite-size representations of Fermi liquids Lin et al. (2001), which have analogs in strongly correlated electron systems in magnetic fields Shao et al. (2015). One should base any PIMC simulation on the single-particle thermal density matrix (equivalently, kinetic action) that is exact under the chosen boundary conditions. We show that the free propagation of a charged particle (equivalently, the thermal density matrix) on a flat torus subjected to a perpendicular magnetic field already exhibits a rather rich structure, although these patterns lose their significance for small imaginary times or large system sizes. This result lets us define the kinetic action in a way that is compatible with the torus.

The PIMC method is applicable beyond toy models only because the sampling of paths could be made efficient by the introduction of multi-slice moves. These replace entire segments of the path Pollock and Ceperley (1984) according to the PDF $\prod_{m=1}^{M}\rho(R_{m-1},R_{m};\tau)$. If, however, the density matrix is complex-valued and the probability density of paths is determined by its magnitude, the familiar bisection method Ceperley (1995) that relies on the Lévy construction of a Brownian bridge, runs into difficulties because the convolution property in Eq. (2) is not applicable to magnitudes. We elaborate on a modification of the multi-slice move algorithm that takes the external magnetic field and the periodicity of the torus into account.

Finally, we demonstrate the use of phase-fixed PIMC for bulk systems in a case study of rotating two-dimensional Yukawa gases. Yukawa bosons arise either in type-II superconductors, where the Abrikosov vortex lines interact by a repulsive modified-Bessel-function potential $\propto K_{0}(r)$ Nelson and Seung (1989); Magro and Ceperley (1993); Nordborg and Blatter (1997), or in strongly interacting Fermi-Fermi mixtures of ultracold atoms, if the mass ratio of the two species, $M/m$, is very far from unity and the motion of both species is confined to two dimensions Petrov et al. (2007). A flux density can be introduced to cold atomic systems by rotating the gas, a technique that has been applied frequently in the past two decades Madison et al. (2000); Abo-Shaeer et al. (2001); Hodby et al. (2001); Haljan et al. (2001). In the model we consider particles that interact via a modified-Bessel-function potential $\propto K_{0}(r)$. This is a good approximation also to the inter-atomic interaction in a Fermi-Fermi mixture at sufficiently long range Petrov et al. (2007). We do not claim, however, to represent either problem faithfully: we do not include the nonuniversal short-range repulsion between Fermi-Fermi bound states, and the inclusion of additional flux density would be difficult to justify for Abrikosov vortices. We have deliberately chosen this system for computational convenience in order to demonstrate the adequacy of our methodology. On the one hand, $K_{0}(r)$ is mildly divergent at short range, thus even the simplest approximation to the high-temperature density matrix, the primitive action, is a reasonable starting point. On the other, as $K_{0}(r)$ decays exponentially at large range, the intricacies of Ewald summation can be avoided.

As a first approach, we use the density matrix of the free Bose and Fermi gases to fix the phase of the many-body density matrix. We are encouraged in this by the fact that in the case of the node fixing problem, which arises analogously for time-reversal invariant fermionic systems, significant progress was possible both for 3He Ceperley (1992) and the hydrogen plasma Pierleoni et al. (1994); Magro et al. (1996) using the nodal surfaces of either the noninteracting system or some well-tested variational ground state wave function. (The two approaches are somewhat complementary.) Simple as it is, we demonstrate that phase-fixed PIMC captures the crystallization of rotating Yukawa bosons and fermions as a function of interaction strength, flux density, and temperature. We emphasize that unlike for the diffusion Monte Carlo or Green’s Function Monte Carlo methods, no trial wave function of the proper symmetry serves as input to such a calculation; but we do choose the aspect ratio of the unit cell so that it can accommodate a finite piece of a triangular lattice.

The paper is structured as follows. In Sec. II we present the density matrix for a single particle in a magnetic field on the torus, with some mathematical details of the derivation delegated to Appendix A, and the considerations of its efficient computation to Appendix B. The adaptation of the multi-slice sampling algorithm is discussed 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 method is applied to rotating systems of two-dimensional Yukawa gases under periodic boundary conditions. In Sec. V we summarize our results and discuss further research directions. Appendix C presents the technical details of the phase-fixing methodology for PIMC.

II The thermal density matrix

We consider a flat torus pierced by a perpendicular magnetic field. Consider the parallelogram spanned by two nonparallel vectors $\mathbf{L}_{1}=(L_{1},0)$ and $\mathbf{L}_{2}=(L_{2}\cos\theta,L_{2}\sin\theta)$. A torus is obtained by identifying the opposite sides of this unit cell; cf. Fig. 1(a).

We will refer to a similar parallelogram that has the origin as its center as the principal domain.

We use the Landau gauge $\mathbf{A}=-By\mathbf{\hat{x}}$ throughout this article. Electrons are characterized by complex coordinates $z=x+iy$, and we define

[TABLE]

so that $L_{1}$ and $L_{1}\tau$ span the parallelogram on the complex plane. In the presence of a perpendicular magnetic field, magnetic translations Zak (1964) are useful:

[TABLE]

where $\mathbf{p}=\frac{\hbar}{i}\nabla-e\mathbf{A}$. In the current gauge, these act as $t(\mathbf{L})\psi(\mathbf{r})=\exp(\frac{ix\mathbf{\hat{y}}\cdot\mathbf{L}}{\ell^{2}})\psi(\mathbf{r}+\mathbf{L})$. We will require each state and the implied density matrix to obey twisted boundary conditions with twist angles $\phi_{1,2}$,

[TABLE]

The two conditions are mutually compatible only if the parallelogram is pierced by an integral number of flux quanta,

[TABLE]

Then the principal domain is also a magnetic unit cell.

If $N_{\phi}\Re\tau=k$ is an integer, i.e.,

[TABLE]

straightforward but tedious algebra yields the single-particle density matrix

[TABLE]

where we have factored out $\rho^{\text{open}}$, the density matrix for open boundary conditions:

[TABLE]

where $\ell=\sqrt{\frac{\hbar}{eB}}$ is the magnetic length, $u=e^{-\beta\hbar\omega_{c}}$, and $\omega_{c}=\frac{eB}{m}$ is the cyclotron frequency com (a). Above, we have used Jacobi elliptic functions with characteristics Mumford (1987); com (b)

[TABLE]

The arguments in Eq. (9) are defined as

[TABLE]

and the constants related to boundary conditions are

[TABLE]

The derivation of Eq. (9) is delegated to Appendix A.

The behavior of the density matrix is shown in Fig. 2 for the most general case, an oblique unit cell. For small imaginary time (high temperature) $|\rho^{\text{PBC}}(\mathbf{r},\mathbf{r}^{\prime};\beta)|$ has a small Gaussian peak around $\mathbf{r}^{\prime}$, which is fixed at the origin in the figure. This peak spreads out by diffusion as $\beta$ is increased, and eventually the Gaussians from neighboring unit cells start to overlap appreciably. However, the density matrix also has a phase due to the external magnetic field, which gives rise to an interference pattern in this time range. There is destructive interference at certain points, which effectively arrests the diffusion. [We will analyze the zeros of $\rho^{\text{PBC}}(\mathbf{r},\mathbf{r}^{\prime};\beta)$ below.] Beyond a certain value of $\beta$, the picture is essentially stationary.

We note that $|\rho^{\text{PBC}}(\mathbf{r},\mathbf{r}^{\prime};\beta)|$ is not invariant for a simultaneous displacement of both $\mathbf{r}$ and $\mathbf{r}^{\prime}$ by the same vector $\mathbf{d}$, which corresponds to choosing a shifted magnetic unit cell on the plane for compactification by periodic boundary conditions, except for special choices of $\mathbf{d}$. This is understood easily by noting that the second characteristic $b_{m}$ appears in Eq. (11) as a simple additive constant to the variable $z$, letting us rewrite Eq. (9) as

[TABLE]

Then it is clear that the arguments of the $\vartheta$ functions depend on the coordinate differences only, and the displacement of the center of mass can be incorporated in the characteristics as

[TABLE]

These in turn correspond to fluxes Aharonov and Bohm (1959); Byers and Yang (1961), and the shift of the center of mass corresponds to a change in the twist angles according to Eq. (13):

[TABLE]

Thus the twisted boundary conditions in Eq. (6), and, consequently, $|\rho^{\text{PBC}}(\mathbf{r},\mathbf{r}^{\prime};\beta)|$, are invariant only if

[TABLE]

for integral $n_{1}$ and $n_{2}$.

In the $\beta\to 0$ limit the density matrix must satisfy $\rho(\mathbf{r},\mathbf{r}^{\prime};\beta)\to\delta(\mathbf{r}-\mathbf{r}^{\prime})$, and this holds for the density matrix appropriate for open boundary conditions in Eq. (10). Using Eq. (9) and the identities of the traditionally defined Jacobi elliptic functions com (b)

[TABLE]

one can check that

[TABLE]

which complies with the discrete magnetic translation symmetries

[TABLE]

which hold for any $\beta$.

In the low-temperature limit, $\beta\to\infty$ ($u\to 0$), the analytic structure of $\rho^{\text{PBC}}(\mathbf{r},\mathbf{r}^{\prime};\beta)$ simplifies significantly. Notice that both for open and periodic boundary conditions, the value of the density matrix goes to zero at any fixed coordinates $\mathbf{r}$ and $\mathbf{r}^{\prime}$. This is an artifact of the zero-point energy $\frac{\hbar\omega_{c}}{2}$, and it does not appear in averages as they involve normalization by the partition function $Z(\beta)=\sum_{n=0}^{\infty}u^{n+1/2}=\frac{\sqrt{u}}{1-u}$. We study the analytic structure in the low-temperature limit by factoring out the nonzero factor $\rho^{\text{open}}(\mathbf{r},\mathbf{r}^{\prime};\beta)$ for convenience:

[TABLE]

where

[TABLE]

where $\tilde{\tau}_{1}=\frac{i}{\pi}\left(\frac{L_{1}}{2\ell N_{\phi}}\right)^{2}$ and $\tilde{\tau}_{2}=i\pi\left(\frac{2\ell N_{\phi}}{L_{1}}\right)^{2}$. $f_{\infty}(z,z^{\prime})$ is holomorphic in $z$, and antiholomorphic in $z^{\prime}$, on the entire complex plane. Fixing $z^{\prime}$, the zeros of $f_{\infty}(z,z^{\prime})$ can be counted by the argument principle of complex calculus. Consider the quadruple domain $Q$ with corners $z^{\prime}+L_{1}(\pm 1\pm\tau)$; cf. Fig. 1(b). We have

[TABLE]

which, exploiting the periodicities in Eq. (19) and the fact that $\rho^{\text{open}}(\mathbf{r},\mathbf{r}^{\prime};\beta)$ is nonzero, implies that the thermal propagator $\rho^{\text{PBC}}(\mathbf{r},\mathbf{r}^{\prime};\beta\to\infty)$ has $N_{\phi}$ zeros in the principal domain in Fig. 1(a). At nonzero temperature, the analytic structure of $\rho^{\text{PBC}}(\mathbf{r},\mathbf{r}^{\prime};\beta)$ is not simple. Nevertheless, we have found numerically that the number of zeros in the principal domain is the same at any finite $\beta$, and the zeros very quickly reach their final location. See Fig. 3 for illustration. If $N_{\phi}$ is odd, there are zeros that do not move at all. For $\phi_{1}=\phi_{2}=0$, in particular, one of them is located in the corners of the principal domain (which are identical by periodicity).

Fig. 4 shows the structure of zeros for different geometries. Multiple zeros occur in regular cases, as for the square unit cell in panel (b).

In Fig. 5 we show the motion of the zeros of the thermal density matrix as we tune the twist angles. Qualitatively, the motion of the zeros shows an interesting analogy with the Hall current: tuning $\phi_{1}$ moves them in the $\mathbf{L}_{2}$ direction–the direction of the electromotive force on a charged particle induced by the change of flux–, and conversely. As a deeper explanation of the motion of the zeros is not crucial to the present work, we leave the analysis of this issue as an open problem.

III Multi-slice sampling

For noninteracting particles and open boundary conditions, the familiar construction of multi-slice moves Pollock and Ceperley (1984) by the bisection method Ceperley (1995) builds a Brownian bridge $R_{L+1},R_{L+2},\dots,R_{R-1}$ 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 from the ideal gas case can be taken into account either at each or just the last level of this recursive procedure. At each level of this recursive construction we need to know the PDF of configuration $R_{i}$, which is to be inserted between $R_{i-s}$ and $R_{i+s}$ at time distances $\pm s\tau$ on the path. If the ideal gas density matrix $\rho_{0}(R,R^{\prime};\beta)$ is real, this is simply

[TABLE]

the convolution property in Eq. (2) ensures that this is a normalized PDF. If we can sample $p(R_{i})$ directly, we implement the heat-bath rule for noninteracting particles. (In fact, with open boundary conditions and zero external magnetic field, $p(R_{i})$ is a Gaussian.) On the other hand, if the free density matrix $\rho_{0}(R,R^{\prime},\tau)$ is complex, paths must be sampled from the PDF $\prod_{m=1}^{M}|\rho(R_{m-1},R_{m};\tau)|$. As $|\rho_{0}(R,R^{\prime},\tau)|$ does not satisfy a convolution property analogous to Eq. (2),

[TABLE]

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 moves to periodic boundary conditions in the absence of a magnetic field in one dimension. The single-particle density matrix is Ceperley (1995)

[TABLE]

where $L$ is the period. (The second equality involves a modular transformation of the function $\vartheta_{3}(z|\tau)$). Optimal sampling could be achieved by the heat-bath rule on slice $m$

[TABLE]

Sampling $x^{\prime}_{m}$ from this PDF results in moves that are always accepted for noninteracting particles. With straightforward algebra,

[TABLE]

where

[TABLE]

$T^{\ast}(x^{\prime}_{m}|x_{m-1},x_{m+1})$ has a very simple structure: the first term is a collection of the periodic copies of the Gaussian peak centered at $(x_{m+1}+x_{m-1})/2$, the second term collects peaks at periodic copies of $(x_{m+1}+x_{m-1}+L)/2$. This suggests a very simple algorithm: with probability $p=\alpha_{0}/(\alpha_{0}+\alpha_{1})$ we sample a Gaussian of variance $\lambda\tau$ at $(x_{m+1}+x_{m-1})/2$, with probability $1-p$ we sample a similar Gaussian at $(x_{m+1}+x_{m-1}+L)/2$. (With no loss of generality we can choose any of the equivalent peaks, and map $x^{\prime}_{m}$ back to the interval $(-L/2,L/2)$.) Further, $T^{\ast}$ in Eq. (26) can be applied on any level of the bisection method to construct a free-particle trajectory between two slices separated by imaginary time $2^{l}\tau$. With interactions present, the deviation of the high-temperature density matrix that defines the PDF of paths from $\rho^{\text{PBC}}_{0}$ could be taken into account by a rejection step on the last level of recursion. (For alternative approaches to periodicity in zero magnetic field, see Ref. Cao, 1994.)

In the presence of an external magnetic field, the density matrix in Eq. (9) is complex-valued. We sample paths by the product of the magnitudes of the density matrices that connect subsequent slices. If we consider moving a bead $z_{m}$ on slice $m$ with all other beads fixed.

[TABLE]

is not a normalized PDF, but this would not impair the Metropolis algorithm. As in the $\beta\to 0$ limit $|\rho^{\text{PBC}}(z,z^{\prime};\tau)|$ with fixed $z^{\prime}$ tends to a system of Gaussian peaks centered at $z^{\prime}+nL_{1}+mL_{1}\tau$, just like in the nonmagnetic case, we try the following. We choose the a priori sampling PDF $T(z^{\prime}_{m}|z_{m-1},z_{m+1})$ as a collection of four Gaussian peaks centered at

[TABLE]

The height of these peaks is proportional to

[TABLE]

for $0\leq i\leq 3$. We choose peak $i$ with probability $p_{i}=\alpha_{i}/(\sum_{j=0}^{3}\alpha_{j})$. We take into account the fact that the diffusive motion described by both $|\rho^{\text{open}}(R,R^{\prime},\tau)|$ and $|\rho^{\text{PBC}}(R,R^{\prime},\tau)|$ is different from the diffusion in the absence of magnetic field. Thus the sampled Gaussian has variance $\frac{1-u}{1+u}\ell^{2}$ with $u=e^{-\hbar\omega_{c}\tau}$. Notice that $\frac{1-u}{1+u}\ell^{2}<\lambda\tau$.

As the heat-bath rule is not obeyed, the acceptance probability is less than unity even for noninteracting particles in single-slice moves:

[TABLE]

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^{\prime}_{(L+R)/2}$ from one of four Gaussian peaks at $Z^{z_{L},z_{R}}_{i}$ of variance $\frac{1-u_{1}}{1+u_{1}}\ell^{2}$, where $u_{1}=e^{-\hbar\omega_{c}\tau_{1}}$ and $\tau_{1}=2^{l-1}\tau$. Then we sample $z^{\prime}_{L+2^{l-2}}$ from one of four Gaussian peaks at $Z^{z_{L},z^{\prime}_{(L+R)/2}}_{i}$ and $z^{\prime}_{R-2^{l-2}}$ from one of four Gaussian peaks at $Z^{z^{\prime}_{(L+R)/2},z_{R}}_{i}$, all having variance $\frac{1-u_{2}}{1+u_{2}}\ell^{2}$, where $u_{2}=e^{-\hbar\omega_{c}\tau_{2}}$ and $\tau_{2}=2^{l-2}\tau$. We continue on subsequent levels, until the trial path $z^{\prime}_{L+1},\dots z^{\prime}_{R-1}$ is complete. During this construction, the ratio of the a priori sampling PDFs

[TABLE]

is stored.

(ii) Once the trial path is available, the ratio of the PDF of the new and the old paths is calculated,

[TABLE]

The constructed trial path is then accepted with probability $A(z\to z^{\prime})=P_{1}P_{2}$.

For testing the efficiency of the above algorithm, in Fig. 6 we show the acceptance ratio for the simplest possible case, namely the simulation of a single free particle on the torus. The phase was fixed to the density matrix in Eq. (9); we set $\beta\hbar\omega_{c}=2$, and there are $N_{\phi}=2$ flux quanta through a rectangular torus. (For the computational advantage of choosing $N_{\phi}$ even, see Appendix B.) For $N$ particles, the acceptance ratio is roughly raised to the $N$-th power; this is the baseline that interactions are expected to reduce further. We have checked systematically that the acceptance ratio depends only weakly on the aspect ratio or the twist angles.

IV Application: rotating Yukawa gases

We consider particles that interact by a repulsive modified-Bessel-function interaction. The system rotates about the $z$-axis with angular velocity $\Omega$. In the co-rotating frame it is described by the Hamiltonian

[TABLE]

where $\epsilon$ and $a$ characterize the strength and the range of the interaction, respectively. The correspondence between $\Omega$ and the formerly defined cyclotron frequency and magnetic length scales is

[TABLE]

We consider both Bose and spinless Fermi systems.

In cold atomic experiments a confinement potential is also present, which is weakened by the centrifugal force in the co-rotating frame. We do not include these terms; we describe a homogeneous portion of the gas. As is apparent from Eq. (32), the Coriolis force couples to momenta just like a uniform magnetic field does for charged particles Wilkin et al. (1998); Cooper (2008).

For $\Omega=0$, a mathematically equivalent system arises in type-II superconductors, where the bosons correspond to Abrikosov vortex lines Nelson and Seung (1989). Both the ground state Magro and Ceperley (1993) and the finite-temperature Nordborg and Blatter (1997) phase diagram of this time-reversal invariant system have been explored by quantum Monte Carlo techniques.

There are four energy scales in the problem: the temperature $k_{\text{B}}T\equiv\beta^{-1}$, the cyclotron energy $\hbar\omega_{c}$, the interaction strength $\epsilon$, and the energy that corresponds to the interaction length scale, $\hbar^{2}/(2ma^{2})$. We introduce the dimensionless parameters

[TABLE]

where $\rho$ is the particle density and $\ell$ is the magnetic length. $\Lambda$ is the de Boer interaction strength parameter. We could also have used

[TABLE]

to turn the inverse temperature dimensionless; the two dimensionless temperature parameters are related as $\tilde{\beta}=\beta^{\ast}/(2\kappa^{2}\Lambda)$. The dimensionless density can be related to the filling factor $\nu$ of Landau levels as $\rho^{\ast}=\kappa^{2}\nu/2\pi$.

With time-reversal symmetry, the system orders in a triangular lattice for strong interaction (small $\Lambda$) Magro and Ceperley (1993); Nordborg and Blatter (1997). With this prior knowledge, we choose the aspect ratio of the rectangular simulation cell so that it can accommodate a finite piece of a triangular lattice with periodic boundary conditions. This means $\sqrt{3}/2$ for $N=4$, 12 and 16 particles, and $\sqrt{3}$ for $N=8$ particles. We emphasize that this choice is the only a priori input to our simulation. The ideal Bose and Fermi gas, respectively, that we use for phase fixing is not ideal either for a crystal or a correlated liquid.

In analogy to free-particle nodes, we fix the phase to the density matrix of the ideal gas,

[TABLE]

for fermions, and

[TABLE]

for bosons; Perm stands for the permanent. As we will see, such an ansatz is sufficiently nonrestrictive for reasonable predictions com (c). (Computationally, of course, the Fermi case is easier.) As phase-fixing for PIMC has already been discussed in the literature Akkineni (2008), we are content with summarizing the technicalities in Appendix C.

The pair-correlation function for $N=12$ bosons at $\beta^{\ast}=0.5$ is shown in Fig. 7. Qualitatively, the transition to the crystalline structure is captured. Due to computational limitations, however, we cannot simulate more than 12 bosons. The pair-correlation for a larger Fermi system is shown in Fig. 8. The qualitative behavior is similar. Notice that the small $\beta^{\ast}$ means that while temperature destroys magnetic effects, it is still small on the interaction energy scale; $\tilde{\beta}$ is on the scale of $10^{2}$. (In the absence of flux, Ref. Nordborg and Blatter (1997) finds essentially ground-state behavior at $\tilde{\beta}\approx 300$.)

It is customary to characterize the crystalline order by the Lindemann ratio $\gamma=\sqrt{\frac{1}{N}\sum_{i=1}^{N}\left\langle(\mathbf{r}_{i}-\mathbf{R}_{i})^{2}\right\rangle}/d$, where $d$ is the lattice constant and $\mathbf{R}_{i}$ is the lattice point nearest to particle $i$. In our case, however, we cannot hold the center of mass fixed during Monte Carlo, because the simultaneous shift of all beads by the same vector is not a symmetry, except for some discrete values, as discussed in Sec. II. One could locate the lattice points with reference to the instantaneous center of mass assuming the lattice is triangular with the lattice constant implied by the density. But this procedure underestimates $\gamma$. Hence, we decided to infer the qualitative behavior from the pair-correlation function instead.

By inspecting the difference of the pair-correlation functions of systems that differ only by one parameter, we have checked that in the $\beta^{\ast}<1$ range our method reproduces the tendencies known for the nonrotating system: the crystalline tendency becomes stronger with increasing $\beta^{\ast}$ at fixed $\Lambda$ and $\rho^{\ast}$, as seen in the related panels of Figs. 7, 8 and 10, and it becomes stronger when decreasing $\Lambda$ at fixed $\beta^{\ast}$ and $\rho^{\ast}$. Also, Fermi systems show stronger peaks in the pair-correlation than Bose systems at identical temperature, density, and de Boer parameter $\Lambda$. It is not possible to go beyond qualitative statements now, as neither finite-size scaling nor a $\tau\to 0$ extrapolation has been performed. With the prior knowledge that the melting transition is first-order, it will be necessary to perform simulations with the particle density as a dynamical variable Nordborg and Blatter (1997). As our goal is to demonstrate the applicability of PIMC to bulk systems in the absence of time-reversal symmetry, and not an in-depth analysis of the Yukawa system, we delegate such a quantitative analysis to future work.

For $\Omega=0$, Yukawa bosons are known to exhibit nonmonotonic behavior as a function of density: at fixed interaction strength $\Lambda$ the system first crystallizes with increasing density, then at sufficiently high density it melts again. Due to computational limitations, we have only been able to verify this for the Fermi system. Fig. 9 shows the evolution of the first peak of the pair-correlation function as the density changes at fixed $\beta^{\ast}$ and $\Lambda$ values for fermions. Apparently, crystalline order prevails only for intermediate densities, just like for bosons at zero temperature in the absence of rotation Magro and Ceperley (1993). Determining the phase boundary will require more extensive simulations.

In the $\beta^{\ast}>1$ range the strength of the crystalline correlations apparently starts to weaken as a function of the inverse temperature for fermions. Such an evolution is shown in Fig. 10 for various de Boer interaction parameters $\Lambda$ as the temperature is tuned from $\beta^{\ast}=0.1$ to 1.2. The pair correlation becomes more crystalline in the $\beta^{\ast}\lesssim 0.6$ range, then stagnates, and seems to weaken again above $\beta^{\ast}\approx 1$. Clearly, more comprehensive calculations in the large-$\beta$ region are necessary to ascertain that this tendency is robust. If so, it indicates the competition of the homogeneous integer quantum Hall liquid state (the ground state candidate for this particular density) and the density-wave ordering, which requires thermal excitations above the cyclotron gap that the interaction can organize in a crystalline order. This competition is, of course, not expected for bosons or bolzmannons; for the latter we have checked the monotonic evolution up to $\beta^{\ast}=1.8$.

It is also interesting to review the evolution of the pair-correlation as a function of flux density (magnetic field or Coriolis-force) when the particle density $\rho^{\ast}$ is held fixed. Again, we could study this only for fermions and bolzmannons; some of the results are shown in Fig. 11. (Notice that while $\beta^{\ast}$ is kept constant, the system becomes colder on the interaction energy scale as $\tilde{\beta}=\beta^{\ast}\nu/(4\pi\Lambda\rho^{\ast})$ with $\nu=N/N_{\phi}$; the ratio of the interaction and the magnetic length scale also changes as $\kappa=\sqrt{2\pi\rho^{\ast}/\nu}$.) We see that the system becomes more crystalline 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 tendency is qualitatively the same for fermions and bolzmannons, but it is stronger for fermions. Note that the flux density would localize particles on the scale of the magnetic length, which is greater than the lattice constant for $\kappa<1$. On the other hand, it is more difficult to obtain converged results for smaller flux densities, which is no doubt related to the shortening of the length scale on which the change of the phase of the many-body wave function can be considered smooth for the phase-fixing procedure; in the limit of vanishing magnetic field, we approach the sudden sign changes that are treated by node fixing in time-reversal-symmetric simulations.

We note that the PIMC calculations for $N=12$ bosons in Fig. 7 required about one day of thermalization and two days of data collection on a single Intel Xeon X5660 CPU core at 2.8 GHz, while the calculations for $N=16$ fermions in Fig. 8 were about half that long. With increasing inverse temperature the number of slices also has to be increased; the most expensive calculation we performed was for $\beta=1.1$ in Fig. 10, with three days of thermalization and eleven days of data collection. The number of flux quanta hardly affects the resources needed: each of the calculations compared in Fig. 11(a) required about three plus six days; the calculations for distinguishable particles in Fig. 11(b) were about a factor of 3 cheaper. As the computing requirement of PIMC scales as a moderate power, typically $N^{3}$, of the system size, and no attempt has yet been made to parallelize the code, we expect we can routinely simulate dozens of particles using the method we elaborated.

V Conclusion and Outlook

We have explored the feasibility of the path-integral Monte Carlo simulation of systems that do not obey time-reversal symmetry under periodic boundary conditions. Technically, this requires the use of the single-particle thermal density matrix that is appropriate for the boundary conditions in the presence of a magnetic field. We have derived several equivalent closed-form expressions for this purpose. The multi-slice sampling algorithm was modified for the case in which the weight of a path is determined by the magnitude of the density matrix, which does not obey a convolution property. We have illustrated the use of these techniques in the simulation of two-dimensional Yukawa systems, where time-reversal symmetry is broken by the Coriolis-force, as commonly done in experiments on cold atomic systems. We have shown that in spite of the crudeness of the phase-fixing we used, the interaction-driven transition between a crystalline phase and a correlated liquid can be captured qualitatively by a PIMC simulation. A comprehensive quantitative study of this system is delegated to future work. Eventually, fermions that interact by the Coulomb potential are of more fundamental interest. For such systems the primitive approximation to the action is clearly not an adequate starting point. More sophisticated approximations exist, but in their current form they rely upon the consequences of time-reversal invariance. The development of suitable approximations for the non-time-reversal-invariant case is underway and is delegated to future publications.

Acknowledgements.

This research was supported by the National Research Development and Innovation Office of Hungary within the Quantum Technology National Excellence Program (Project No. 2017-1.2.1-NKP-2017-00001), and by the Hungarian Scientific Research Funds No. K105149. We are grateful to the HPC facility at the Budapest University of Technology and Economics. We thank Péter Lévay and Balázs Hetényi for useful discussions. C. T. was supported by the Hungarian Academy of Sciences. H. G. T. acknowledges support from the “Quantum Computing and Quantum Technologies” PhD School of the University of Basel.

Appendix A The derivation of the single-particle density matrix

A.1 Single-particle states on the torus

In the gauge $\mathbf{A}=-By\mathbf{\hat{x}}$ the states in the lowest Landau level assume the form Haldane and Rezayi (1985)

[TABLE]

where $f(z)$ is a holomorphic function. We seek the holomorphic part $f(z)$ of the lowest Landau level eigenstates in terms of Jacobi elliptic functions, see Eq. (11). The twisted boundary conditions we impose in Eq. (6) yield $N_{\phi}$ distinct states Lévay (1995); Read and Rezayi (1996),

[TABLE]

$m=0,1,\dots,(N_{\phi}-1)$, and $a_{m},b_{m}$ defined in Eq. (13). We note that $\psi_{0m}(z)$, together with its higher Landau level descendants that follow later, is normalized for the magnetic unit cell,

[TABLE]

This particular basis corresponds to a string arrangement Haldane and Rezayi (1985) of zeros of the holomorphic function $f(z)$ in the principal domain.

The orbitals in higher Landau levels are obtained by the application of the Landau level ladder operators,

[TABLE]

where

[TABLE]

with $\partial_{z}=\frac{1}{2}(\partial_{x}-i\partial_{y})$ and $A_{z}=\frac{1}{2}(A_{x}-iA_{y})$. In our particular gauge, $\hat{a}^{\dagger}=\frac{\ell}{\sqrt{2}}\left(i\partial_{x}+\partial_{y}-\frac{y}{\ell^{2}}\right)$. The degeneracy of each Landau level is $N_{\phi}$. Straightforward algebra yields

[TABLE]

where $C_{p,m}=\frac{2\pi N_{\phi}\ell^{2}}{L_{1}}(p+a_{m})=L_{1}(p+a_{m})\Im\tau$.

A.2 The thermal density matrix

If we substitute Eq. (43) in the definition of the density matrix, Eq. (1), the summation over $n$ can be performed by Mehler’s formula, and we get

[TABLE]

Introducing new summation variables $n_{1}=p+p^{\prime}$ and $n_{2}=p-p^{\prime}$, double-counting is avoided if $n_{1},n_{2}$ are either both even or both odd. This decouples the summation variables in all terms except for a factor of $\exp(i\pi N_{\phi}n_{1}n_{2}\Re\tau)$. This can be omitted if Eq. (8) holds. As $L_{1}/N_{\phi}$ is the separation of the guiding centers of orbitals in the $\mathbf{L}_{1}$ direction, this condition simply means that a translation by $\mathbf{L}_{2}$ should be compatible with these guiding center positions. By simple algebra and the application of $\vartheta$ functions in Eq. (11) we obtain

[TABLE]

where we have used the definitions in Eq. (12), and

[TABLE]

The density matrix in Eq. (45) can be cast in a different form by the application of a modular transformation $\tau^{\prime}_{1}\to\tau_{1}=-\frac{1}{\tau^{\prime}_{1}}$, $z^{\prime}_{1}\to z_{1}=\frac{z^{\prime}_{1}}{\tau^{\prime}_{1}}$ in the corresponding $\vartheta$ functions. The result is Eq. (9). The structure of Eq. (9) is more transparent perhaps because the $x$- and $y$-components of the difference vector $\mathbf{r}-\mathbf{r}^{\prime}$ and the center-of-mass vector $\frac{\mathbf{r}+\mathbf{r}^{\prime}}{2}$ appear on the same footing in the $\vartheta$ functions.

Appendix B Computational considerations

While our first formula for the thermal density matrix, Eq. (45), and the one we obtain by a modular transformation, Eq. (9), are mathematically equivalent, they do differ from a computational point of view. As each $\vartheta$ function is computed as a sum of Gaussians with subsequently shifted arguments, it is essential that those Gaussians should be narrow. This is ensured if the parameters ($\tau_{1}$, $\tau^{\prime}_{1}$, $\tau_{2}$) of those $\vartheta$’s have a large magnitude. Notice that $\tau_{1}$ and $\tau_{2}$ are pure imaginary, and

[TABLE]

Hence it is advantageous to use Eq. (45) for large $\beta$ and Eq. (9) for small $\beta$. Spelling out the summations implicit in the Jacobi $\vartheta$ functions,

[TABLE]

where

[TABLE]

and

[TABLE]

Here, the $A^{\prime}$, $B^{\prime}$ terms come from Eq. (9) and the unprimed ones are from Eq. (45). Notice that $A^{\prime}_{dmn_{1}}\neq A_{dmn_{1}}$ and $B^{\prime}_{dmn_{1}}\neq B_{dmn_{1}}$, the primed and unprimed expressions are interchangeable only within the summation over $n_{1}$ and $n_{2}$, respectively. We have found it convenient to use Eq. (49) in the low-temperature range $\tanh\left(\frac{\beta\hbar\omega_{c}}{2}\right)>\frac{L_{1}}{2N_{\phi}L_{2}\sin\theta}$, and Eq. (50) otherwise (high temperature). For the other term, $B_{dmn_{2}}$ in Eq. (51) is almost always preferable to $B^{\prime}_{dmn_{2}}$ in Eq. (52), except if $N_{\phi}$ and $\theta$ are small and $\beta$ large. Using Eq. (13), $B_{dmn_{2}}$ is independent of $m$ iff $\Re\tau$ is an integer, i.e.,

[TABLE]

is an integer. Notice that this condition is stricter than Eq. (8). (Both conditions hold trivially for a rectangular torus.) Then, using $A^{\prime}_{mn_{1}d}$ in Eq. (50) and $B_{mn_{2}d}$ in Eq. (51), the summation over $m$ can be performed. If, furthermore, $N_{\phi}$ is even, an extremely compact formula is obtained:

[TABLE]

Notice that Eq. (54) amounts to obtaining the density matrix for twisted periodic boundary conditions from the corresponding object for the infinite plain [Eq. (10)] as the sum

[TABLE]

However, the two infinite summations in this formula do not decouple unless the condition in Eq. (53) holds and $N_{\phi}$ is even.

Appendix C Phase fixing

As phase fixing for PIMC has already been described in the literature Akkineni (2008), we just review the relevant formulas for completeness. The thermal density matrix satisfies Bloch’s equation

[TABLE]

where

[TABLE]

is the Hamiltonian that acts on the unprimed coordinates, and $\lambda=\frac{\hbar^{2}}{2m}$. We let $\nabla\equiv(\nabla_{1},\dots,\nabla_{N})$ and $A(R)\equiv(\mathbf{A}(\mathbf{r}_{1}),\dots,\mathbf{A}(\mathbf{r}_{N}))$. Separating the magnitude and the phase of the density matrix as

[TABLE]

Eq. (56) maps to two coupled partial differential equations

[TABLE]

where we have suppressed the arguments $(R,R^{\prime};\beta)$ for $\rho$ and $\varphi$, and $(R)$ for $V$ and $A$, respectively. Consider some variational many-body density matrix $\rho_{T}(R,R^{\prime};\beta)=|\rho_{T}(R,R^{\prime};\beta)|e^{i\varphi_{T}(R,R^{\prime};\beta)}$. We seek the density matrix $\rho(R,R^{\prime};\beta)$ under the assumption that $\varphi(R,R^{\prime};\beta)=\varphi_{T}(R,R^{\prime};\beta)$, i.e., with its phase fixed. Then Eq. (59) is formally equivalent to a Bloch equation for $|\rho(R,R^{\prime};\beta)|$ with effective potential ($R^{\prime}$ is fixed)

[TABLE]

Thus PIMC with phase fixing samples paths with real and nonnegative weight, using a fixed-phase dependent effective interaction.

If we know $\varphi_{T}(R_{m},R_{m-1};\beta)$ and its gradient $\nabla_{R_{m}}\varphi_{T}(R_{m},R_{m-1};\beta)$, we can apply the following approximation. The gradient is decomposed into components parallel and perpendicular to the semiclassical path between $(R_{m-1},0)$ and $(R_{m},\tau)$:

[TABLE]

The perpendicular component is taken into account by the primitive action. On the other hand, the evolution of the phase is approximated by a cubic polynomial on the semiclassical trajectory, and the contribution of the parallel component of the gradient of $\varphi_{T}$ is integrated on this trajectory as in the semiclassical approximation to the action. Technically, we assume the following quantities are known:

[TABLE]

and $g_{1}^{\perp}$, $g_{1}^{\parallel}$, $g_{2}^{\perp}$, $g_{2}^{\parallel}$ are magnitudes of the perpendicular and parallel components of $g_{1}$ and $g_{2}$, respectively, in the sense of Eq. (61). (If the phase is fixed to a single-particle density matrix, $g_{1}=-y^{\prime}\mathbf{\hat{x}}/\ell^{2}$ both for open and periodic boundary conditions. If the phase of the free Fermi or Bose gas is used, cf. Eqs. (36-37), $g_{1}=-\sum_{i}y^{\prime}_{i}\mathbf{\hat{x}_{i}}/\ell^{2}$.)

The perpendicular component is taken into account by the primitive action:

[TABLE]

The next contribution is the line integral of $(G^{\parallel})^{2}$ on the straight path between $R_{m-1}$ and $R_{m}$, if $\varphi_{T}$ is approximated by a cubic polynomial on this route.

[TABLE]

We proceed in the same way for the dot product of the phase gradient and the vector potential. $A\cdot G^{\perp}$ contributes at the end points:

[TABLE]

and for $A\cdot G^{\parallel}$ we again use the semiclassical action with the cubic approximation for $\varphi_{T}$:

[TABLE]

where $c_{2}=\frac{3(\varphi_{2}-\varphi_{1})}{\delta R^{2}}-\frac{2g^{\parallel}_{1}+g^{\parallel}_{2}}{\delta R}$ and $c_{3}=\frac{g^{\parallel}_{1}+g^{\parallel}_{2}}{\delta R^{2}}-\frac{2(\varphi_{2}-\varphi_{1})}{\delta R^{3}}$. Finally, the semiclassical contribution of the $A^{2}$ term is

[TABLE]

The total contribution is the sum of Eqs. (63) to (67).