Zonal Estimators for Quasiperiodic Bosonic Many-Body Phases.

In this work, we explore the relevant methodology for the investigation of interacting systems with contact interactions, and we introduce a class of zonal estimators for path-integral Monte Carlo methods, designed to provide physical information about limited regions of inhomogeneous systems. We demonstrate the usefulness of zonal estimators by their application to a system of trapped bosons in a quasiperiodic potential in two dimensions, focusing on finite temperature properties across a wide range of values of the potential. Finally, we comment on the generalization of such estimators to local fluctuations of the particle numbers and to magnetic ordering in multi-component systems, spin systems, and systems with nonlocal interactions.

1. Introduction

Path-integral Monte Carlo (PIMC) methods [1] are of great importance for the simulation of strongly correlated systems where other techniques fail, especially in two and three spatial dimensions. Over the last thirty years, this has been amply demonstrated on quantum fluids [2,3,4] and, more recently, in ultra-cold gases like, for instance, dipolar systems [5,6,7,8,9] and Rydberg atoms [10,11,12,13]. For strongly-interacting quantum fluids, there is at present considerable interest towards the exploration of patterns owing to peculiar symmetries such as quantum-cluster crystals [14,15], stripe phases [16], or cluster quasicrystals [17,18,19], with the aim of understanding fundamental physical phenomena. In this regard, also thanks to the increase in computational capabilities, advancements in PIMC methods continue to play a key role.

In this work, we detail the numerical techniques required to investigate the quantum properties of interacting trapped bosons in an external quasiperiodic potential at finite temperature, with specific attention to superfluidity. Quasicrystals are a fascinating state of solid-state matter exhibiting behaviors halfway between a periodically ordered structure and a fully disordered system. They were first synthesized in 1982 (and their discovery later announced in 1984) by Shechtman et al. [20]. Later, Bindi et al. demonstrated that quasicrystals can also originate naturally, in the presence of extreme conditions such as collisions between asteroids [21]. Their properties have already been the subject of extensive theoretical investigation [22], motivated in part by the discovery of aperiodic tilings that can cover the plane without being bounded by the symmetries of classical crystallography, such as the Penrose tiling [23]. At finite temperature, the thermodynamic features of quasicrystals can be established in terms of the interplay between different length and energy scales pertaining to the inter-particle potentials [24]. These classical systems were found to remain stable even at zero temperature [25]. Quasicrystalline properties have been observed in a variety of physical systems, for instance, in nonlinear optics [26,27,28], on twisted bilayer graphene [29] and in ultra-cold trapped atoms [30,31]. In the latter case, quasicrystalline structures generated by means of optical lattices are employed to experimentally investigate remarkable effects such as many-body localization in one and two-dimensions [32], and have been suggested as a candidate to probe the existence of two-dimensional Bose glasses [33]. In this regard, recent PIMC simulations support the existence of a Bose glass phase, fully stable and robust at finite temperature, in a region of parameters suitable for experimental setups [34,35]. Other works have delineated zero-temperature phase diagrams, in the mean-field approximation as well as for a strong interactions using ab-initio techniques [36,37,38,39].

Here, we summarize the derivation of the pair-product approximation for particles interacting through hard-core interactions in two or three dimensions, and we present the details of our implementation. Then, we explore a new zonal estimator, which gives access to local information about the superfluidity in finite regions of trapped systems, and is therefore well-suited to the study of spatially inhomogeneous potentials. Zonal estimators can be relevant to the detection of correlated phases, such as the Bose glass phase, which is characterized by rare regions where superfluidity and finite compressibility coexist [34].

This paper is organized as follows. In Section 2, we present and discuss the PIMC methodology for ensembles of interacting bosons through the pair-product approximation. In Section 3, we introduce a model Hamiltonian describing interacting trapped bosons subjected to a quasiperiodic potential. Structural properties such as density profiles and diffraction patterns are shown in Section 4, whereas we examine global quantum features in Section 5. The zonal estimator of the superfluid fraction is explored in Section 6. To conclude, Section 7 is devoted to the discussion of our findings, drawing some conclusions.

2. Methodology

