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

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 multislice 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.

Figure 1

(a) The principal domain of the torus. We also depict ${\mathbf{L}}_1$, ${\mathbf{L}}_2$, and $\theta$ as defined in the text; we identify the plane with the complex plane, and we indicate the corners of the principal domain using the complex parameter $\tau$ defined in Eq. (4). (b) The quadruple domain used for finding the zeros of the density matrix in the low-temperature limit.

Figure 2

The dependence of $|{\rho }^{\text{PBC}}(\mathbf{r},{\mathbf{r}}^{\prime };\beta )|$ on imaginary time $\beta$. There are $N_{\varphi }=6$ flux quanta in the principal domain, $L_2/L_1=1.17,\theta \approx {55}^{\circ }$, ${\varphi }_1={\varphi }_2=0$, and we have fixed ${\mathbf{r}}^{\prime }$ at the origin. The panels correspond to $\beta \hslash {\omega }_c=0.3$, 0.7, 1.1, and 5, respectively.

Figure 3

The low-temperature limit $\beta \to \infty$ of the thermal density matrix. As no change is discernible beyond $\beta =100$, the density plot has been generated using this value. (a) $N_{\varphi }=4$ particles, $\theta \approx {65}^{\circ },|{\mathbf{L}}_2|/|{\mathbf{L}}_1|=1.2,{\varphi }_1={\varphi }_2=0$. In the zoomed area, we show how the zeros move to their asymptotic position as a function of inverse temperature $\beta$. (b) The same for $N_{\varphi }=5,\theta \approx {50}^{\circ },|{\mathbf{L}}_2|/|{\mathbf{L}}_1|=1.25,{\varphi }_1={\varphi }_2=0$. Note that one of the zeros is fixed at the corner of the principal region, which is the generic behavior when $N_{\varphi }$ is odd.

Figure 4

The structure of zeros of the thermal density matrix for $\beta \hslash {\omega }_c=200$, where the picture is stationary for different geometries and flux quanta. We show $|{\rho }^{\text{PBC}}(\mathbf{r},{\mathbf{r}}^{\prime };\beta )|$; the zeros are the darkest spots. We set ${\varphi }_1={\varphi }_2=0$ and fix ${\mathbf{r}}^{\prime }$ at the origin. (a) Generic torus with $N_{\varphi }=6,L_2/L_1=1.13$, and $\theta \approx {75}^{\circ }$; (b) square principal domain ($\theta ={90}^{\circ },L_2/L_1=1$) with $N_{\varphi }=7$; (c) generic torus with $N_{\varphi }=11,L_2/L_1=1.19$, and $\theta \approx {72}^{\circ }$; (d) hexagonal principal domain ($\theta ={60}^{\circ },L_2/L_1=1$) with $N_{\varphi }=12$.

Figure 5

The trajectories of the zeros of the thermal density matrix as we tune (a) the twist angle ${\varphi }_1$ and (b) the twist angle ${\varphi }_2$ between 0 and $2\pi$. We set $\beta \hslash {\omega }_c=200,{\mathbf{r}}^{\prime }=\mathit{0},N_{\varphi }=2,L_2/L_1=1.19$, and $\theta \approx {56}^{\circ }$. The marks at specific points on the trajectories correspond to multiples of $\pi /5$. The speed of the zeros is not uniform, as is visible from the distance between adjacent labeled points.

Figure 6

Acceptance ratios for sampling the motion of a single particle on a rectangular torus pierced by $N_{\varphi }=2$ flux quanta. The inverse temperature of the system is $\beta \hslash {\omega }_c=2$, and the number of slices ranged between $M=8$ and 256. Level $l$ means that $2^l-1$ slices are updated in each multislice move.

Figure 7

