Abnormal superfluid fraction and structural properties of electrons in 2D and 3D quantum dots: an ab initio path-integral Monte Carlo study

We present extensive new direct path-integral Monte Carlo results for electrons in quantum dots in two and three dimensions. This allows us to investigate the nonclassical rotational inertia (NCRI) of the system, and we find an abnormal negative superfluid fraction [Phys. Rev. Lett. 112, 235301 (2014)] under some conditions. In addition, we study the structural properties by computing a sophisticated center-two particle correlation function. Remarkably, we find no connection between the spatial structure and the NCRI, since the former can be nearly identical for Fermi- and Bose-statistics for parameters where the superfluid fraction is diverging towards negative infinity.


Introduction
The study of fermionic many-body systems at finite temperature constitutes a highly active field of research, offering a plethora of interesting physical effects. Prominent examples include ultracold atoms [1,2] such as 3 He [3,4,5,6,7], which undergoes a pairing-induced superfluid phase transition at low temperature [8] and exhibits an intriguing, non-trivial phonon-roton-maxon dispersion relation in the spectrum of collective excitations [5,6,7]. A second relevant example is given by exotic warm dense matter [9,10,11,12,13], which naturally occurs in astrophysical objects like brown dwarfs [14,15] and giant planet interiors [16,17], and has been predicted to manifest on the pathway towards inertial confinement fusion [18]. A particularly interesting system is given by electrons in quantum dots [19], which are studied in the present work. First and foremost, we mention that these exhibit a number of interesting physics, like the quantum breathing mode response to a monopole excitation [20,21], and the formation of Wigner molecules [22] and even crystallization [23] at strong coupling. A further important research topic regarding electrons in quantum dots is the investigation of addition energies [24,25,26], which can be directly compared to experimental measurements [27]. Moreover, we mention the formation of vortices and their dependence on the respective spatial structure [28,29], and the rich interplay of the rotational inertia with a strong external magnetic field [30].
From a theoretical perspective, the accurate description of all aforementioned systems requires to simultaneously take into account i) coupling effects due to the interactions, ii) thermal excitations as a result of the finite temperature, and iii) quantum degeneracy effects like Pauli blocking since identical fermions are indistinguishable. In addition, confined few-body systems [31] additionally require to take into account the interplay between i)-iii) and iv) the external potential, which, in the present study, is modelled as a harmonic oscillator. Due to these challenges, electrons in quantum dots constitute a rigorous benchmark that is often used to test novel many-body methods, e.g., Refs. [32,33,34,35,36,37,38,39,40,41].
In the present work, we use the direct path-integral Monte Carlo (PIMC) method [42,43] without any nodal constraints [44] to simulate electrons in both 2D and 3D harmonic confinements. Therefore, our simulations are afflicted with the notorious fermion sign problem (FSP) [45,46,47], but are exact within the given Monte Carlo error bars, which allows for a rigorous treatment of the intriguing interplay between the effects i)-iv) mentioned above. From a practical perspective, the PIMC method gives us direct access to the superfluid fraction of the system [48,49], which exhibits a remarkable abnormal divergence towards negative infinity in the case of fermions for some conditions [50,51]. This effect was first reported by Yan and Blume [50] for a short-range pair potential and subsequently found by Dornheim [51] for quantumdipole systems, and has been explained in terms of the symmetry of the thermal density matrix. Here, we extend these considerations to electrons in quantum dots and find analogous behaviour. In addition, we study the structural properties of the system using a recently suggested center-two particle (C2P) correlation function [52], which efficiently takes into account the rotational symmetry of the system. This allows us to resolve the interesting interplay between quantum statistics and the Coulomb interaction between pairs of electrons, and to explore the connection between the aforementioned abnormal superfluid fraction and the static structure.
The paper is organized as follows: In Sec. 2, we introduce the relevant theoretical background including the model system (1), the utilized PIMC simulation approach (2.2), the estimation of the superfluid fraction (2.3), and finally the center-two particle correlation function (2.4). Sec. 3 is devoted to the presentation of our extensive new simulation results, which we start by shortly revisiting the noninteracting system in Sec. 3.1. Subsequently, we study both the superfluid fraction and the structural properties of electrons in quantum dots, focusing on the dependence on the temperature (3.2), on the system-size (3.3), and on the coupling strength (3.4). The paper is concluded with a concise summary and outlook in Sec. 4.

