Path-integral Monte Carlo study of electronic states in quantum dots in an external magnetic field

TL;DR

This study uses path-integral Monte Carlo simulations to investigate correlated electron states in quantum dots under magnetic fields, revealing complex behaviors beyond mean-field predictions and providing insights into ground-state properties.

Contribution

It introduces a novel application of path-integral Monte Carlo to analyze correlated states in quantum dots under magnetic fields, highlighting states that differ from mean-field theory.

Findings

Discovery of diverse correlated states in quantum dots

Correlations sometimes saturate at low temperatures

Results provide insights into ground-state properties

Abstract

We explore correlated electron states in harmonically confined few-electron quantum dots in an external magnetic field by the path-integral Monte Carlo method for a wide range of the field and the Coulomb interaction strength. Using the phase structure of a preceding unrestricted Hartree-Fock calculation for phase fixing, we find a rich variety of correlated states, often completely different from the prediction of mean-field theory. These are finite temperature results, but sometimes the correlations saturate with decreasing temperature, providing insight into the ground-state properties.

Equations

$$\begin{gathered} \mathcal{H}=\sum _{i=1}^N\frac{({\mathbf{p}}_i-e\mathbf{A}({\mathbf{r}}_i){)}^2}{2m^{\ast }}+\frac{1}{2}{m}^{\ast }{\omega }_0^2\sum _{i=1}^N{r}_i^2 \\ +\frac{e^2}{4\pi {ϵ}_r{ϵ}_0}\sum _{i<j}\frac{1}{|{\mathbf{r}}_i-{\mathbf{r}}_j|}, \end{gathered}$$

$${\mathcal{H}}_0=-\frac{{\hslash }^2}{2m^{\ast }}{\nabla }^2+\frac{1}{2}{m}^{\ast }{\omega }^2{r}^2-\frac{{\omega }_c}{2}{L}_z,$$

$${\eta }_{nl}(\mathbf{r})=\frac{e^{il\varphi }}{\sqrt{2\pi }}\sqrt{\frac{n!}{(n+|l|)!}}{e}^{-r^2/4}{\left(\frac{r}{\sqrt{2}}\right)}^{|l|}{L}_n^{|l|}\left(\frac{r^2}{2}\right),$$

$${ϵ}_{nl}=\hslash \omega (2n+|l|+1)-\frac{\hslash {\omega }_c}{2}l.$$

$$\lambda =\frac{e^2/(4\pi {ϵ}_r{ϵ}_0\ell )}{\hslash \omega }=\frac{\ell }{a_{\text{B}}^{\ast }},$$

$$\gamma =\frac{{\omega }_c}{\omega }.$$

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

$$\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}$$

$$V_{\text{eff}}(R)=\frac{{\hslash }^2}{m^{\ast }}{\left(\nabla {\phi }_T(R,R^{\prime };\beta )-\frac{e}{\hslash }A(R)\right)}^2$$

$$\begin{gathered} {\rho }_0(\mathbf{r},{\mathbf{r}}^{\prime };\beta )=\frac{1}{2\pi {\ell }_c^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 }_c^2}+\frac{i(x^{\prime }y-x{y}^{\prime })}{2{\ell }_c^2}\right), \end{gathered}$$

$$U_{\text{C}}(R_0,R_1,\tau )=\tau \int dR{\mu }_{\text{C}}(R|R_0,R_1,\tau )V(R),$$

$${\mu }_{\text{C}}(R|R_0,R_1,\tau )=\frac{1}{\tau }{\int }_0^{\tau }dt\mu (R,t|R_0,R_1,\tau ),$$

$$\begin{gathered} \mu (R,t|R_0,R_1,\tau )= \\ \frac{\left({\prod }_{i=1}^N{\rho }_0({\mathbf{r}}_{0,i},{\mathbf{r}}_i;t)\right)\left({\prod }_{i=1}^N{\rho }_0({\mathbf{r}}_i,{\mathbf{r}}_{1,i};\tau -t)\right)}{{\prod }_{i=1}^N{\rho }_0({\mathbf{r}}_{0,i},{\mathbf{r}}_{1,i};\tau )}. \end{gathered}$$

$$\{{\eta }_i^{\text{HF},\uparrow }(\mathbf{r}){\}}_{i=1,\dots ,N_{\text{basis}}}\text{and}\{{\eta }_i^{\text{HF},\downarrow }(\mathbf{r}){\}}_{i=1,\dots ,N_{\text{basis}}}$$

$$\begin{gathered} {\Psi }_{{\Omega }_{\uparrow },{\Omega }_{\downarrow }}^{\text{HF}}(R)=\text{Det}\left[{\eta }_i^{\text{HF},\uparrow }({\mathbf{r}}_j{)}_{i\in {\Omega }_{\uparrow },1\le j\le N_{\uparrow }}\right] \\ \times \text{Det}\left[{\eta }_i^{\text{HF},\downarrow }({\mathbf{r}}_{N_{\uparrow }+j}{)}_{i\in {\Omega }_{\downarrow },1\le j\le N_{\downarrow }}\right]. \end{gathered}$$

$${\rho }_{\text{T,0}}(R,R^{\prime };\beta )={\Psi }_0^{\text{HF}}(R){\Psi }_0^{\text{HF}}(R^{\prime }{)}^{\ast }.$$

$${\rho }_{\text{T}}(R,R^{\prime };\beta )=\sum _{{\Omega }_{\uparrow },{\Omega }_{\downarrow }}{e}^{-\beta E_{{\Omega }_{\uparrow },{\Omega }_{\downarrow }}}{\Psi }_{{\Omega }_{\uparrow },{\Omega }_{\downarrow }}^{\text{HF}}(R){\Psi }_{{\Omega }_{\uparrow },{\Omega }_{\downarrow }}^{\text{HF}}(R^{\prime }{)}^{\ast },$$

$$g(\theta )=⟨\sum _{i<j}\delta (\theta -{\theta }_{i,j})⟩,$$

$$\rho (r)=\frac{1}{2\pi }{\int }_0^{2\pi }\rho (r,\varphi )d\varphi ,$$

$$\begin{gathered} U_{\text{C}}^{\text{conf.}}(r_{i,0},r_{i,1},\tau )=\frac{m^{\ast }{\omega }_0^2{\ell }_c^2}{4(1-u)}[\frac{2(r_{i,0}^2+r_{i,1}^2)}{{\ell }_c^2(1-u)}{I}_2 \\ +\left(4+\frac{i}{2{\ell }_c^2}(x_0{y}_1-x_1{y}_0)+\frac{1+u}{1-u}\frac{2{\mathbf{r}}_{i,0}\cdot {\mathbf{r}}_{i,1}}{{\ell }_c^2}\right)I_1], \end{gathered}$$

$$\begin{array}{rl}{I}_1& =\tau (1+u)+\frac{2}{\hslash {\omega }_c}(u-1),\\ I_2& =\frac{1-u^2}{\hslash {\omega }_c}-2\tau u.\end{array}$$

$$U_{\text{C}}(R_0,R_1,\tau )=\sum _{i<j}{u}_{\text{C}}^{i,j}$$

$$\begin{gathered} u_{\text{C}}^{i,j}={\int }_0^{\tau }dt\int d^2{r}_i{d}^2{r}_j\frac{e^2}{4\pi ϵ|{\mathbf{r}}_i-{\mathbf{r}}_j|} \\ \times \mu ({\mathbf{r}}_i,{\mathbf{r}}_j,t|{\mathbf{r}}_{0,i},{\mathbf{r}}_{0,j},{\mathbf{r}}_{1,i},{\mathbf{r}}_{1,j},\tau ) \end{gathered}$$

$$\begin{gathered} \mu (\dots )=\frac{{\rho }_0({\mathbf{r}}_{0,i},{\mathbf{r}}_i;t){\rho }_0({\mathbf{r}}_i,{\mathbf{r}}_{1,i};\tau -t)}{{\rho }_0({\mathbf{r}}_{0,i},{\mathbf{r}}_{1,i};\tau )} \\ \times \frac{{\rho }_0({\mathbf{r}}_{0,j},{\mathbf{r}}_j;t){\rho }_0({\mathbf{r}}_j,{\mathbf{r}}_{1,j};\tau -t)}{{\rho }_0({\mathbf{r}}_{0,j},{\mathbf{r}}_{1,j};\tau )}. \end{gathered}$$

$$\mu (\dots )=\mu (\mathbf{R},t|{\mathbf{R}}_0,{\mathbf{R}}_1,\tau )\mu (\mathbf{r},t|{\mathbf{r}}_0,{\mathbf{r}}_1,\tau ),$$

$$\begin{gathered} u_{\text{C}}^{i,j}=\frac{1}{2}{\int }_0^{\tau }dt\int \frac{d^{2}k}{(2\pi {)}^2}V(k) \\ \times \int d^{2}r\mu (\mathbf{r},t|{\mathbf{r}}_0,{\mathbf{r}}_1,\tau )e^{-i\mathbf{k}\cdot \mathbf{r}}. \end{gathered}$$

$$\begin{array}{r}{\sigma }_t^2=\frac{(1-u_1)(1-u_2)}{(1-u_0)},\Delta \mathbf{r}={\mathbf{r}}_1-{\mathbf{r}}_0,\\ {\mathbf{r}}_t={\mathbf{r}}_0+\frac{(1-u_1)(1+u_2)}{2(1-u_0)}\Delta \mathbf{r},\end{array}$$

$$\begin{gathered} u_{\text{C}}^{i,j}={\int }_0^{\tau }dt\int \frac{d^{2}k}{(2\pi {)}^2}V(k) \\ \times \exp \left(-{\sigma }_t^2{\ell }_c^2{k}^2-\frac{{\sigma }_t^2}{2}\hat{\mathbf{z}}\cdot \left(\mathbf{k}\times \Delta \mathbf{r}\right)-i\mathbf{k}\cdot {\mathbf{r}}_t\right). \end{gathered}$$