(a)–(c) The pair-correlation function for $N=12$ bosons at density $\rho a^2=0.02$ ($\kappa =0.25066$) at filling factor $\nu =2$ (i.e., $N_{\varphi }=6$ flux quanta piercing the torus) and interaction strength $\mathrm{\Lambda }=0.035$, 0.04, and 0.045, respectively. $M=32$ slices were used, the imaginary time step is $\tau =0.015625$. The temperature is low on the scale of interactions, as $\tilde{\beta }=114$, 99, and 88 in panels (a)–(c). Panels (d) and (e) show the differences of the pair-correlation functions, $g_{\mathrm{\Lambda }=0.04}-g_{\mathrm{\Lambda }=0.035}$ and $g_{\mathrm{\Lambda }=0.045}-g_{\mathrm{\Lambda }=0.04}$, respectively, as $\mathrm{\Lambda }$ is changed for systems shown in the top row. The triangular lattice of dark spots shows the decreasing crystalline correlation as $\mathrm{\Lambda }$ is increased. The small deviations from perfect $C_6$ symmetry in panels (b) and (c) can be attributed to imperfect thermalization, and could be reduced by longer Monte Carlo runs. Taking the differences between pair-correlation functions in panels (d) and (e) amplifies these small errors.

Figure 8

Second row (d)–(f): the pair-correlation function for $N=16$ fermions at density $\rho a^2=0.02$ ($\kappa =0.25066$) at filling factor $\nu =2$ (i.e., $N_{\varphi }=8$ flux quanta piercing the torus) and interaction strength $\mathrm{\Lambda }=0.035$ at inverse temperature ${\beta }^{*}=0.4$, 0.5, and 0.6, respectively. In (d) $M=16$ slices were used, $\tau =0.025$, and $\tilde{\beta }=132$; in (e) $M=16,\tau =0.03125$, and $\tilde{\beta }=114$; in (f) $M=24,\tau =0.025$, and $\tilde{\beta }=99$. Panels (a) and (c) show the differences of the pair-correlation functions $g_{\beta =0.4}-g_{\beta =0.5}$ and $g_{\beta =0.6}-g_{\beta =0.5}$, respectively, between colder and warmer systems shown in consecutive panels in the second row. The triangular 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 ${\beta }^{*}=0.5$, but the de Boer parameter is tuned from $\mathrm{\Lambda }=0.03$ in panel (h) to $\mathrm{\Lambda }=0.04$ in panel (b). Panels (g) and (i) show the differences of the pair-correlation functions $g_{\mathrm{\Lambda }=0.035}-g_{\mathrm{\Lambda }=0.03}$ and $g_{\mathrm{\Lambda }=0.04}-g_{\mathrm{\Lambda }=0.035}$, respectively, as $\mathrm{\Lambda }$ is tuned for systems shown in the second column. The triangular lattice of dark spots shows the decreasing crystalline correlation as $\mathrm{\Lambda }$ is increased.

Figure 9

The height of the first peak of the pair-correlation function for $N=12$ fermions at ${\beta }^{*}=0.5$ and $N_{\varphi }=6$ for various $\mathrm{\Lambda }$ de Boer parameter values as a function of density $\rho a^2$. The nonmonotonic evolution indicates that crystalline order exists only for a limited range of densities.

Figure 10

The evolution of the first peak of the pair-correlation function for $N=16$ fermions at flux $N_{\varphi }=8$ at density $\rho a^2=0.02$ ($\kappa =0.25066)$, as a function of the inverse temperature for some values of the de Boer parameter $\mathrm{\Lambda }$ for which crystalline structure is manifest at intermediate temperatures. A small horizontal shift has been applied to the last two curves to make the overlapping error bars visible.

Figure 11

The difference $g_{N_{\varphi }=2}\left(\mathbf{r}\right)-g_{N_{\varphi }=6}\left(\mathbf{r}\right)$ between pair-correlation functions at different flux densities (6 and 2 flux quanta through the torus) at ${\beta }^{*}=0.5$ for $N=16$ and $\rho a^2=0.02$ for fermions (a) and bolzmannons (b). The total area of the simulation call is scaled to unity, thus the peak locations may coincide. (The filling factor corresponding to $N_{\varphi }=6,2$ is $\nu =\frac{8}{3},8$, respectively.) The triangular lattice of bright peaks corresponds to stronger crystalline correlations at smaller flux density. Bolzmannons in panel (b) are still liquidlike; a small rotation of the hexagonally distorted rings from directions where crystalline structure will emerge can be attributed to imperfect thermalization.