Model system
We consider N spin-polarized electrons in both two-and three-dimensional quantum dots, which we model by a simple harmonic potential [19]. This gives the Hamiltonian where we assume oscillator units, corresponding to the characteristic length l 0 = /mΩ (with Ω being the trap frequency) and energy scale E 0 = Ω. The first term corresponds to the kinetic contributionK and the last two terms to the external potential and the Coulomb interaction,V ext andŴ , respectively. Moreover, the coupling constant λ can, in principle, be tuned in experiments [53].

Path-integral Monte Carlo
We use the direct path-integral Monte Carlo method [54,55], see Ref. [42] for an extensive review article. In a nutshell, the basic idea behind PIMC is to stochastically sample the thermal density matrix in the canonical ensemble whereĤ is the Hamiltonian, β = 1/k B T the inverse temperature, and R = (r 1 , . . . , r N ) T contains the coordinates of all N particles. While a direct evaluation of Eq. (2) is not possible as different contributions to the full Hamiltonian do not commute, one makes use of a Trotter decomposition [56] into P terms that are each given at P -times the original temperature. The resulting scheme becomes exact in the limit of large P and the respective convergence has been carefully checked; see Ref. [57] for a detailed discussion of different factorization schemes. The PIMC method thus allows for the quasi-exact evaluation of thermal expectation values within the given Monte Carlo error bar, which scales as ∆A ∼ 1/ N MC with the number of Monte Carlo samples N MC . In practice, we use a canonical adaption [58] of the worm algorithm presented by Boninsegni et al. [59,60]. An additional obstacle regarding the PIMC simulations of fermions (i.e., particles obeying Fermi-Dirac statistics such as electrons or protons) is given by the antisymmetry of Eq. (2) under the exchange of particle coordinates [43]. First and foremost, negative terms cannot be interpreted as a probability distribution, and the expectation value from Eq. (3) gets modified to [47] where . . . denotes the expectation value computed from the absolute values of the density matrix, andŜ measures the respective sign. In particular, the denominator of Eq. (2) is called the average sign and constitutes a straightforward measure for the amount of cancellation within the fermionic PIMC simulation. This can easily be seen by considering the first-order term to the relative statistical uncertainty [61], which directly implies that the Monte Carlo error bar increases when the sign gets small. In fact, S exponentially decreases both with the inverse temperature β and the system size N , which can be directly translated to ∆A as where ∆f denotes the difference in the free energy density of the original ( . . . ) and the modified ( . . . ) systems. Eq. (6) is the root of the infamous fermion sign problem [47], as the exponential increase of A can only be compensated by increasing the compute time as 1/ N MC , which quickly becomes infeasible. In fact, Eq. (6) constitutes the central computational bottleneck in our simulations and limits the feasible system size to N ≤ 10.

Superfluid fraction
The superfluid fraction of a finite system is typically defined in terms of its nonclassical rotational inertial (NCRI), i.e., its response to an infinitesimal rotation. More specifically, one makes use of the Landau two-fluid model, where the total density is decomposed into a normal part reacting to the rotation, and a superfluid part that remains unaffected, n = n n + n sf .
The corresponding superfluid fraction is then defined as the ratio of the superfluid to the total density, with I and I cl denoting the actual and classical moments of inertia. Within the pathintegral picture, the former can be straightforwardly estimated as [62] γ sf = which is often being referred to as the area estimator in the respective literature. In particular, Eq. (9) depends on the area that is enclosed by the paths in the plane perpendicular to the rotational axis, Figure 1. Illustration of the center-two particle correlation function ρ 2 (r 1 , r 2 , ϑ). Due to the rotational symmetry of the Hamiltonian Eq. (1), the two-particle correlations only depend on the relative angle ϑ (denoted as α throughout this work) and the respective distances to the center of the trap. Taken from Ref. [65]. Reprinted with permission of WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim.
In addition, we note that a non-zero superfluid fraction as it is defined in Eq. (9) does not necessarily indicate the usual phenomenon of superfluidity that manifests for example as a fricitionless flow in ultracold 4 He [42]. In fact, superfluidity naturally emerges as a long-range order in the off-diagonal density matrix [64], which is not possible for the small system sizes considered in this work. Therefore, one should interpret it strictly in terms of NCRI defined above. Still, as the term superfluid fraction has been used throughout the literature [50,51,53,49], we will stick to it here as well.