$$\begin{gathered} u_{\text{C}}^{i,j}={\int }_0^{\tau }dt{\int }_0^{\infty }\frac{kdk}{(2\pi {)}^2}V(k)e^{-{\sigma }_t^2{\ell }_c^2{k}^2} \\ \times {\int }_0^{2\pi }d\varphi \exp \left(-\frac{{\sigma }_t^2}{2}k\Delta r\cos ({\beta }_t+\varphi )-ik{r}_t\cos (\varphi )\right). \end{gathered}$$

$$\frac{{\ell }_{\tau }^2}{{\ell }_c^2}\ll 1,$$

Full text

I Introduction

Quantum dots established in gated submicron regions of a two-dimensional electron gas, contacted either vertically or laterally, have been studied intensively in the last decades Chakraborty (1992); Kouwenhoven et al. (2001); Reimann and Manninen (2002); Alhassid (2000). Such structures trap a few to a few hundred electrons, and are often referred to as artificial atoms. Unlike in real atoms, where the electron-electron interaction plays a quantitatively important but qualitatively minor role, in weakly confined quantum dots the electron density can be tuned to small values, and the Coulomb repulsion among electrons can play a prominent role. Many experiments are adequately described in terms of weak-coupling theories (Hartree-Fock, density functional theory, or the constant interaction model). This is no doubt due to the strong screening of the Coulomb interaction, most importantly, by the leads in vertical quantum dots. With lateral gates, however, the intermediate coupling regime is accessible, while the strongly coupled regime requires weaker confinement in cleaner samples, or the use of semiconductors that have an effective mass larger than GaAs.

Electrons in quantum dots show an especially intriguing behavior in the presence of an external magnetic field, as the relevant length scales—the confinement length, the magnetic length, and the effective Bohr radius of the interaction—can be tuned to comparable values, letting confinement, interaction and magnetic effects compete. In the presence of an external magnetic field, the ground-state properties have been studied by exact diagonalization (in truncated Hilbert spaces) Maksym and Chakraborty (1990); Wagner et al. (1992); Mikhailov and Savostianova (2002); Cha and Yang (2003); Tavernier et al. (2003, 2006), Hartree-Fock mean-field theory Yannouleas and Landman (1999, 2000a); Reusch et al. (2001); Reusch and Grabert (2003); Yannouleas and Landman (2003a), post-Hartree-Fock projection techniques Yannouleas, C. and Landman, U. (2001); Yannouleas and Landman (2002a, b, 2003b); Giovannini et al. (2007); De Giovannini et al. (2008), and the diffusion Monte Carlo method with phase fixing Bolton (1996); Güçlü and Umrigar (2005); Güçlü et al. (2005). Moreover, at finite temperature, the energy scale $k_{\text{B}}T$ of thermal excitations brings a fourth player into the game. Finite temperature theoretical studies rely upon either the finite temperature Hartree-Fock theory Dean et al. (2001) or the exact diagonalization of the complete truncated Hamiltonian Maksym and Chakraborty (1990) (as distinct from finding the ground state and the lowest excitations by the Lánczos method).

The path-integral Monte Carlo (PIMC) method has been applied to quantum dots only for zero external magnetic fields Egger et al. (1999a, b); Harting et al. (2000); Filinov et al. (2001); Lozovik et al. (2003); Reusch and Egger (2003); Weiss and Egger (2005). No doubt this is due to the sign problem, which arises in quantum Monte Carlo simulations of fermionic systems, and the typically more severe phase problem, which occurs whenever the presence of an external magnetic field or spin-orbit coupling prevents us from using real-valued wave functions and density matrices, irrespective of the statistics of the simulated particles. For zero magnetic field, the multilevel blocking method has been shown to mitigate the sign problem Mak et al. (1998); Egger et al. (1999a); Egger and Mak (2000), but the generalization of this method to magnetic problems does not seem to be straightforward. It is also possible to use phase fixing in path-integral Monte Carlo; this requires the use of a guiding many-body density matrix. Recently we have shown Haidekker Galambos and Tőke (2018) that a relatively crude way of fixing the phase may nevertheless yield valuable results in PIMC. Encouraged by this finding, here we apply PIMC for characterizing the electronic structures that arise in harmonically confined dots. We fix the phase by using the outcome of the unrestricted Hartree-Fock (UHF) approximation to the many-body problem. For tiny systems we could use the phase of the density matrix of the Hartree-Fock effective Hamiltonian. But we can study much larger systems by using the density matrix built up of the Hartree-Fock ground state only; this can be regarded as the zero-temperature limit of the former approximation. As we will see, a phase-fixing ansatz still leaves sufficient freedom in the PIMC method to let the simulated system discover strongly correlated phases that qualitatively differ from the guide density matrix that is applied for fixing the phase. This situation is quite different from the use of phase fixing in diffusion Monte Carlo Ortiz et al. (1993); Bolton (1996); Güçlü et al. (2005); Ghosal et al. (2006), a method that becomes very inefficient unless the trial wave function is a good approximation to the true ground state, as the trial wave function is also used for importance sampling during the random walk.

We emphasize that we are interested both in the performance of the simulation method and the correlated behavior of the electron system. We are particularly interested in methods that make the small magnetic field range accessible. This range is especially important in experiments Zhitenev et al. (1997). But most efficient correlated methods work better in the large magnetic field limit, e.g., exact diagonalization (configuration interaction) as restricted to the lowest Landau level Maksym and Chakraborty (1990), or composite fermion diagonalization Jeon et al. (2004, 2007). We will see that PIMC is still efficient when the cyclotron energy is about 20% of the confinement energy, which for typical experimental parameters corresponds to less than 0.4 T Zhitenev et al. (1997), and the lowest Landau level approximation is definitely invalid. However, below this range PIMC also becomes increasingly difficult.

Preliminary information on the ordering within the quantum dot is obtained by inspecting the spatial structure of the UHF ground state. Within PIMC, on the other hand, only the correlation functions yield structural information, as the spatial profiles are subject to a random reorientation during the random walk process involved in the simulation, apparently not hindered by phase fixing. We find that the predictions of the PIMC method are often at odds with the UHF theory, even at small temperatures, where essentially ground-state properties are expected to dominate. For example, at the coupling strength where spin polarization becomes complete, UHF predicts the reentrance of featureless rotationally invariant states, but PIMC discerns the incipient structure of the Wigner molecule in the correlations. In some cases, PIMC completely revises the spatial structure suggested by UHF. For example, for $N=6$ electrons PIMC detects a transition from a hexagonal Wigner molecule to a pentagonal one at strong coupling at some magnetic fields. PIMC also finds correlated structures that never occur in UHF, as we will see below. These findings also affect the validity of the post-Hartree-Fock wave functions derived for the Wigner molecules as obtained by projecting the symmetry-breaking UHF state to a fixed angular momentum subspace Yannouleas, C. and Landman, U. (2001); Yannouleas and Landman (2002a, b, 2003b); Giovannini et al. (2007); De Giovannini et al. (2008), as these wave functions inherit the symmetry of the UHF ground state.

In spite of the valuable structural information, however, the PIMC method with phase fixing is unable to yield a complete revised phase diagram at present. The reason is that PIMC is performed in subspaces of fixed $z$ component of the total spin, and in the absence of a reliable estimator of the (conditional) free energy in these subspaces, the information on the spin remains limited.

The paper is structured as follows. In Sec. II, we introduce the quantum dot model we study, and define the dimensionless parameters that are relevant in the analysis. In Sec. III, we review the path-integral Monte Carlo method, the unrestricted Hartree-Fock method, and the use of phase fixing in PIMC. In Sec. IV, we present the predictions of these methods in detail for small dots, containing $N=3$ (quantum dot lithium) to $N=8$ electrons (quantum dot oxygen). We summarize our findings and conclude on the phase-fixed PIMC method in Sec. V. The cumulant approximation to the action in the presence of an external magnetic field, which is the main technical invention of this paper, is derived in the Appendix both for the Coulomb interaction and the harmonic confinement potential.

II Model

We study a two-dimensional quantum dot of rotational symmetry; we assume that the confinement potential is parabolic. The Hamiltonian is

[TABLE]

where we use the symmetric gauge $\mathbf{A}(\mathbf{r})=\frac{B}{2}(-y,x,0)$; $\hbar\omega_{0}$ is the confinement energy, and $\epsilon_{r}$ and $m^{\ast}$ are material-specific parameters (e.g., 12.7 and 0.067 me in GaAs, respectively). The noninteracting part of the Hamiltonian can also be written, for a single particle for simplicity, as

[TABLE]

where $\omega_{c}=\frac{eB}{m^{\ast}}$ is the cyclotron frequency, $\omega=\sqrt{\omega_{0}^{2}+\omega_{c}^{2}/4}$, and $L_{z}$ is the angular momentum perpendicular to the plane of the electron system. We will use $\ell=\sqrt{\frac{\hbar}{m^{\ast}\omega}}$ as the unit of length throughout this paper. The eigenstates of $\mathcal{H}_{0}$ are the well-known Fock-Darwin orbitals,

[TABLE]

where $l$ is the eigenvalue of $L_{z}$. The corresponding energy is

[TABLE]

In the limit $\omega_{0}\to 0$, $\epsilon_{nl}$ tends to the $l$-independent quantized energies of Landau levels, while in the $B\to 0$ limit it becomes identical to the energy of the two-dimensional harmonic oscillator.

The interacting electron system in the quantum dot is characterized by two dimensionless parameters: the ratio of the Coulomb energy scale to the confinement,

[TABLE]