In this section, we review the implementation of the PIMC to the study of an interacting Bose gas in an external potential at finite temperature. This methodology aims to sample the partition function of a quantum system at finite temperature. In line with Feynman’s path integral theory [40,41], thermodynamic properties are addressed by considering an equivalent classical system, in which each quantum particle is represented by a classical polymer. As a result, quantum quantities, like for instance superfluidity or Bose–Einstein condensation, can be mapped across the equivalence as properties of the polymers themselves [42]. The evaluation of those quantities takes place via a standard classical Monte Carlo procedure such as the Metropolis algorithm [43], allowing us to sample thermodynamic properties within a precision limited only by numerical and statistical errors. At present, one of the most efficient ways of sampling configurations of connected polymers is operated through the so-called worm algorithm [2,44]; originally developed for the grand-canonical ensemble, we routinely use the worm algorithm in its canonical version to sample superfluid fraction, condensate fraction, or ground state energy.

In the following, we recall the formalism and derivations at the core of PIMC. The partition function, , is defined as the trace of the equilibrium density matrix operator, , at temperature T Ref. [45]: being the inverse temperature parameter, .

For N distinguishable particles, denoting with the position of the i-th particle, and introducing = , we can project the density matrix operator on the basis of spatial coordinates , obtaining

For an ensemble of bosons, taking into account permutations, we arrive at where = denotes a permutation of the particle coordinates.

Introducing a decomposition of the density matrix operator into a convolution of density matrices at a higher temperature, Equation (4) yields with breaking up into M smaller intervals , and the coordinates of particles on a given time slice. To each particle corresponds, then, a classical polymer made of beads, connected with each other through harmonic springs [45]. Errors introduced by the equivalence are reduced as M increases. Likewise, the same decomposition can be applied to the evaluation of observables by Monte Carlo sampling.

For a generic diagonal observable, such that , it follows that

Equation (5) is evaluated through a stochastic process, consisting of the generation of random configurations from the probability distribution

The thermodynamic average then is measured as an average of over the sampled configurations [46].

Having to employ a finite number of time slices, M, the most sensitive step of the procedure lies in finding a good approximation of the high-temperature density-matrix elements in (4) [46,47,48]. Due to the nature of the two-body interaction potential between the bosons in Hamiltonian (see (37) for an application), and to the density regime of interest, in the present work we apply a pair-product approximation (PPA) ansatz [45]. The rest of this section is devoted to the treatment and implementation of contact interactions in this context; similar derivations and other details can be found, e.g., in [48,49] and references therein, and the supplemental material of [37].

We express the density-matrix terms as (to keep Equation (7) simple, we omit the indices m). Here, with , is the density matrix of the non-interacting Hamiltonian of N particles, .

For two particles, labeled i and j, we can decompose the Hamiltonian into a center-of-mass term and a relative term, with the relative term being ; for free particles, the relative Hamiltonian is only . Here we have introduced the relative coordinates , , and the reduced mass . We can then write the propagators where .

We use a standard Metropolis procedure, which consists of generating new configurations according to the free particle distribution, and then accepting or rejecting them according to a statistical weight, which takes external potentials and interactions into account. The form (7) is best suited for this procedure, as long as we can efficiently determine the terms under the product symbol. For ease of notation, we now take and define

In order to estimate for the model proposed in Equation (37), we can expand it on the eigenfunctions of the relative Schrödinger equation

For central potentials, which only depend on , like the one considered in this study, the equation splits into an angular part, giving rise to spherical harmonics in d dimensions, and a radial part:

In terms of these wavefunctions, we can expand the relative density matrix into

The coefficients are defined as , for . The functions are the Legendre polynomials of degree l [50]. Finally, the are the bound states of the potential, if any, and their energies; they play no role in the study of repulsive potentials, such as the one we are considering.

The free problem has straightforward solutions that, in the two-dimensional case, yield

In three dimensions, we get

are the Bessel functions of the first kind, and are the spherical Bessel functions of the first kind; the factor appearing in the two-dimensional case is due to the normalization and orthogonality relations. In both cases, the sum can be computed analytically by employing tabulated integrals, leading back to the simple form of (9).

More generally, it is necessary to solve Equations (13) or (14) numerically to find the eigenfunctions. If the interaction is a short-range potential, so that when , for some value of , the eigenfunctions in the region are a generalization of the free case:

and are, respectively, the Bessel and spherical Bessel functions of the second kind. For , it is still necessary to solve the Schrödinger equation. The phase shifts are determined by imposing smoothness conditions on the wavefunction at . In the particular case of a hard-core potential of radius , the requirement is

In order to implement the above formalism efficiently in our simulations, we recast Equation (11) as follows: where

The interacting propagator cannot be calculated analytically; the integrals must be computed numerically and tabulated before running the simulations. While, in principle, we could tabulate the entire propagator as a function of r, , and , the decomposition of (25) has several advantages. First of all, it cleanly separates the contributions from the free propagator, which are relevant at any m at large enough distances, so that we only need to write tables for those values of m that actually present a variation with respect to the free case. Moreover, since the angular variable is explicitly considered in the sum, the integrals only need to be tabulated as a function of r and , reducing computational time and memory usage considerably.

In our two-dimensional simulations, at all temperatures, we have found the contributions from the harmonics to be negligible, so that we only use

Having established the form of the propagator, we can now use it to compute values of thermodynamic observables. Some special care must be devoted to the thermal estimator of the total or kinetic energy, for which the effective potential leads to a contribution of the form with

For a generic of the form with we can introduce it is then possible to show that d being the dimensionality of the system. In particular, this applies to the propagator (28) used in the present work.

As a final technical consideration, we note that the above treatment of the interaction has been carried out in the absence of external potentials, which represent, instead, a crucial component of our physical problem. The Trotter decomposition of (4) allows us to treat different potential and interaction terms independently, or to group them together as needed, as long as we pick a fittingly small value of . For a problem in the harmonic trap, it is most efficient to sample configurations from the harmonic propagator which can be computed analytically [42,51]. We obtained satisfying results by employing this distribution together with the pair-product propagator in the absence of external potentials (11), as displayed in Figure 1, so that the complete form of our propagator is

Figure 1

Comparison between exact density profiles from PIMC simulations (solid lines) and those obtained from the 2D Thomas–Fermi approximation (dashed lines). The shades of purple correspond to different values of the interaction, from small (light) to strong (dark), as indicated in the legend.

3. Application: Trapped Bosons in a Quasicrystal Potential

As motivated in Section 1, we aim to discuss the utilization of PIMC methods, implementing zonal estimators for the quantum properties, in systems displaying a non-periodic patterning. We introduce a model of N identical bosons in continuous two-dimensional space described by the Hamiltonian where m is the particle mass, and are the momentum and position operators of the i-th particle, is the frequency of the two-dimensional harmonic trap confining the bosons. is an external potential defined by with setting the spatial modulation, and the strength of the external potential. The structure of maxima and minima of this potential displays the geometry of the aperiodic, eightfold-symmetrical Ammann–Beenker tiling [52], therefore underlying a quasicrystalline structure. is a contact potential with scattering length a. Figure 2 depicts the configurations of interacting bosons described by Hamiltonian (37) with only the harmonic trap (a) and in the presence of the quasicrystalline potential (b).

Figure 2

Snapshots of particle configurations at , for two values of . Red lines represent the links between successive beads in the equivalent polymer system, as described in Section 2. The shaded background represents the potential, with darker to brighter areas corresponding to lower to higher values. On the left, we show a configuration in the absence of the quasiperiodic potential (), where the only external potential is the harmonic trap. On the right, the presence of the quasiperiodic potential is reflected in the distribution of the polymers, which tend to localize at the minima.

In Figure 3, we show the convergence of the thermodynamic average of the total energy as we increase the number of slices M, and reduce the time step . Varying the strength at , convergence is already established around . Other observables, such as the superfluid fraction, are found to converge even faster.

Figure 3

Trotter limit. The circles correspond to the measured values of the average energy , at decreasing values of the imaginary time step . Different colors represent different values of the potential: (dark blue), (grey blue), (light blue).

We now briefly review the ground-state mean-field approach, which we use to plot density profiles and display the validity of the pair-product approximation. The system is described by a two-dimensional time-independent Gross–Pitaevskii equation [53],