Structural properties
The second central aim of the present work is the investigation of the impact of quantum statistics on the structural properties of electrons in quantum dots. To this end, we have implemented the center-two particle (C2P) correlation function [52], which has already been successfully applied to harmonically confined quantum systems in earlier works [65,66,51]. The basic idea is illustrated in Fig. 1, which shows an exemplary configuration of N = 13 particles in a 2D harmonic trap. Here the shaded red circles illustrate the quantum smearing of the electrons, distinguishing them from classical point particles. To characterize two-particle correlations, one can make use of the rotational symmetry of the problem and use the modified coordinates (distance to the center of the trap of particles I and II, and the angle between them, called ϑ in Fig. 1 and α throughout the present work) shown in Fig. 1 without any loss of information. More specifically, we estimate the corresponding two-particle density matrix ρ 2 (r I , r II , α) as a histogram in our PIMC simulation.
In addition, Thomsen and Bonitz [52] have suggested that it would be advantageous to filter out effects that are caused by the inhomogeneous density profile n(r) and not by correlation effects themselves. To this end, the actual C2P is normalized by the two-particle density matrix of a hypothetical uncorrelated system with the same radial density profile as the interacting one, and is thus given by We note that Eq. (11) does not depend on the angle α by design, as noninteracting particles are not affected by their relative orientation towards each other. While the C2P defined in Eq. (12) already constitutes a sophisticated tool for the investigation of spatial correlations, the visualization and analysis of a 3D quantity is often difficult. As a final step, we therefore define the integrated C2P which can be interpreted in the following way: Eq. (13) constitutes a measure for the probability to find one particle at distance r I to the center of the trap, with an angular difference of α to a second particle somewhere in the interval r II,min ≤ r II ≤ r II,max , where the min and max values are typically chosen as the limits of a shell in the density profile n(r). All results for the C2P shown in this work have been obtained according to Eq. (13).

Noninteracting system
Let us start the discussion of our new simulation results with a brief study of the temperature dependence of the superfluid fraction of noninteracting bosons and fermions in a 3D harmonic trap. While similar analyses have been presented in 2D in Ref. [51] and in 3D in Ref. [50], here we select somewhat different parameters. In addition, the availability of exact analytical results makes it a valuable test case to verify the implementation of our simulation scheme. Finally, the abnormal superfluid fraction of interacting electrons presented in a later section follows the same basic mechanism as the ideal case, which makes it a good starting point. In  fraction is small and almost equal for bosons and fermions for large temperatures T = β −1 , whereas differences emerge for smaller T . In fact, it always holds as the systems becomes purely classical in this limit, I → I cl . Further, the first correction to I cl is due to the finite extension λ β = √ 2πβ of quantum particles, which is similar for both bosons and fermions. With decreasing temperature, the paths of different particles in a PIMC simulation start to overlap, and the formation of permutation cycles, i.e., paths containing more than a single particle in them, becomes increasingly likely; see Ref. [43] for a recent extensive discussion. In practice, this leads to an increased impact of quantum statistics, which explains the stark qualitative and quantitative differences between the red and green data sets.
More specifically, bosons always attain a superfluid fraction of one in the limit of T = 0 at these conditions, independent of the particle number N . Yet, we stress that this does certainly not mean that ideal bosons become an actual superfluid, as the Landau criterion requires a small value of the coupling parameter λ in this case [67]. Fermions, on the other hand, exhibit a substantially more complicated behaviour. In particular, the superfluid fraction diverges towards negative infinity for both N = 2 and N = 6, and attains one for N = 4, although for lower temperatures compares to  bosons. This remarkable effect has been explained by Yan and Blume [50] in terms of the structure of the density matrix. Specifically, it is the presence of an energetically low-lying state with a finite circulation for some N that leads to a diverging moment of inertia and to the observed behaviour in γ sf . Finally, we stress the excellent agreement between theory and PIMC data for both bosons and fermions everywhere, which is a strong verification of our implementation and analysis. In addition, we note that simulations are restricted to above a minimum temperature in the case of fermions due to the fermion sign problem (see Sec. 2.2 above), as can be seen most clearly for the leftmost red data point for N = 4. In contrast, there is no such issue for bosons, where simulations are possible even at extremely low temperatures.