where $a^{\ast}_{\text{B}}=\frac{\hbar^{2}4\pi\epsilon_{r}\epsilon_{0}}{m^{\ast}e^{2}}$ is the effective Bohr radius of the host semiconductor, and the dimensionless measure of the magnetic field,

[TABLE]

[Sometimes the ratios $\lambda_{0}=\frac{e^{2}/(4\pi\epsilon_{r}\epsilon_{0}\ell_{0})}{\hbar\omega_{0}}$ and $\gamma_{0}=\frac{\omega_{c}}{\omega_{0}}$ are used, where $\ell_{0}=\sqrt{\frac{\hbar}{m^{\ast}\omega_{0}}}$ is the oscillator length of the confinement potential. Experimentally, neither $\omega$ nor $\omega_{0}$ is constant, as the latter depends on the gate voltage and the former also depends on the magnetic field.]

III Methods

The path-integral Monte Carlo method Ceperley (1995) evaluates physical quantities as derivatives of the partition function; the latter is expressed as an imaginary-time path integral, and path integration is performed by standard Monte Carlo methods. We consider the thermal density matrix of the system,

[TABLE]

where $\beta=\frac{1}{k_{\text{B}}T}$ is the inverse temperature, $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. We express $\rho(R,R^{\prime};\beta)$ in terms of density matrices that correspond to a higher temperature:

[TABLE]

Here, $\tau\equiv\beta/M$ is the imaginary time step between adjacent slices. There are several plausible schemes for approximating $\rho(R_{m-1},R_{m};\tau)$ Ceperley (1995). The simplest of these, the so-called primitive action, is based on the Suzuki-Trotter formula, and it results in the product of a purely kinetic and an interaction part. Unfortunately, the primitive approximation to the action is well-known to be inadequate for the Coulomb interaction because of the sharp divergence of the latter at short range Ceperley (1995). Here, we adopt the cumulant approximation to the action Ceperley (1983). Its standard formulation, however, has to be changed in the presence of an external magnetic field, as sketched below and discussed in detail in Appendix.

In the presence of an external magnetic field, the density matrix unavoidably becomes complex-valued. The estimators of physical quantities are then obtained as sums of complex values. This results in an extremely low signal-to-noise ratio in the Monte Carlo integration, even for small systems. We avoid this problem by phase fixing Ortiz et al. (1993); Bolton (1996); Haidekker Galambos and Tőke (2018). In the case of PIMC, we have to fix the phase of the many-body density matrix in Eq. (7), not of a wave function. The paths are sampled by the probability density function $\prod_{m=1}^{M}|\rho(R_{m-1},R_{m};\tau)|$, where we assume that the diagonal density matrix is the integration kernel of the estimator, and we let $R=R^{\prime}\equiv R_{0}=R_{M}$. Using a trial (or guide) many-body density matrix $\rho_{T}(R,R^{\prime};\beta)=|\rho_{T}(R,R^{\prime};\beta)|e^{i\varphi_{T}(R,R^{\prime};\beta)}$ for fixing the phase, this means that an additional effective interaction

[TABLE]

appears in the action, besides the physical potential terms. Notice that $|\rho_{T}(R,R^{\prime};\beta)|$ plays no role in this procedure. While for the effective interaction a semiclassical approximation is used Haidekker Galambos and Tőke (2018), i.e., $V_{\text{eff}}$ is evaluated on a straight line that connects $R$ and $R^{\prime}$, for the harmonic confinement and for the Coulomb repulsion among electrons, we use the cumulant approximation to the action. Single-particle propagation in the presence of an external magnetic field is described by the free density matrix

[TABLE]

where $\ell_{c}=\sqrt{\frac{\hbar}{eB}}$ is the magnetic length and $u=e^{-\beta\hbar\omega_{c}}$. The cumulant approximation is

[TABLE]

where the complex amplitude of visiting configuration $R$ sometime during the propagation of the electron configuration from $R_{0}$ to $R_{1}$ is

[TABLE]

with

[TABLE]

It is easy to check that $\mu$ is gauge-independent, which transfers to $\mu_{\text{C}}$ and to the cumulant action $U_{\text{C}}$. The evaluation of Eqs. (11)-(13) is described in detail in Appendix.

We will use the outcome of a preliminary unrestricted Hartree-Fock calculation for phase fixing. The UHF method is quite standard and it has been applied to quantum dot problems many times Pfannkuche et al. (1993); Palacios et al. (1994); Fujito et al. (1996); Yannouleas and Landman (1999, 2000b, 2003a); Reusch et al. (2001); Reusch and Grabert (2003). We just note that we use UHF in the basis of Fock-Darwin states, Eq. (3), calculate the Coulomb matrix elements following Ref. Anisimovas and Matulis, 1998. The number of basis states was $N_{\text{basis}}=55$ for $N\leq 6$ and $N_{\text{basis}}=66$ for $N=7,8$. As the result of the self-consistent solution of the Pople-Nesbet-Berthier equations in each $S_{z}=(N_{\uparrow}-N_{\downarrow})/2$ sector, where $N_{\uparrow}+N_{\downarrow}=N$, two sets of orbitals

[TABLE]

are obtained for each spin. We assume these orbitals are ordered in increasing order of the corresponding eigenvalues $\{\epsilon^{\uparrow}_{i}\}$ and $\{\epsilon^{\downarrow}_{i}\}$ for $i=1,\dots,N_{\text{basis}}$. The UHF wave function for the choice $\Omega_{\uparrow}$, $\Omega_{\downarrow}$ of orbitals is

[TABLE]

In particular, the UHF ground state $\Psi^{\text{HF}}_{0}(R)$ is specified by $\Omega_{\uparrow}=\{1,2,\dots,N_{\uparrow}\}$ and $\Omega_{\downarrow}=\{1,2,\dots,N_{\downarrow}\}$. These mean-field wave functions are unrestricted in the sense that they are not necessarily eigenstates of either the magnitude of the total spin $S^{2}$ or the total angular momentum $L_{z}=\sum_{i=1}^{N}L_{z,i}$. At the same time, the UHF ground state typically has lower energy than the restricted Hartree-Fock ground-state, which has well-defined total spin and angular momentum quantum numbers. UHF correctly predicts the “magic numbers” $N=2,6,12,20,\dots$ in the addition energies, and it justifies Hund’s rules for deciding the ground-state spin at zero magnetic field. A UHF calculation often informs us about the internal structure of the symmetry breaking state even at strong interactions, where the validity of mean-field theories may be questionable. Notably, the post-Hartree-Fock projection techniques Yannouleas, C. and Landman, U. (2001); Yannouleas and Landman (2002a, b, 2003b); Giovannini et al. (2007); De Giovannini et al. (2008) that restore the correct symmetry of the wave function take an unrestricted Hartree-Fock ground state as their starting point.

We stress that a UHF calculation is just a preliminary step to our quantum Monte Carlo approach. Even though a large number of UHF studies of quantum dots exist in the literature, a thorough exploration of the charge and spin structures in terms of the magnetic field ($\gamma$ or $\gamma_{0}$) and the coupling strength of the Coulomb interaction ($\lambda$ or $\lambda_{0}$) is apparently not available, even for small dots. This is understandable, as UHF has important limitations and its predictions should always be handled with care. Nevertheless, as the output of UHF is an input to our correlated PIMC calculations, we need a more or less complete UHF phase diagram, even if it has known shortcomings. For example, UHF overestimates the energy of spin unpolarized states more than that of fully spin polarized ones Pfannkuche et al. (1993), hence it locates the spin transitions at too small $\lambda$’s. At high values of $\lambda$, the lowest UHF energy in distinct $S_{z}$ sectors differ only slightly, hence the prediction of the spin by UHF therefore becomes somewhat ambiguous. The spin-singlet–spin-triplet oscillations for $N=2$ are not reproduced by UHF, nor are the more complicated spin recurrence patterns for higher electron numbers Bolton (1996); Mikhailov and Savostianova (2002); Tavernier et al. (2003). And we will see cases where the spin $S_{z}$ of the UHF ground state in the strong coupling limit decreases from the highest value $S_{z}=N/2$, which occurs when the spatial structure of the ground state also shows an unusual symmetry.

One easy choice for the trial density matrix is in terms of the UHF ground state,

[TABLE]

For small systems, one can try to build a trial density matrix from all UHF wave functions:

[TABLE]

where $E_{\Omega_{\uparrow},\Omega_{\downarrow}}$ is the energy of state $\Psi^{\text{HF}}_{\Omega_{\uparrow},\Omega_{\downarrow}}$.

Note that both the UHF and the PIMC calculations are performed in fixed $S_{z}$ subspaces. As UHF is a ground-state method, the true UHF ground state is found by the comparison of the lowest energy states in each subspace. Thus jumps in $S_{z}$ yield partial information on the total spin $S\geq S_{z}$ of the ground state.

Regarding PIMC, we note that permutations are sampled for up and down-spins separately.

IV Results

We will focus on the effect of the electron-electron interaction and neglect the Zeeman energy completely. The latter is known to be small in GaAs structures ($g=-0.44$), but it is significant in InSb quantum dots ($g\approx-50$). The Zeeman energy can be reintroduced trivially if necessary; it always favors states of greater total spin, but leaves the spatial structure intact within the spin-related phase domains.

Unrestricted Hartree-Fock calculations are performed at fixed $\gamma$ as $\lambda$ is increased from 0 to 9 in steps of 0.1. We use the converged orbitals of the previous $\lambda$ value as the starting point of the iterative solution of the UHF problem at $\lambda+0.1$. When the symmetry of the ground state changes during such an upward sweep, we perform a downward sweep starting with the orbitals of the new symmetry. The goal is to avoid metastable solutions of Pople-Nesbet-Berthier equations. The UHF ground state is found eventually by energy comparison for the two sweeps. The $B=0$ limit is included in the phase diagrams for comparison with the literature, but our primary focus is on $B>0$. As Fig. 1 shows for the case of $N=8$, the ground state is not predominantly in the lowest Landau level, at least in the UHF approximation, for moderate values of $\gamma$.