The two-dimensional mean-field parameter g satisfies, up to logarithmic corrections, with the particle density at the center of the two-dimensional trap, and a the scattering length of the repulsive interaction [54,55,56]. Contrary to the three-dimensional case, where the relationship between mean-field constant and scattering length is linear, the logarithm in (41) implies that a must be vanishingly small for small but finite values of g.

We introduce the harmonic oscillator length , and rewrite (41) as where the last approximate equality applies in the limit . We use this approximate form to describe the interaction [37], defining

In Figure 1, we plot the density profiles obtained for interacting bosons in a two-dimensional harmonic trap, sliced across the center. The measured profiles (solid lines) are compared with those predicted by a two-dimensional Thomas–Fermi approximation

In practice, since in our simulations we work at fixed particle number , we first derive as . We find good agreement between the numerical profiles and analytical ones, even for large values of , supporting the validity of the approach used to treat the propagator (36). In the following, we express lengths in units of , and energies in units of ; we employ a weak harmonic trapping with . The connection between the values relative to the harmonic trap and those pertaining to the lattice is the ratio .

4. Density Profiles and Diffraction Patterns

In the present section, we introduce the physical estimators obtained via PIMC, and we discuss the results achieved for trapped bosons in a quasiperiodic potential in two dimensions.

There are several complementary ways to display spatial configuration of the quantum system and its classical-polymer equivalent system. One is to select a system’s configuration at a given simulation step, plotting the position of the beads : the resulting snapshots provide a first graphical estimate of the particle’s probability distribution. In Figure 2 we show snapshots of particle configurations at , for two values of . Red lines represent the links between successive beads in the equivalent polymer system, as described in Section 2. The shaded background represents the potential, with darker to brighter areas corresponding to lower to higher values. On the left, we show a configuration in the absence of the quasiperiodic potential (), where the only external potential is the harmonic trap. On the right, the presence of the quasiperiodic potential is reflected in the distribution of the polymers, which tend to localize at the minima.

Our approach give access to density profiles, which are obtained as averages over simulation steps, as well as over the positions of all different beads associated to each particle. In continuous space, the average is performed by separating the simulation area into bins, and counting the number of beads in each at every simulation step. In Figure 4, we show two-dimensional density patterns for the system at and . Figure 4a is the density profile in the harmonic trap, the same shown in Figure 3. As the value of increases, the quasicrystalline structure appears initially as a modulation (b) and then as localization at the deepest minima (c–d).

Figure 4

Density patterns at and . In each picture, the normalized two-dimensional density takes on values going from 0 (black) to 1 (white). The density patterns are shown at four values of : (a), (b), (c), and (d).

Diffraction patterns can be investigated employing the structure factor that, for example, is observed experimentally in scattering experiments. For a particle density distribution the structure factor reads where is defined as following which is the Fourier transform of the particle distribution (45) [57].

Figure 5 displays some examples of diffraction patterns. As we should expect, the structure factor evolves from a single peak in the fluid phase to a typical quasicrystalline pattern as increases. The three rows in the figure correspond to three different temperatures; we can see that, aside for some smearing of the peaks due to thermal fluctuations, the structure remains essentially unchanged. This is expected since, even above , strong intensities of the lattice lead to localization, and the formation of a classical insulator.

Figure 5

Diffraction patterns at . In each picture, the density in k-space is shown, taking on values from a lower cutoff (black) to the maximum (white). Each column displays a value of : (a,f,k), (b,g,l), (c,h,m), (d,i,n), (e,j,o). Each row displays a value of T: (a–e), (f–j), (k–o).

5. Global Quantum Estimators

A first impression of the importance of quantum effects in the simulation can be obtained by looking at particle permutations [58]. As mentioned above, due to the bosonic nature of the particles, the polymers of the equivalent classical system connect to each other, forming long cycles containing varying number of particles. We call the number of particles in a cycle . At each simulation step, we can construct a histogram, counting how many cycles contain exactly particles for each value of . We can then average the histogram over simulation steps, and normalize it, so that for each value of we obtain the probability of finding a cycle with exactly particles. The histogram of is shown in Figure 6 for different values of . The distribution of in the superfluid phase shows that permutations entail cycles comprising almost all particles in the trap. Intermediate values of lead to the disappearance of permutation cycles comprising of all particles in the system, but still allow for cycles of a few hundred particles. Even larger values lead to a sharp drop in as a function of , with only cycles of the order of ten particles remaining relevant.