Abnormal superfluid fraction of electrons: dependence on temperature
Let us begin our investigation of interacting electrons by studying the two-particle Hamiltonian. The relative and center of mass coordinates are separable; the centerof-mass energy spectrum is solved analytically and the relative energy spectrum is obtained by the finite-element method. Eigen energies subject to a relative energy cutoff of around 100 Ω are computed. We calculate the actual moment of inertia I through I = 2 M 2 th / (k B T ), where M is the angular momentum projection quantum number and indicates the thermal average [50]. The classical moment of inertia I cl is obtained from the thermal average N k=1 r 2 k,⊥ . Here r 2 k,⊥ is the square of distance of the k-th particle to the z axis. Due to spherical symmetry, it is related to the trap energy  converting it to the thermal average of the derivative of the energy with respect to the trap depth through the Hellmann-Feynman theorem. The solid blue, dotted red, and dashed green lines in Fig. 3 show the numerical results for λ = 0, λ = 1, and λ = 10. First and foremost, we note that these results are in excellent agreement to our new PIMC results, which are depicted by the corresponding symbols. Secondly, we find that the coupling strength has a comparatively small impact on the behaviour of γ sf even for λ = 10; indeed, the results for moderate coupling (λ = 1) can hardly be distinguished from the ideal curve. This is a strong indication that the response of electrons in a quantum dot to an external rotation is largely determined by the properties of the corresponding noninteracting system.
Let us next consider our new PIMC data for the temperature dependence of γ sf for larger systems that is shown in Fig. 4. The left panel shows results for N = 4 particles, and the squares and circles have been obtained for Fermi-and Bose-statistics. The green curves show results for the ideal case (λ = 0) already discussed in Fig. 2 above and have been included as a reference. The dash-double-dotted grey curves have been obtained for a small yet finite value of the coupling parameter, λ = 0.5, and closely follow the ideal case. In particular, the fermionic results cannot be distinguished from the noninteracting curve with the bare eye, whereas the gap for Bose-statistics is significant.
Going to moderate coupling strength, λ = 3 (dashed red), shifts the curves substantially to the left compared to smaller λ. Heuristically, this can be understood as follows: due to the stronger Coulomb repulsion, the particles are more clearly separated. Consequently, it takes lower temperatures until paths of different particles overlap, which is the driving mechanism of a large superfluid fraction [49,53]. Again, we find that the effect of the Coulomb repulsion is more pronounced in the case of bosons compared to fermions. The reason for this effect has been reported by Dornheim [51] for the case of dipole-dipole interaction: for fermions, weak interactions are effectively masked by the Pauli exclusion principle which separates even ideal fermions from each other; in contrast, bosons tend to cluster next to each other, and are thus more strongly affected by the Coulomb repulsion. A more detailed analysis of this effect is presented in Sec. 3.4 below.
Finally, the blue curves correspond to relatively strong coupling, λ = 10. Interestingly, the curves for the two different types of quantum statistics are quite close to each other in this case, and both approach unity for similar temperatures.
Let us next consider the same property, but for N = 6 shown in the right panel of Fig. 4. For bosons, we find a quite similar picture compared to the left panel. Evidently, the particular particle number does not substantially shape the qualitative physical behaviour at these conditions. In stark contrast, all fermionic curves attain a maximum in γ sf in the range 0.5 ≤ β ≤ 5 until they diverge towards negative infinity for T → 0. Still, the impact of the coupling parameter λ is only quantitative, and the main feature is already present in the ideal case. Furthermore, we again find that fermions are less affected by the Coulomb repulsion compared to bosons.
In Fig. 5, we present a similar analysis, but for a purely 2D system. The results are very similar to the 3D case, but the T = 0 limit is flipped in the case of fermions: for N = 4, γ sf diverges towards negative infinity, and for N = 6, it attains unity just as in the Bose-case. Still, the phenomenological reason is the same for both 2D and 3D: it is the presence of a state with finite circulation [50] in the respective cases that causes the negative divergence of the superfluid fraction.
Let us briefly postpone our analysis of the connection of the abnormal superfluid fraction of fermions to the structural properties of the system to touch upon the manifestation of the fermion sign problem. In Fig. 6, we show the PIMC expectation values for the denominator in Eq. (4), i.e., the average sign. The top and bottom rows have been obtained for 2D and 3D systems, and the left and right columns correspond to N = 4 and N = 6 spin-polarized electrons. Furthermore, the different data points distinguish data for various values of the coupling parameter λ. First and foremost, we note that all depicted curves exhibit the same qualitative behaviour: for small β (large temperatures), the extension of the paths of different particles is substantially smaller than the average inter-particle separation. Consequently, permutation cycles only rarely appear in the PIMC simulation, most configuration weights are positive, and the average sign is large. With decreasing temperature, the thermal wavelength increases as λ β ∼ √ β and, eventually, the paths of all particles in the system overlap. This makes the formation of permutation cycles significantly more likely, configurations with positive and negative weights appear with a similar frequency in the PIMC simulation, and the average sign drops. The straight lines in Fig. 6 have been obtained from exponential fits according to S(β) = ae −bβ (15) above an empirical minimum value of β, with a and b being the free parameters. Evidently, the fits are in excellent agreement to the PIMC data, which confirms the exponential nature of the fermion sign problem reported in earlier studies [47]. A related interesting point is the PIMC estimation of the corresponding energies. More specifically, we consider the well-known virial theorem [68], which gives a relation between the kinetic (K), interaction (W ), and external potential (v ext ) contributions to the total energy, In practice, we directly estimate K in our PIMC simulation using the standard thermodynamic estimator (see, e.g., Ref. [42]) and compare it to the right-hand side of Eq. (16). The results are depicted in Fig. 7, where we show the relative deviation between the two estimations of K in per cent. We note that the virial theorem as it is given in Eq. (16) is a general property of the Hamiltonian and holds independent of the particular type of quantum statistics. The left panel has been obtained for Bosestatistics, and the red circles and blue diamonds show results for the 3D and 2D case, where the latter have been shifted by 0.2% for better visibility. Evidently, all data points fluctuate around zero deviation within the given level of statistical uncertainty. Further, we find that the error bars are only weakly dependent on β, and are of the order of ∆K/K ∼ 0.01%. The right panels shows results from the same simulation computed for Fermistatistics by evaluating Eq. (4). Again, we find that the data points fluctuate around zero deviation between the two different estimations of K, which is a strong verification of our simulation scheme and its implementation. Yet, we find a drastic increase of the error bars with β as predicted by Eq. (6). This nicely illustrates the relative difficulty of fermionic PIMC simulations compared to the Bose-case, where this problem is absent. In particular, simulations of as few as four electrons become unfeasible even on modern supercomputers when the temperature is decreased.
Let us next consider the structural properties of electrons in a 2D quantum dot, and their relation to the abnormal superfluid fraction reported above. To this end, we estimate the integrated C2P defined in Eq. (13) in Sec. 2.4. The results are shown in Fig. 8 for N = 4 at λ = 10 and β = 10. The top row shows the integrated C2P itself in the α-r I -plane, with the reference particle being located somewhere in the range 1.5 ≤ r II ≤ 2.5 (horizontal dashed grey lines), i.e., around the maximum of the radial density (see the bottom right panel of Fig. 8). Further, the left and right panels show results for fermions and bosons, respectively. Naturally, the C2P is fully symmetric in α around 180 degrees and we restrict ourselves to this half-range. Most importantly, we find that the type of quantum statistics has only very mild effects on the structural properties on the system at these conditions. This is a direct consequence of the strong coupling strength, which effectively separates the paths of individual particles from each other in the PIMC simulation and suppresses the formation of permutation cycles [43].
For small angular separations α, the integrated C2P exhibits a pronounced exchangecorrelation hole, followed by a maximum at α = 90, a subsequent minimum which is substantially less pronounced than the exchange-correlation hole itself, and a second maximum at the opposite end of the system, i.e., around α = 180. This can be seen particularly well in the bottom left panel of Fig. 8, where we show a scan-line over the C2P at r I = 2, see the solid red line in the top row. The red circles and green squares show results for fermions and bosons, but the deviations between the two data sets can barely be resolved with the naked eye. From a physical perspective, this indicates that the most common configuration of particles is given by a single ring around the center of the trap containing all four particles. This is further substantiated by the radial density profile n(r) shown in the bottom right panel of Fig. 8. Here the main impact of quantum statistics can be found around r = 0, where the probability to find a fermionic particle is reduced compared to the case of bosons.
Let us now connect these findings regarding the structural properties of the system to the abnormal superfluid fraction investigated above. In particular, it can be seen in Fig. 5 that γ sf strongly depends on the type of quantum statistics at these conditions. Remarkably, this is not reflected in the structural properties of the system, which, indeed, can hardly be distinguished at all.