In PIMC simulations, we use a time step $\tau\hbar\omega=0.01$, hence the number of slices is between 80 and 320 for $\beta^{\ast}=\beta\hbar\omega=0.8$ to 3.2. We have checked that the predictions stabilize at this time step. We update 7, 15, or 31 slices in a multi-slice move, as explained in Ref. Haidekker Galambos and Tőke, 2018.

While ordering is manifest in the charge and spin densities by UHF, we cannot expect this to be the case by PIMC. We gain insight into the ordering through the angular correlation function (ACF)

[TABLE]

where $\theta_{i,j}=\cos^{-1}\left(\frac{\mathbf{r}_{i}\cdot\mathbf{r}_{j}}{r_{i}r_{j}}\right)$ is the angle between two electron positions, taken as nonnegative. We will also use the spin-up ACF $g_{\uparrow\uparrow}(\theta)$, the spin-down ACF $g_{\downarrow\downarrow}(\theta)$, the different-spin ACF $g_{\uparrow\downarrow}(\theta)$, and the same-spin ACF $g_{\text{same}}(\theta)=g_{\uparrow\uparrow}(\theta)+g_{\downarrow\downarrow}(\theta)$. We also often plot the angle-averaged radial density,

[TABLE]

and its spin-specific versions $\rho_{\uparrow}(r)$ and $\rho_{\downarrow}(r)$.

Notice that the quantum dot problem can also be discussed classically Bedanov and Peeters (1994). For some particle numbers, e.g., $N=3,7,10,12,19$, we get especially stable configurations, which can be interpreted as finite pieces of a triangular lattice; while for others like $N=5,6,8$, there are several classical configurations that are very close in energy Reusch and Grabert (2003). The quantum mechanical picture is definitely more complicated, primarily because of zero point motion and the presence of the spin. Nevertheless we will see below that the correlated behavior is more complex in the latter cases.

IV.1 Quantum dot lithium, $N=3$

For three particles there are four relevant ground-state structures. We denote them by the combination of the $z$-component of the spin and the symmetry group of the spin density. (Note that this may not be the same as the symmetry of the wave function.) (i) The fully polarized, rotationally invariant state is called $(\frac{3}{2},C_{\infty})$. (ii) The fully polarized, broken symmetry state is called $(\frac{3}{2},C_{3v})$. (iii) The partially polarized, rotationally invariant state is called $(\frac{1}{2},C_{\infty})$. (iv) The partially polarized, broken symmetry state is called $(\frac{1}{2},C_{s})$. Figure 2 shows examples for the symmetry-breaking states $(\frac{3}{2},C_{3v})$ and $(\frac{1}{2},C_{s})$.

For $B>0$, all four structures are relevant in the UHF approximation. Figure 3 gives the UHF phase diagram. At fixed $\gamma$, the four states typically occur in the sequence $(\frac{1}{2},C_{\infty})$, $(\frac{1}{2},C_{s})$, $(\frac{3}{2},C_{\infty})$, $(\frac{1}{2},C_{s})$ in increasing order of $\lambda$, with two exceptions. (i) For $0.4\lesssim\gamma\lesssim 0.5$, the $S_{z}=\frac{3}{2}$ system is already symmetry-broken when the ground state becomes spin-polarized, thus the spin-flip at $\lambda=2.8$ connects two symmetry-broken states. (ii) For $\gamma\geq 1.5$ the range of $(\frac{1}{2},C_{s})$ disappears, i.e., the symmetry-breaking occurs only within the fully spin-polarized range; this behavior is typical for large $\gamma$ for all system sizes.

IV.1.1 The $S_{z}=\frac{3}{2}$ subspace of $N=3$

The reentrance of rotational symmetry with the spin-flip transition deserves further attention. Hence we performed systematic PIMC simulations for all $\gamma>0$ where this occurred.

Figure 4 shows that the angular correlation is strongly peaked near $\theta=2\pi/3$ already at the smallest coupling $\lambda$ where the spins are fully polarized, i.e., in the vicinity of the $(\frac{1}{2},C_{s})$-$(\frac{3}{2},C_{\infty})$ phase boundary. This correlation peak grows further with increasing $\lambda$, and it is accompanied by the relative depletion of the density near the origin, in comparison with the UHF ground state. The correlated motion of electrons also allows for a smaller diameter of the cloud. The angular correlation evolves continuously across the $\lambda$ value where the $C_{3v}$ structure becomes the UHF ground state: no sudden change is observed, even though for larger $\lambda$’s we take the guiding phase from the symmetry-breaking UHF ground state. We conclude that the threefold ordering that corresponds to a triangular Wigner molecule is present throughout the spin-polarized regime; the restored rotational symmetry is a misleading result of the UHF approximation. Note that $g(\theta)$ is independent of the inverse temperature $\beta^{\ast}$ for $\gamma=0.2$, but for $\gamma=0.6$ it saturates for $\beta^{\ast}\geq 1.6$ only. Hence the energy scale associated with threefold ordering in the later case can be estimated. Using $\hbar\omega_{0}=3.32$ meV, commonly assumed when analyzing Ref. Ashoori et al., 1992, this corresponds to $T=\frac{2\hbar\omega_{0}}{\beta^{\ast}k_{\text{B}}\sqrt{4-\gamma^{2}}}\lesssim 25$ K, much larger than the experimental temperature.

As increasing temperature washes out correlations, and mean-field theory neglects correlations, the high temperature PIMC results become more similar to the UHF results, as see in Fig. 4(b). This observation holds in all cases we have checked.

IV.1.2 The $S_{z}=\frac{1}{2}$ subspace of $N=3$

We performed PIMC on both sides of the $(\frac{1}{2},C_{\infty})\to(\frac{1}{2},C_{s})$ transition. For moderate magnetic fields the correlations characteristic for the $C_{s}$ state — identical spin at large angle about the origin, opposite spins at moderate angles, testifying an isosceles triangle configuration — already appear at small coupling, in the $C_{\infty}$ range; see Fig. 5(a) for an example. These correlations intensify without any sudden transformation when the UHF ground state changes to $(\frac{1}{2},C_{s})$; cf. Fig. 5(b). On the other hand, at larger fields $\gamma=1.4$ and 1.6 we find that the minority spin electron is localized off the more tightly bound majority spin molecule, and tends to be on the same line as the majority spin electrons, see the $\theta=0$ and $\theta=\pi$ peaks of the different-spin ACF in Fig. 5(c). In these cases, we have not been able to reach the asymptotic behavior at the highest inverse temperature $\beta^{\ast}$.

IV.2 Quantum dot beryllium, $N=4$

The relevant UHF ground-state structures are the rotationally invariant $(0,C_{\infty})$, $(1,C_{\infty})$, $(2,C_{\infty})$ states, and the symmetry-broken states $(0,C_{2v})$, $(1,C_{s})$, $(2,C_{4v})$, whose particle and spin density configurations are sketched in Fig. 6. A fourth symmetry-breaking structure $(0,C_{s})$ appears in PIMC only.

For $B>0$, the weak-coupling limit has $S_{z}=0$, as the degenerate shell of Fock-Darwin states is destroyed by any magnetic field, and Hund’s rule no longer selects the ground state. For relatively weak fields, $\gamma\leq 0.6$, the subsequent transitions follow the path $(0,C_{\infty})\to(0,C_{2v})\to(2,C_{\infty})\to(2,C_{4v})$. The $S_{z}=1$ states play no role here. For $\gamma\geq 0.8$, there is an intervening $S_{z}=1$ state between the unpolarized and the fully polarized states. For $\gamma=0.8$ and 1, this is symmetry-broken, while for $\gamma\geq 1.2$, it is rotationally invariant. The interval where $(0,C_{2v})$ is the ground state shrinks with increasing $\gamma$ and completely disappears for $\gamma\geq 1.2$. See Fig. 6.

IV.2.1 The $S_{z}=2$ subspace of $N=4$

We performed PIMC simulations near the coupling strength where UHF predicts the last spin flip, in order to check the UHF prediction that the electron system recovers rotational invariance in this range. As Fig. 7 demonstrates, the angular correlations characteristic of a $C_{4v}$ type are clearly present at the lowest coupling strength where UHF predicts full spin-polarization. At higher $\lambda$’s these correlations increase slightly (not shown). For $\gamma=0.2$ and $0.4$, the temperature does not affect the result, indicating that it is below the energy scale of ordering, while for $\gamma\geq 0.6$ this is not the case. Again, no sudden change is seen in the correlation when the UHF ground state changes to $(2,C_{4v})$. Similar correlations are present for somewhat higher magnetic fields. For $\gamma=1.6$, however, the picture changes. At small couplings we see peaks in the angular correlation near $\pi/2$ and $\pi$, but the latter peak is definitely sharper; cf. Fig. 8(a). The particle density in real space [Fig. 8(b)] shows a unidirectional modulation, but we note that such a picture is unavoidably washed out in longer simulations because of the orientational random walk of the system. At larger couplings, Figs. 8(c) and 8(d), we observe an elongated system with two clearly separated localized peaks. This behavior is also present in the angular correlation, where the peak $\theta=\pi/2$ splits into two nearby peaks. This gradual reordering does not coincide with the $(2,C_{\infty})\to(2,C_{4v})$ transition of the UHF ground state; it occurs at higher $\lambda$.

IV.2.2 The $S_{z}=0$ subspace of $N=4$