Figure 6

Probabilities of permutation cycles as a function of cycle length, at and . Each color corresponds to a different value of , as indicated in the legend. Although histograms are generally more appropriate, we use solid lines for ease of comparison between different cases. All plots are in logarithmic scale.

We further characterize the quantum regime by considering the superfluid fraction. We recall that, in the context of the two-fluid model [57,59,60], the density of a quantum system displaying superfluidity at low temperature can be described by decomposing its density into a sum of two fields: , where the first component describes the normal density, whereas the second is the superfluid one. In this context, the superfluid fraction is defined as the ratio of the superfluid density to :

The two components exhibit contrasting behavior in terms of, for example, flow transport and entropy [61]. The superfluid component displays zero viscosity, and is therefore unresponsive to the application of external velocity fields. In particular, when subjected to an angular velocity field, the superfluid exhibits a reduction of the total moment of inertia compared to a classical fluid in the same conditions. This phenomenon is encoded in a fundamental relation, which links to the system’s moments of inertia:

Here, we indicated I as the moment of inertia related to , representing the moment of inertia, which is the one the same mass of fluid would have if it behaved classically.

By applying linear response theory on Equation (49), one can extract through PIMC. In fact, it is possible to show that the expectation value of the angular momentum in the quantum system is given in terms of the area encircled by tangled paths in the classical system of polymers [3,62]. This way to evaluate Equation (49) results as particularly appropriate for all trapped and finite-size bosonic systems. Since we are dealing with a pure two-dimensional system, we are interested in studying infinitesimal rotations around an axis perpendicular to the plane. Following the formulation given in Ref. [62], in its complete form reads where is the component of the total area enclosed by particle paths on the plane, defined as with , the position of the bead corresponding to the i-th particle on the m-th time slice, the same as in Section 2, whereas the classical moment of inertia in Equation (50) is computed as

Usually, the term is set to zero, or ignored altogether [45]. The most convincing argument is that, from an energetic point-of view, and therefore as far as equilibrium probabilities are concerned, any configuration is equivalent to a symmetric one where the directions of all links between nearest-neighbor beads have been reversed; if both configurations can be accessed, they should be visited with equal probability. Since reversing all links also changes the sign in the area corresponding to the configuration, the immediate consequence is that the average value of the area will be zero, leading to . However, in configurations where bosons are localized into clusters (such as the deepest minima of the quasiperiodic optical lattices, e.g., see Figure 2), this term does not necessarily average to zero, and it must be kept into account. This observation is justified by the fact that the system might spend a long period of time (compared to simulation times) in a region of configuration space where , leading to a manifestation of ergodicity breaking.

In Figure 7, we show plots of as a function of , at different values of temperature and interaction strength. In all cases, deeper quasiperiodic lattices bring about a reduction of the global superfluidity, up to a critical value at which it is completely depleted, and the system transitions to a localized phase.

Figure 7

Measured values of the global superfluid fraction at different values of . Colored markers correspond to different values of the interaction: (teal circles), (black triangles), (green diamonds). Lines are a guide for the eye, and also serve to distinguish different temperatures: (solid lines), (dashes), (dashes and dots), (dots).

6. Zonal Superfluid Estimators

In dealing with inhomogeneous systems, it also worthwhile to extrapolate superfluid features that may be spatially dependent, introducing density fields , , and , which extend the uniform quantities introduced above. We can then combine Equation (48) and Equation (49), and the definition of the classical moment of inertia to write with r representing the distance from the center of the coordinate system. The definition of a local superfluid fraction follows from (48), as ; with this, we can rewrite (53) as meaning that the global superfluid fraction is an average of the local superfluid fraction, weighted by the local moment of inertia. While the introduction of inhomogeneous density fields is natural, their expression in terms of of the classical polymers requires some elaboration; following Kwon et al. [63], who have introduced a physically motivated and consistent definition of , we write where