Particle number dependence
A further interesting topic of investigation is the dependence of the properties of small electronic systems on the number of particles. In this context, we mention that the properties of confined few-particle systems often sensitively depend on N , which can be seen particularly well in the investigation of addition energies [24,25,26,27]. Furthermore, Filinov et al. [53] have found that the nonclassical rotational inertia of strongly coupled charged bosons exhibits certain "magic numbers", where the superfluid fraction is most pronounced. In this work, we provide a similar analysis for the superfluid fraction γ sf of electrons in quantum dots.
The results are shown in Fig. 9, where we show the N -dependence of γ sf both in 2D (left panel) and 3D (right panel) at moderate coupling, λ = 3. The circles and squares show results for Bose-and Fermi-statistics, and the red and green curves have been obtained for β = 5 and β = 3. For all depicted cases, the bosonic curves are quite smooth and exhibit a monotonous increase (the small exception being N = 3 for β = 3 in 2D) in the superfluid fraction with N . Heuristically, this can be understood as follows: in the harmonic trap, the average density increases with N . This, in turn, leads to an increased value of the degeneracy parameter χ = λ β /r, where r is the average inter-particle distance [49], and thus to an increasing value of γ sf .
For fermions, on the other hand, the situation is significantly more complicated and interesting. In the 2D-case, the superfluid fraction is rapidly oscillating with N , which is especially pronounced for β = 5. In particular, the negative and positive values of γ sf can be traced back to the ideal case, see Ref. [51] for such data in two dimensions. A second trend, which can be seen most clearly for β = 3, is given by the overall decrease of the superfluid fraction with N . This, too, is expected by simple physical intuition. An increase in N leads to a transition from a finite to a bulk system. For bosons, the crossover from a normal to a superfluid system then attains the character of a true phase transition under the right conditions [49]. In contrast, a single species of fermions such as spin-polarized electrons cannot exhibit the required off-diagonal long-range order [64], and a corresponding bulk system will thus not be superfluid. This is reflected in the observed decrease in γ sf with N . For completeness, we mention that superfluidity in Fermi-systems is still possible, but requires a pairing mechanism such as the interaction with phonons leading to superconductivity [69].
We find similar trends in the case of fermions for 3D systems, although both the oscillations with N and the decrease with increasing system size are somewhat less pronounced.
Let us conclude this investigation of the N -dependence by briefly revisiting the fermion sign problem. To this end, we show the average sign obtained for the simulations shown in Fig. 9 in Fig. 10. The points show the PIMC data, and the straight lines correspond to exponential fits of the form and have been obtained by introducing an empirical minimum value of N = 7 (vertical dashed grey line). Strictly speaking, Eq. (17) is only expected to hold for a bulk system such as the uniform electron gas [47,12], whereas the dependence of a finite, trapped system is expected to be more complicated. Still, we do find at least a qualitative agreement with the exponential decay. The particularly small value of S for the 3D case at β = 5 is directly reflected by the large error bars in γ sf observed in Fig. 9. This again strongly illustrates the severe limitations of fermionic PIMC simulations due to the notorious sign problem.