Here we focus on the small-$\lambda$ region, where UHF predicts a transition from $(0,C_{\infty})$ to $(0,C_{2v})$. In most cases we observe a gradual build-up of different-spin correlations near $\theta=\pi/2$, a feature of the $C_{2v}$-type spin modulation. No sudden change occurs at the $\lambda$-values where UHF predicts the $C_{\infty}$ to $C_{2v}$ transition, c.f. the first two panels in Fig. 9. Notice that for $\lambda=0.6$, the same-spin correlation has a step-like fine structure. We cannot offer an interpretation of these steps, but note that they are very well reproducible and temperature-independent as seen in Fig. 9(a).

A very interesting situation arises at $\gamma=0.6$. For small $\lambda$, we see ring-like structures with slightly increased correlations with respect to the UHF state. For $\lambda\geq 1.2$, within the $(0,C_{2v})$ region of UHF, we see a correlation of like spins at $\theta\approx\pi/2$, which corresponds to a dipolar spin modulation, not a quadrupolar one of $C_{2v}$ symmetry. See Fig. 9(c). Schematically, this state would correspond to the $(0,C_{s})$ configuration in Fig. 6, slightly distorted, a molecule that never occurs as a UHF ground state.

IV.2.3 The $S_{z}=1$ subspace of $N=4$

This subspace is relevant in UHF in small intervals for $\gamma\geq 0.8$. No symmetry-breaking is predicted by UHF for $\gamma\geq 1.2$. PIMC is performed at the smallest $\lambda$ where the UHF ground state has $S_{z}=1$. In all cases we find that at sufficiently low temperature the minority spin electron increasingly avoids the cloud of majority spin electrons, cf. Fig. 10. At $\gamma=1.2$, $\lambda=0.8$ and $\gamma=1.6$ and $\lambda=0.4$, the majority spins prefer opposite angles. In real space histograms this appears as a striped pattern in the majority spin cloud, cf. Figs. 10(c) and 10(d). The strongly localized minority spin electron is located in an interstitial. At $\gamma=1.2$, $\lambda=0.8$, there are two such locations, as the peak near $\pi/2$ in the different-spin angular correlation demonstrates. At $\gamma=1.6$, $\lambda=0.4$ there are four such locations, the different-spin angular correlation having peaks near $\theta=\pi/4$ and $3\pi/4$.

IV.3 Quantum dot boron, $N=5$

Apart from the three rotationally invariant states $(S_{z},C_{\infty})$ (where $S_{z}=\frac{1}{2}$, $\frac{3}{2}$, or $\frac{5}{2}$), the symmetry-breaking states that occur as a mean-field ground state or a correlated state are sketched in Fig. 11.

As seen in Fig. 11, UHF predicts different behavior for small magnetic fields $0<\gamma\leq 0.8$ and larger ones $\gamma\geq 1$. In the former case the spatial symmetry breaks early, then full polarization is reached through a sequence of spin flips, which restore rotational symmetry before eventually giving way to a spin-polarized symmetry-broken state. In the latter case, spin is polarized first, and the only symmetry-broken state is fully polarized. Note the exceptional nature of the strong coupling limit at $\gamma=0.2$: here the symmetry of the particle density becomes $C_{4v}$ instead of $C_{5v}$.

IV.3.1 The $S_{z}=\frac{5}{2}$ subspace of $N=5$

Figure 12 shows PIMC simulations in the fully polarized sector on both sides of the symmetry-breaking transition as predicted by UHF. Again, we find that the angular correlation evolves continuously on both sides; the angular correlation is already strong at the $\lambda$ where the last spin flips. Most notably, PIMC confirms the UHF prediction that the electrons are correlated in a centered square structure at $\gamma=0.2$.

IV.3.2 The $S_{z}=\frac{3}{2}$ subspace of $N=5$

We focus on the $(\frac{3}{2},C_{\infty})$ UHF region. As Fig. 13(a) demonstrates, ordering into pentagonal structure is clear at $\lambda=1.8$ for a weak magnetic field $\gamma=0.4$: the center of the dot is depleted, while there are correlation peaks near $\theta=2\pi/5$ and $\theta=4\pi/5$ both for identical and different spins. This state, which we call $(\frac{3}{2},C_{s})$, never occurs as an UHF ground state. The same structure emerges at higher couplings (tested at $\lambda=2.3$). At a lower coupling, $\lambda=1.3$, the spin-down electron is localized at a somewhat greater radius than the ring of majority spin electrons; the same happens at $\gamma=0.8$ and $\lambda=1.5$, 1.8. On the other hand, at $\gamma=1$, the spin-down electron is near the origin, surrounded by a square of spin-up electrons; cf. Fig. 13(b). We will call this state $(\frac{3}{2},C_{4v})$. Note that the down-spin electron actually occupies a small ring around the center of the trap. The spin-up square and the down-spin electron find it advantageous to revolve around the center of mass in a correlated manner; a situation we see quite often. Increasing the magnetic field further, the spin-down ring widens and approaches the ring of spin-up electrons, as seen in Fig. 13(c). At the same time, the angular correlation peaks of spin-up electrons move to slightly smaller angles. We interpret this scenario as a gradual distortion of $(\frac{3}{2},C_{4v})$ into $(\frac{3}{2},C_{s})$. Finally, at $\gamma=1.6$, we see in Fig. 13(d) that the radius of the spin-down ring becomes larger than that of the up-spins; the latter no longer shows strong angular ordering.

IV.3.3 The $S_{z}=\frac{1}{2}$ subspace of $N=5$

Here PIMC demonstrates the emergence of correlated structures that have no precursor in the UHF calculation. For a small magnetic field $\gamma=0.2$, we find that a double ring structure emerges, see the inset of Fig. 14(a). For majority spins, we observe a weak avoidance, only slightly stronger than by UHF, but the angular correlation peak is near $\theta=\pi/2$. The minority spin angular correlation is also peaked at $\theta=\pi/2$. We identify this as the new structure $(\frac{1}{2},C_{s}^{\prime\prime\prime})$ in the lower row of Fig. 11. A similar effect is observed at $\gamma=0.4$, c.f. Fig. 14(b) for $\lambda=1.2$. Here, the minority spins are correlated at a greater angle, suggesting a distorted structure. The central spin-up electron is slightly off the origin.

Quite different correlations arise at a greater field, $\gamma=0.6$, as shown in Fig. 14(c). The center of the trap is depleted, and only a single ring is manifest. The different-spin angular correlation is peaked near $\theta=2\pi/5$, and a minor peak is discernible near $\theta=4\pi/5$. The same-spin angular correlation also shows a double-peak structure, but now the dominant peak is at the greater, the minor peak at the smaller angle. This is consistent with the structure of the $(\frac{1}{2},C_{s})$ state, now encoded in the correlations. Note that the pentagon of localized electrons does not need to be regular, as the presence of two spins unavoidably reduces the symmetry. These findings at $\gamma=0.6$ are at odds with the UHF prediction of a $(\frac{1}{2},C_{2v})$ state.

At a still greater field, $\gamma=0.8$, we find yet another behavior. As Fig. 14(d) demonstrates, the total charge density has a dip at the center, the minority spins are strongly correlated near $\theta=2\pi/5$ and the majority ones near multiples of $\theta=2\pi/5$. We consider this as a modification of the $(\frac{1}{2},C_{s}^{\prime\prime})$ structure, now encoded in correlations. The different sharpness of the correlation peaks for the two spins is a surprise to us.

IV.4 Quantum dot carbon, $N=6$

Apart from the four rotationally invariant states $(S_{z},C_{\infty})$ (where $S_{z}=0,1,2,3$), the relevant symmetry-breaking states are sketched in Fig. 15.

Sections of the UHF phase diagram at fixed magnetic field $\gamma$ are also shown in Fig. 15. We find several remarkable features which should be critically assessed. At $\gamma=0.6$, rotational symmetry is recovered after the spin changes from 0 to 2 at $\lambda=1.7$; the same happens at $\gamma=1$ and $\lambda\geq 1.3$ and at $\gamma=1.2$ and $\lambda\geq 1$. At $\gamma=0.4$, spatial ordering is predicted with a subsequent jump in the total spin; but this evolution is interrupted by a rotationally symmetric $S_{z}=1$ phase at $\lambda=0.9$ and 1. Most remarkably, the strong coupling limit is not fully polarized in the UHF approximation at $\gamma=0.8$ and 1.4. At these magnetic field strengths, the fully polarized state has sixfold symmetry with no density peak at the center of the trap. For $\gamma\geq 1$, gradual spin polarization precedes spatial ordering.

IV.4.1 The $S_{z}=3$ subspace of $N=6$

To determine the spatial structure in the large $\lambda$ limit, we performed PIMC simulations for $\gamma=0.8$ and 1.4. The results are shown in Fig. 16. At $\gamma=0.8$, PIMC demonstrates that the hexagonal ring configuration gives way to a pentagonal structure $(3,C_{5v})$, although the peaks in the ACF are not exactly at multiples of $2\pi/5$. By contrast, at $\gamma=1.4$ no such transition is observed, and the symmetry remains sixfold up to the highest couplings studied. The radial density shows two nearby rings, cf. the inset of Fig. 16(c,d); the two radii approach each other with increasing $\lambda$. We interpret this as a distortion of the hexagon in the last item of Fig. 15; the symmetry is $C_{3v}$, but the strong coupling limit approaches $C_{6v}$. It would be interesting to know if these optimized structures make the large-$\lambda$ limit fully polarized; unfortunately, our current method is unable to decide this issue.

IV.4.2 The $S_{z}=0$ subspace of $N=6$