Equations (53), (55) and (56) allow us to investigate the superfluid behavior locally, usually by sampling (55) on a grid. However, and especially for strongly inhomogeneous systems such as (37), the amount of detail proves to be excessive, and the estimators too noisy. In order to overcome this issue, and to extract some degree of spatial information about the system, we act on a middle level by introducing a zonal estimator.

We divide the system into K regions, labeled by , and introduce a zonal superfluid fraction by specializing (53) to where is the fraction of the total classical moment of inertia corresponding to region k, so that . The zonal superfluid fractions can be recombined, using (54), to give

In the present case, we exploit the circular symmetry provided by the trap to divide the system into three concentric belts. It is important to stress that, for instance, the decomposition (58) may display sectors with a finite superfluid fraction, but still give a negligible contribution to the global , if the associated moment of inertia is small. This is the case for regions close to the trap center, as we show in Figure 8, which displays the superfluid fraction in different regions against temperature.

Figure 8

Temperature behavior of the superfluid fraction in different regions, at . The regions are depicted on the right, where they are superimposed to the density profiles of Figure 4, for (e) and (f). The inner radius is (purple line) and the outer radius is (red line). Colored markers in the plots (a–d) correspond to values of measured in different regions: (blue circles), (purple squares), (red diamonds). The same colors are used to shade the regions in (e,f). We also report the global superfluid fraction, which was already displayed in Figure 7 (black triangles). The four plots each correspond to a different value of : (top left), (top right), (bottom left), (bottom right).

Most familiar is the behavior in the case of , in Figure 8a, where we can see the global superfluidity drop from nearly at low temperature to at the critical temperature . The depletion of superfluidity proceeds at different rates across the trap: the inner regions display a value of , which is still close to 1, even at high values of T, such as . The global superfluidity, however, is dominated by the contribution of bosons in the outer region, due to their larger classical moment of inertia.

As increases, the effects of the quasiperiodic lattice become prominent; the global superfluidity is depleted by localization, as already shown in Figure 7, but, for some values of , the zonal superfluidity remains finite in the central region. In [34], we used this information, coupled with the fluctuations in particle number, to characterize a Bose Glass phase induced by the quasiperiodic potential.

7. Discussion and Conclusions

The present work we detailed PIMC methods to explore both local and global quantum properties of interacting bosons confined in external quasiperiodic potential. Those properties have been probed in a finite temperature regime with specific attention to superfluidity. In detail, we summarized the derivation of the pair-product approximation for particles interacting through hard-core interactions in two or three dimensions, and we presented the details of our implementation. Then, we explored a new zonal estimator, which gives access to local information about the superfluidity in finite regions of trapped systems, and it is therefore well-suited to the study of spatially inhomogeneous potentials. For the example presented in this work, zonal estimators are relevant to the detection of correlated phases, such as the Bose glass phase, which is characterized by rare regions where superfluidity and finite compressibility coexist [34]. Similar zonal estimators can be applied to other quantities, such as regional fluctuations of particle number, associated with density compressibility, and magnetic ordering in multi-component systems or spin systems related to the spin compressibility.

Moreover, one might apply such estimators to the characterization of local properties of self-assembled quasicrystalline phases in free space generated by two-body nonlocal interactions [64]. A prime example of a two-body model potential leading to quasiperiodic patterns is in the paradigmatic hard-soft corona potential, which is largely used to investigate purely classical systems [65]. The same model has also been applied to bosonic systems, where the effects of zero point motion, as well as quantum exchanges, disclose rich phase diagrams including quantum quasicrystal with 12-fold rotational symmetry [17]. Additionally, quantum properties of self-assembled cluster quasicrystals revealed that, in some cases, quantum fluctuations do not jeopardize dodecagonal structures, showing a small but finite local superfluidity [19]. Cluster quasicrystals display peculiar features, not exhibited by simpler quasiperiodic structures. By increasing quantum fluctuations, in fact, a structural transition from quasicrystal to cluster triangular crystal featuring the properties of a supersolid is observed [19,66,67]. We point out that the discussed methodology is also useful to analyze the superfluid character of further peculiar inhomogeneous systems such as, for instance, bosons enclosed within spherical traps or subject to a polyhedral-symmetric substrate potential [68,69,70].