Dependence on coupling strength
A highly interesting question is the precise nature of the intricate interplay of the Coulomb repulsion with the effect of quantum statistics. In this regard, our computationally expensive PIMC simulations are uniquely suited to provide reliable answers as no approximations are involved. In Fig. 11, we show results for the integrated C2P [cf. Eq. (13) above] for N = 4 spin-polarized electrons in a 2D harmonic trap at a moderate temperature, β = 1. The left and right columns have been obtained for Fermi-and Bose-statistics, and the reference particle is always located in the interval 1.5 ≤ r II ≤ 2.5 (horizontal dashed grey lines in the top left panel). The top row corresponds to a very weakly coupled system (λ = 10 −3 ), where the Coulomb repulsion has a negligible effect. For fermions, the Pauli exclusion principle nevertheless results in a pronounced exchange-correlation hole around α = 0, whereas the rest of the system appears to be relatively featureless. In stark contrast, a corresponding Bose-system at the same conditions exhibits the opposite behaviour, as bosons tend to cluster around each other. Increasing the coupling strength to λ = 0.1 leaves the Fermi-system unchanged, whereas the positive exchange-correlation hill of the bosons is substantially decreased. From a physical perspective, this means that  bosons are strongly influenced even by weak repulsive interactions, whereas the latter are effectively masked by the Pauli principle in the case of fermions. For completeness, we mention that the same point has recently been reported by Dornheim [51] in the case of dipole-dipole repulsion. The same trend has also been noted in our discussion of Figs. 5 and 4 above, where the superfluid fraction γ sf was substantially more affected by the coupling parameter λ for bosons compared to fermions.
Going back to the integrated C2P shown in Fig. 11, a further increase of the coupling strength to λ = 1 (middle row) again only leads to very mild changes for Fermi-statistics, whereas there appears a qualitative change for bosons: the exchange-correlation hill is changed to a hole, although it is still significantly less pronounced compared to the left panel. At λ = 3, i.e., at moderate coupling, the importance of quantum statistics is drastically reduced and both panels give a nearly identical picture. Finally, the strong Coulomb repulsion at λ = 10 effectively separates individual particles, and no differences between bosons and fermions can be seen with the naked eye.
Let us next consider Fig. 12, where we show a scan-line of the integrated C2P at r I = 2 (solid red line in the top left panel of Fig. 11). The solid and dashed lines distinguish fermions and bosons, and the different colours correspond to values of the coupling parameter λ. The arguably most striking feature of this plot is the exchangecorrelation hill for bosons at λ ≤ 0.1. In addition, it clearly illustrates the small impact of the coupling strength on fermions for the same λ. For λ = 3, there is still a distinct impact of quantum statistics, whereas no such difference can be resolved for λ = 10. The latter case also constitutes the only curve with a second maximum around α = 180. Let us close this discussion of the integrated C2P with a practical remark. Obviously, the probability to find two identical fermions at the exact same position must always be zero, ρ 2 (r, r, 0) = 0. Yet, by integrating over a finite range of the reference particle r II,min ≤ r II ≤ r II,max , the integrated C2P can be nonzero for all values of r I even for α = 0.