In this sector PIMC predicts remarkable deviations from the UHF structure. For small fields, $\gamma=0.2$ and 0.4 shown in Figs. 17(a) and 17(b), a multiple ring structure gradually emerges in the density with increasing $\lambda$. The radius of the greater spin-down ring is approximately the same as the radius of the single spin-up ring, but there is a smaller spin-down ring as well. The same-spin angular correlation suggests triangular ordering, and the different-spin correlations are peaked near $\theta=\pi/3$. Apparently, this is a symmetry-reduced version of the $(0,C_{3v})$ state, with the spin-down triangle being isosceles, and possibly reduced in size. As these may involve some loss of exchange energy, UHF never favors such a low-symmetry structure, but PIMC appreciates the energy gain from better correlations. At a higher field $\gamma=0.8$, we find cases where the correlations are characteristic of a $(0,C_{s})$ state, as seen from the density peak at the center, the almost opposite angle of spin-up electrons, and the correlation of spin-down electrons at two angles defining an isosceles triangle; different spins are correlated near $\theta\approx 2\pi/5$, see Fig. 17(c). On the other hand, still at $\gamma=0.8$ but at higher $\lambda$, PIMC detects a $C_{3v}$ structure where the UHF ground state was $(0,C_{s})$, see Fig. 17(d). At intermediate coupling values, the behavior is quite mixed, precluding a transparent interpretation.

IV.4.3 The $S_{z}=1$ subspace of $N=6$

Here PIMC finds correlated structures where UHF predicts $(1,C_{\infty})$ in the $0.4\leq\gamma\leq 1.2$ regime. The four up-spin electrons are located in a square configuration, as seen in the angular correlation function and the spin-resolved radial density in Fig. 18. One of the two down-spin electrons is located at the center of the trap, the other is off-center, avoiding the vertices of the up-spin square structure, as seen in the peaks at $\theta=\pi/4$ and $3\pi/4$ in the different spin angular correlation function. At higher $\gamma$, there is an additional feature: the down-spin electrons tend to occur at opposite sides of the trap, and, at the same time, the density has a dip at the center.

IV.4.4 The $S_{z}=2$ subspace of $N=6$

In this sector, whose $C_{\infty}$, $C_{s}$, and $C_{5v}$ states are relevant in UHF in intervals for the magnetic fields we studied, PIMC often and unequivocally overrides the UHF structures. At $\gamma=0.2$ the correlations of the $C_{s}$ state are present as a weak modulation of the angular correlation. At $\gamma=0.4$, however, PIMC shows a transition from the $C_{s}$ to the $C_{5v}$ state: as Fig. 19(a) demonstrates for $\lambda=1.9$, the same-spin angular correlation shows a pentagonal structure but the single down-spin electron is not correlated with the majority spins, and it is at the center of the dot in the spin density (inset). This applies in the whole interval of $S_{z}=2$ at $\gamma=0.4$, although for higher $\lambda$’s only at the lowest temperature we studied. At $\gamma=0.6$, on the other hand, a ring structure emerges in the density, and both the same-spin and the different-spin angular correlations show peaks at multiples of $\theta=\pi/3$ [see Fig. 19(b)]; this overrides both the $C_{\infty}$ and $C_{5v}$ predictions of UHF. We identify this state as $(2,C_{s}^{\prime})$, electrons localized at the vertices of a hexagon, one of them with a down spin. On the other hand, PIMC confirms the $(3,C_{5v})$ structure in the high-$\lambda$ limit. At $\gamma=0.8$, PIMC predicts $C_{s}$ symmetry instead of $C_{5v}$, as seen from the similarity of the same-spin and the different-spin angular correlations, c.f. Fig. 19(c). In the intervals of the $(2,C_{\infty})$ state at $\gamma\geq 1$ PIMC clearly shows the emergence of the correlations of the $(2,C_{s})$ structure as seen in Figs. 19(d) and 19(e).

IV.5 Quantum dot nitrogen, $N=7$

Apart from the four rotationally invariant states $(S_{z},C_{\infty})$ (where $S_{z}=\frac{1}{2}$, $\frac{3}{2}$, $\frac{5}{2}$, or $\frac{7}{2}$), the relevant symmetry-breaking states are sketched in Fig. 20.

A hexagon with an electron at the center is just a finite piece of a Wigner crystal, and as Fig. 20 shows, almost all symmetry-breaking UHF ground states are of this type, with different spin structures. The only exception is the $(\frac{1}{2},C_{s})$ state, which has a small interval of relevance for $\gamma=0.2$ and 0.4. For strong coupling, UHF always predicts full spin polarization, a $(\frac{7}{2},C_{6v})$ state.

IV.5.1 The $S_{z}=\frac{7}{2}$ subspace of $N=7$

Again, we can prove that the recurrence of rotational symmetry upon flipping the last spin is a mistake of UHF. As seen in Fig. 21, for $\gamma=0.2$, 1, 1.2, and 1.4 the angular correlation function indicates a hexagonal arrangement, even though for $\gamma=1.4$ this correlation is rather weak at first, but gets stronger with increasing $\lambda$.

IV.5.2 The $S_{z}=\frac{1}{2}$ subspace of $N=7$

In the least polarized sector, which is relevant for the UHF ground state at small coupling, UHF predicts a transition from the rotationally symmetric state to the $C_{s}$, $C_{s}^{\prime}$, or $C_{2v}$ molecules. PIMC, by contrast, discovers a rich variety of correlated structures. For $\gamma=0.2$ and $\lambda=1.3$, still in the $(\frac{1}{2},C_{\infty})$ range by UHF, the spin-resolved radial density finds a spin-down at the center, a ring of four spin-up electrons, and two spin-down ones at slightly greater and at slightly smaller radii than the radius of the ring; cf. Fig. 22(a). The spin-up ACF shows a very weak angular correlation, consistent with unlocalized, liquid-like behavior of the spin-up electrons in the ring. The down-spin ACF has a strong peak at $\theta=3\pi/4$ and a weaker one near $\theta=\pi/4$; this fairly rigid structure of presumably two spin-down electrons (the central electron is typically uncorrelated with the rest) in the background of unlocalized spin-ups in the ring is difficult to identify with a localized structure. Moreover, this behavior is independent of temperature in the $0.8\leq\beta^{\ast}\leq 2.4$ range, and it is quite robust: we find it with only quantitative differences for other $\lambda$’s at $\gamma=0.2$. The change of the UHF ground state from $C_{\infty}$ to $C_{s}$ at $\lambda=1.4$ has no effect. At $\gamma=0.4$ we find something else: for intermediate couplings $\lambda=0.9$ and 1.1, a square arrangement of the spin-up electrons is manifest, but the spin-down electrons form a ring at a slightly greater radius; cf. Fig. 22(b). The spin-down electrons are close to three side centers of the square, as seen form the $g_{\downarrow\downarrow}(\theta)$ and $g_{\uparrow\downarrow}(\theta)$ correlation functions. As seen in Fig. 22(c), at the upper end of the $S_{z}=\frac{1}{2}$ interval, $\lambda=1.7$, a different order emerges: a spin-down electron near the center is surrounded by a slightly distorted pentagon that contains one spin-down electron; the other spin-down electron is outside; we find similar order for $\gamma=0.8$ and $\lambda=1.7$. On the other hand, the $(\frac{1}{2},C_{s})$ structure is clearly discernible in PIMC at $\gamma=0.6$, as seen in Fig. 22(d): there is a central ring of two spin-down electrons, and the remaining spin-down electron and the four spin-up ones are arranged in a pentagon. Such correlations continue in the small-$\lambda$ region where UHF finds a rotationally symmetric ground state. At $\gamma=0.8$ and intermediate couplings $\lambda=1.2$ and 1.4 (on two sides of the $C_{\infty}\to C_{2v}$ UHF transition), there is a central spin-up electron, and the angular correlation functions show a hexagonal arrangement of electrons of alternating spin around it, as the spin-resolved ACF’s are peaked near $\theta=2\pi/3$ but the different-spin ACF is peaked at $\theta=\pi/3$. However, the radial density in Fig. 22(e) demonstrates that one of the spin-down electrons is at a slightly smaller radius than the rest.

IV.5.3 The $S_{z}=\frac{3}{2}$ subspace of $N=7$

We find remarkable correlated states where the $(\frac{3}{2},C_{\infty})$ UHF ground state is relevant. As Fig. 23 shows for $\gamma=1$ and $\lambda=1.3$, the density assumes a three-ring structure. Each of the two inner rings contains a spin-down electron. The outer ring has five spin-up electrons in a pentagonal arrangement at high temperature (not shown). As the temperature is decreased, the angle between spin-up neighbors gradually changes from about $2\pi/5$ to $\pi/3$, while the spin-down electron approaches the empty slot in the outer ring; the figure shows this latter situation. The spin-down electrons are most of the time on opposite sides of the origin. The different-spin angular correlation function shows that the spin-down electrons try to avoid the spin-up pentagon, as far as their correlated motion permits. We identify this state as a distorted variant of the $(\frac{3}{2},C_{s})$ structure. At higher magnetic fields ($\gamma=1.2$, $\lambda=1$ and $\gamma=1.4$, $\lambda=0.7$), the two down-spin rings gradually merge, and the features in the same-spin angular correlation functions weaken, suggesting a liquid-like character.

IV.5.4 The $S_{z}=\frac{5}{2}$ subspace of $N=7$

In this sector we seek for possible ordering in the rotationally invariant $(\frac{5}{2},C_{\infty})$ phase, which is relevant in UHF in small intervals for $\gamma\geq 0.6$. We observe a remarkable phenomenon at $\gamma=0.6$, cf. Figs. 24(a) and 24(b). At the beginning of the $S_{z}=\frac{5}{2}$ interval, $\lambda=1.8$, we have a slightly distorted $(\frac{5}{2},C_{s})$ structure with a spin-up electron at the center, and peaks at multiples of $\pi/3$ in the spin-up and the spin-independent angular correlation function. On the other end of the interval, at $\lambda=2.4$, the center of the dot is depleted and the angular correlation has peaks at multiples of $2\pi/7$; at $\theta=\pi$ it has a dip. We denote this structure $(\frac{5}{2},C^{\prime}_{s})$, which never occurs in UHF. In both cases the down-spin electron is located at a slightly smaller radius than the five or six up-spin electrons in the ring, respectively. At intermediate $\lambda$’s we see a gradual distortion of $(\frac{5}{2},C_{s})$ to $(\frac{5}{2},C^{\prime}_{s})$, the spin-up electrons in the ring gradually opening a slot for the central spin-up electron to occupy. At $\gamma=0.8$ and $\lambda=1.8$ we see a similar intermediate state (not shown). At $\gamma=1$ we find $(\frac{5}{2},C_{s})$ in the entire $1.5\leq\lambda\leq 2.1$ interval. At a still stronger field, $\gamma=1.2$, we observe a further distortion of the $(\frac{5}{2},C_{s})$-type correlated structure: the spin-down electron pulls closer to the center, and the five spin-up electrons in the ring assume a pentagonal coordination with angular correlation peaks at multiples of $2\pi/5$, Fig. 24(c). Finally, at $\gamma=1.4$ and $\lambda=0.8$ [Fig. 24(d)], the spin-down electron is located near the center, surrounded by a hexagon of up spins, as in the $(\frac{5}{2},C_{6v})$ state; but the central electron now occupies a small ring.

IV.6 Quantum dot oxygen, $N=8$

Apart from the five rotationally invariant states $(S_{z},C_{\infty})$ (where $S_{z}=0,1,2,3,4$), the relevant symmetry-breaking states are sketched in Fig. 25.

The UHF phase diagram is shown in Fig. 25. For small $\gamma\leq 0.6$, spatial symmetry breaks first according to UHF, and it is followed by a sequence of spin flips. Symmetry is apparently restored for short intervals, a mistake that PIMC immediately refutes (see below). For $\gamma\geq 1$, the rotationally invariant state first polarizes its spin, and then undergoes a symmetry-breaking transition later according to UHF. At the intermediate value $\gamma=0.8$, UHF predicts an intractably complicated sequence of transitions. Further, it is surprising that the large-$\lambda$ limit at $\gamma=0.2$ is not a centered heptagonal structure, but a hexagonal ring with two additional particle peaks inside, ($4,C_{2v}$).

IV.6.1 The $S_{z}=4$ subspace of $N=8$

We have decided to check the last peculiarity at $\gamma=0.2$ by PIMC. As Fig. 26 shows, the sixfold ordering is always manifest in the correlations; we have checked this in the range $2\leq\lambda\leq 5.2$. At $\lambda=5.2$, the strength of the correlation slightly depends on the temperature, but it is unclear to us what this feature signifies. The radial density shows a double peak-structure, as expected; the structure is always slightly more compact according to PIMC than according to UHF, cf. the insets in Fig. 26. The corroboration of the $(4,C_{2v})$ correlations by PIMC is significant if we recall that for $N=6$ our method does overrule the prediction of the suspicious symmetry by UHF, cf. Fig. 16.

IV.6.2 The $S_{z}=0$ subspace of $N=8$

In the $S_{z}=0$ sector, we find correlations by PIMC that are mostly independent of the UHF picture. For the small field $\gamma=0.2$, UHF predicts a sequence of $(0,C_{\infty})$, $(0,C_{4v})$, and $(0,C_{s})$ states. At the highest coupling of the unpolarized interval, $\lambda=1.7$, we find that the spin-up electrons form a ring, while two spin-down electrons are located around the same radius, one at the center, one at an intermediate radius, cf. Fig. 27(a). The same-spin correlation has a sharp peak near $\theta=\pi/3$ and a broad one near $2\pi/3$, suggesting a hexagonal ring. The down-spin correlation function has a broad peak near $\theta=2\pi/3$, indicating that the down-spin electrons in the ring are second neighbors. Qualitatively, the structure is the same for smaller $\lambda$’s, even in the $(0,C_{4v})$ domain. In the small coupling $(0,C_{\infty})$ range the radial structure of spin-down electrons is visible, as well as a Coulomb hole-type correlation in the different-spin ACF. The high-coupling state at $\gamma=0.4$ is the same as at $\gamma=0.2$, as seen in Fig. 27(b), but the small coupling behavior is different. At $\lambda=0.7$, the $C_{4v}$ order appears, with a square of the spin-up electrons inside, and four spin-down electrons rotated by $\pi/4$. Two of the four spin-down electrons are roughly at the same radius as the spin-up ones, two are farther off. As $\lambda$ is increased, spin-down electrons approach the center, one of them is inside at $\lambda=1.6$ (together with a spin-up electron) [Fig. 27(c)], two of them inside at $\lambda=2.4$ (but the spin-up electron leaves). The angles between the remaining six electrons get distorted gradually to make the above mentioned hexagon at $\lambda=2.4$. At greater field, $\gamma=0.6$ and 0.8, the weak-coupling behavior is similar, but for high coupling $\lambda=1.7$, the two spin-down electrons within the hexagon are equivalent, forming a small ring with vanishing electron density at the center. See cf. Fig. 27(d).

IV.6.3 The $S_{z}=1$ subspace of $N=8$

We find puzzling results here. For $\gamma=0.2$, UHF identifies a wide interval of the $(1,C_{s})$ symmetry breaking state. Using this for phase fixing, PIMC at $\lambda=1.8$, 2.3, and 2.8 does not fully confirm such an order. For the smaller two couplings the expected peaks at multiples of $2\pi/7$ are not present in either the same spin, the different spin, the spin-up or the spin-down angular correlation function. While a first peak can be crudely identified, the function flattens out for larger angles, precluding any interpretation. At $\lambda=2.8$, the upper limit of the $(1,C_{s})$ interval in UHF, a second peak is located near $\theta\approx 4\pi/7$, but the strength of the peaks does not correspond to the count of first, second, etc. neighbors in the ring, with the appropriate spin; it is also unclear if there is a third peak near $\theta\approx 6\pi/7$. This is seemingly one of the few cases where PIMC finds a less ordered state than mean-field theory. We leave the analysis of this regime to future studies.

We obtain clearer results for larger fields. At $\gamma=0.8$ and $\lambda=0.5$, we can identify a pentagonal ring of spin-up electrons; one of the down-spin electrons is near the origin, two are outside of the ring, typically on opposite sides ($\theta\approx\pi$), see Fig. 28(a). At $\lambda=1.7$, where UHF predicts $(1,C_{s}^{\prime\prime\prime})$, PIMC identifies a quite different molecular ordering: the angular correlation function of spin-up electrons has peaks near multiples of $\pi/3$, and the same is true for the different-spin correlation function. There is a hexagonal ring of five spin-up and a spin-down electron. The two remaining spin-down electrons are inside the ring; they are correlated with each other and the electrons in the ring in a complex manner, cf. Fig. 28(b). This structure also emerges at $\gamma=1.2$ and $\lambda=1$, where UHF predicts $(1,C_{\infty})$ (not shown). Yet another behavior is seen at $\gamma=1.4$ and $\lambda=0.7$: the spin-down ACF lets us identify an isosceles triangle, but the majority electrons do not show azimuthal ordering; the different spin correlation function only testifies a weak avoidance; cf. Fig. 28(c). Finally, we have a puzzle at $\gamma=1$ and $\lambda=1.2$: for $\beta^{\ast}\leq 2.8$, we find a hexagonal ring as in Fig. 28(b), but for $\beta^{\ast}=3.2$ the system switches to a pentagonal arrangement as in Fig. 28(a). The change in the radial density is negligible during this angular rearrangement.

IV.6.4 The $S_{z}=2$ subspace of $N=8$

While UHF identifies just narrow intervals of featureless $(2,C_{\infty})$ states in the $S_{z}=2$ sector, PIMC suggests a rich variety of correlated phases. At $\gamma=0.6$ and $\lambda=1.8$, both the same-spin and the spin-up angular correlation function shows peaks at multiples of $2\pi/7$, suggesting a seven-member ring around an electron near the origin; c.f. Fig. 29(a). The radial density confirms the existence of a spin-up ring, but the spin-down electrons are located at slightly greater and slightly smaller radii, as shown in the inset. The spin-down ACF has comparable peaks around $\theta\approx 4\pi/7$ and $6\pi/7$, which indicates that the location of spin-down electrons at these angles is close to degenerate. Further, at $\gamma=1$ and $\lambda=1.5$ the radial density informs us that the central electron has spin down, c.f. Fig. 29(b) Each of $g_{\text{same}}(\theta)$, $g_{\uparrow\downarrow}(\theta)$ and $g_{\uparrow\uparrow}(\theta)$ has peaks at multiples of $2\pi/7$. $g_{\downarrow\downarrow}(\theta)$, on the other hand, indicates that the central electron is slightly off the origin, preferring locations opposite to that of the down-spin electron in the ring. We denote this state by $(2,C_{s})$. A different picture emerges at $\gamma=1.2$, $\lambda=1.1$: a hexagonal ring of up-spin electrons surrounds a smaller ring of two down-spin electrons, c.f. Fig. 29(c). The spin-down angular correlation has peaks near $\theta=\pi/2$ and $5\pi/6$, suggesting that the spin-down electrons are located as in $(2,C_{s}^{\prime})$ in the lower row of Fig. 25. The internal structure of the down-spin ring cannot be resolved with certainty; the preferential angles depend on temperature and no clear tendency is discernible. The same behavior is present at $\gamma=1.4$, $\lambda=0.8$, although the localization of the down-spins is not resolvable (now shown). At the strongest magnetic field we study, $\gamma=1.6$, $\lambda=0.5$, all we can claim is that down-spin electrons are outside of a cloud of up-spin electrons.