Summary and Outlook
In summary, we have presented extensive new ab initio path-integral Monte Carlo simulation results for electrons in 2D and 3D quantum dots. In particular, we have investigated the nonclassical rotational inertia of such systems and have found a divergent, negative abnormal superfluid fraction which manifests under certain conditions due to the interplay of Fermi-statistics and the harmonic confinement [50]. In addition, we have thoroughly investigated the structural properties of such systems by computing the center-two particle correlation function suggested by Thomsen and Bonitz [52]. The comparison of results for Bose-and Fermi-statistics has revealed that the onset of the abnormal superfluid fraction is hardly connected to the static structure, which can be nearly independent of the type of quantum statistics, whereas the respective superfluid fractions diverge. Instead, it is a direct consequence of the structure of the density matrix with respect to the angular momentum. In addition, we have substantiated earlier findings for harmonically confined quantum dipole systems [51], where it has been found that the repulsion between two fermions is effectively masked by the Pauli exclusion principle for weak too moderate values of the coupling parameter λ. In stark contrast, bosons sensitively react even to small changes in the coupling strength.
Future extensions of our work might include the investigation of spin-unpolarized systems, which would allow for the study of a spin-resolved C2P. In particular, it is clear that the observed effective screening of the pair interaction by the Pauli principle will not apply to fermions of a different spin-orientation, which will lead to an interesting mixture of different effects. In addition, we note that the central computational bottleneck of the present study is given by the fermion sign problem, which precludes the simulation of larger systems and lower temperatures. In this regard, we mention the remarkable recent progress in the fermionic quantum Monte Carlo simulation at finite temperature in the context of warm dense matter research [70,71,32,72,73,74,75,76,77,78,79,80,81]. Specifically, the adaption of these novel methods to the simulation of electrons in quantum dots seems promising and rewarding. Finally, we mention that the C2P that has been used in the present work allows to reliably describe crossovers and melting phenomena in finite systems [52]. This is in stark contrast to previously employed Lindemann-type approaches [82] that rely on the underlying Monte Carlo sampling scheme. Thus, the present set-up can be straightforwardly used to study the impact of quantum effects onto the crystallization of trapped few-body systems.
(BMBF) and by the Saxon Ministry for Science, Culture and Tourism (SMWK) with tax funds on the basis of the budget approved by the Saxon State Parliament. The PIMC calculations have been carried out on a Bull Cluster at the Center for Information Services and High Performance Computing (ZIH) at Technische Universität Dresden, and at the Norddeutscher Verbung für Hoch-und Höchstleistungsrechnen (HLRN) under grant shp00026.