IV.6.5 The $S_{z}=3$ subspace of $N=8$

At a moderate magnetic field $\gamma=0.8$, UHF predicts a transition from a featureless $(3,C_{\infty})$ state to a sevenfold symmetric $(3,C_{7v})$ ground state at $\lambda=2.3$. By PIMC, however, we find that the heptagonal ring is already present at $\lambda=1.8$, the lower limit of the $S_{z}=3$ range at $\gamma=0.8$, cf. Fig. 30(a). The central spin-down electron shows hardly any angular correlation with the spin-up ones in the ring, but it is displaced slightly from the center. This behavior is already present at $\beta^{\ast}=1.6$ and hardly changes with decreasing temperature. At $\lambda=2.3$, where the UHF ground state is $(3,C_{7v})$, PIMC confirms the same state (not shown), but only at a low temperature $\beta^{\ast}=3.2$. We see, however, a deviation from the $C_{7v}$ Wigner molecule picture at stronger couplings. As Figs. 30(b) and 30(c) demonstrate for $\lambda=2.9$, one spin-up electron enters the ring, which starts to reorganize into a hexagonal structure.

Very interesting structures emerge at greater fields. At $\gamma=1$, $\lambda=1.7$, we find a spin-up electron near the center, surrounded by a heptagon, cf. Fig. 30(d). The spin-down electron is part of the latter but it is located at a slightly smaller distance from the center. We denote this state $(3,C_{s})$. For a stronger coupling $\lambda=2.4$, the spin-down electron and the electron near the center approach each other, while the ring assumes a hexagonal structure, as seen in Figs. 30(e) and 30(f). Similar order occurs at a higher field $\gamma=1.2$ and $\lambda=1.3$; the spin-down electron moves further inward and it is less correlated with the electrons in the hexagon. Further, at $\gamma=1.4$ and $\lambda=0.9$, the two electrons of different spins occupy two sites symmetrically around the origin, surrounded by a hexagon.

V Conclusion

We have applied the path-integral Monte Carlo method with phase fixing to explore the correlated phases in few-electron quantum dots in a magnetic field. We have found that the method is flexible enough to yield structures that strongly differ, even qualitatively, from the mean-field state we use for phase fixing. It is also applicable in moderate magnetic fields where the lowest Landau level approximation, a frequently used simplifying assumption both in configuration interaction and in variational studies, does not hold. Both temperature and correlation effects are fully taken into account; the most important limitation of the method is that the simulation is performed at fixed $z$ component of the total spin, and a comparison of different spin sectors is challenging. Nevertheless, there are suggestions in the literature for obtaining free energies from PIMC simulations, and this is a possible direction of further studies.

We have found that electrons in the quantum dot very often show correlations characteristic of Wigner molecules at moderate couplings where the (unrestricted) Hartee-Fock method still predicts an unbroken symmetry. Moreover, we have found cases where PIMC completely revises the symmetry-breaking structures found in the mean-field approximation, as seen, for example, in Figs. 5, 8(b), 10, and 16(a). Further, PIMC identifies correlated structures that have no precedent in mean-field theory, e.g., the $(\frac{3}{2},C_{s})$, $(\frac{3}{2},C_{4v})$, and $(\frac{1}{2},C_{s}^{\prime\prime\prime})$ phases of quantum dot boron ($N=5$), the $(2,C_{s}^{\prime})$ state of quantum dot carbon ($N=6$), and the $(2,C_{s}^{\prime})$ and $(2,C_{s})$ states of quantum dot oxygen ($N=8$), and a few others. We have presented finite temperature results, but there are cases where the correlation functions apparently no longer change with decreasing temperature, letting us infer with some degree of confidence the ground-state properties of the system.

Experimental connections are no doubt remote at the moment. This is mainly due to two reasons: the difficulty of accessing the moderate to strongly coupled range, and the scarcity of experimental probes. Using GaAs parameters, $\hbar\omega_{0}=3$ meV, and disregarding screening mechanisms, the coupling constant in current experiments can be estimated as $\lambda<1.9$. The strongly coupled regime might be approached by using a weaker confinement potential, but the larger size of the dot also requires cleaner samples. The coupling parameter can also be increased by using a semiconductor with greater effective mass, e.g., Silicon or holes in GaAs. Transport experiments detect the structural changes in the dot by the modulation of the addition energies, a quantity we cannot obtain with certainty from a path-integral Monte Carlo calculation at the moment. The development of estimators connected to optical properties is an option that we delegate to further studies.

Appendix A The derivation of the cumulant action in an external magnetic field

Here we provide computationally tractable expressions for the cumulant action in Eqs. (11)-(13), both the Coulomb interaction and the harmonic confinement potential.

For the harmonic confinement potential $V(R)=\frac{1}{2}m^{\ast}\omega_{0}^{2}\sum_{i=1}^{N}r_{i}^{2}$ the program sketched in Eqs. (11)-(13) can be performed analytically. Focusing on a single particle $i$, a straightforward calculation yields

[TABLE]

where $u=e^{-\hbar\omega_{c}\tau}$, $\ell_{c}=\sqrt{\frac{\hbar}{eB}}$ is the magnetic length, and

[TABLE]

For the pair interaction, the cumulant action becomes a sum over all pairs:

[TABLE]

where $u_{\text{C}}^{i,j}\equiv u_{\text{C}}(\mathbf{r}_{0,i},\mathbf{r}_{0,j},\mathbf{r}_{1,i},\mathbf{r}_{1,j},\tau)$ is

[TABLE]

with the complex weight

[TABLE]

Here, $\rho_{0}$ is taken from Eq. (10). Let us transform to center-of-mass and relative coordinates $\mathbf{R}=(\mathbf{r}_{i}+\mathbf{r}_{j})/2$ and $\mathbf{r}=\mathbf{r}_{i}-\mathbf{r}_{j}$ for the initial, the final, and the running variables. Further, let us insert $\int\frac{d^{2}k}{(2\pi)^{2}}e^{-i\mathbf{k}\cdot\mathbf{r}}V(k)$ for the pair interaction. The complex weight separates as

[TABLE]

and the center-of-mass integral $\int d^{2}R$ can be performed. Then

[TABLE]

The integral on the relative coordinate can also be performed. Using the notation $u_{0}=e^{-\hbar\omega_{c}\tau}$, $u_{1}=e^{-\hbar\omega_{c}t}$, $u_{2}=e^{-\hbar\omega_{c}(\tau-t)}$, and

[TABLE]

we obtain

[TABLE]

Let us change to polar coordinates as shown in Fig. 31, i.e., measuring the angle of $\mathbf{k}$ from the direction specified by $\mathbf{r}_{t}$. Then

[TABLE]

In order to evaluate Eq. (28), we have to make another approximation. Let us expand the exponent in the integrand to first order in $B$, i.e., $e^{f(B)}\approx e^{(f(B=0))}(1+Bf^{\prime}(B=0))$. Physically, this approximation can be interpreted in terms of a comparison of the de Broglie thermal wave length $\ell_{\tau}=\sqrt{\hbar^{2}\tau/m}$ and the magnetic length:

[TABLE]

which becomes $\hbar\omega\tau\gamma\ll 1$ in terms of our dimensionless parameters. As we fix $\hbar\omega\tau=0.01$ and $\gamma\leq 1.6$, Eq. (29) is fulfilled in our simulations. Substitute $t=\lambda\tau$, $q=k\ell_{\tau}$, $\tilde{r}=r/\ell_{\tau}$, $\Delta\tilde{r}=\Delta r/\ell_{\tau}$, and $\mathbf{\tilde{r}}_{\lambda}=\mathbf{\tilde{r}}_{0}+\lambda\Delta\mathbf{\tilde{r}}$, elementary algebra yields

[TABLE]

Performing the $q$ and the $\phi$ integrations, we obtain that the real and the imaginary parts are zeroth-order and first-order in the magnetic field $B$, respectively. Explicitly,

[TABLE]

and

[TABLE]

where $I_{n}$ are modified Bessel functions of the first kind, and $X=\frac{\tilde{r}_{\lambda}^{2}}{8\lambda(1-\lambda)}$.

Equations (31) and (32) still need to be evaluated numerically. Introducing the variables

[TABLE]

we can factor out scales as

[TABLE]

We then tabulate $T_{0}(y,s)$ and $T_{1}(y,s,x)$ in the interval $0\leq y,s\leq\tilde{r}_{\text{max}}$ and $0\leq x\leq\tilde{r}_{\text{max}}/2$. We also tabulate first derivatives of $T_{0}$ and $T_{1}$, which are necessary for the thermodynamic estimator $\frac{du_{\text{C}}^{i,j}}{d\tau}$. Note that $T_{0}$ has only two variables because it knows nothing about the magnetic field, and for zero magnetic field the Runge-Lenz vector is also conserved. We note the asymptotic behavior in the $y\to 0$ and the $x\to s$ limits

[TABLE]

Because of these divergences, we extrapolate beyond the last grid point where $\frac{\partial T_{1}}{\partial y}$, $\frac{\partial T_{1}}{\partial s}$ and $\frac{\partial T_{1}}{\partial x}$ is finite. $T_{0}$ and its derivatives are regular everywhere. We choose $\tilde{r}_{\text{max}}=6$ and use linear interpolation in the tabulated domain; outside of this domain we use the semiclassical approximation to the action,

[TABLE]

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). This work benefited from the TopMat program of the PSI2 project funded by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02. We are grateful to the HPC facility at the Budapest University of Technology and Economics. T. H. G. acknowledges support from the “Quantum Computing and Quantum Technologies” PhD School of the University of Basel.