Localisation of weakly interacting bosons in two dimensions: disorder vs lattice geometry effects.

We investigate the effects of disorder and lattice geometry against localisation phenomena in a weakly interacting ultracold bosonic gas confined in a 2D optical lattice. The behaviour of the quantum fluid is studied at the mean-field level performing computational experiments, as a function of disorder strength for lattices of sizes similar to current experiments. Quantification of localisation, away from the Bose glass phase, was obtained directly from the stationary density profiles through a robust statistical analysis of the condensate component, as a function of the disorder amplitude. Our results show a smooth transition, or crossover, to localisation induced by disorder in square and triangular lattices. In contrast, associated to its larger tunneling amplitude, honeycomb lattices show absence of localisation for the same range of disorder strengths and same lattice amplitude, while also exhibiting partial localisation for large disorder amplitudes. We also conclude that the coordination number z have a partial influence on how fast this smooth transition occurs as the system size increases. Signatures of disorder are also found in the ground state energy spectrum, where a continuous distribution emerges instead of a distribution of sharp peaks proper to the system in the absence of disorder.

Introduction

Absence of transport in solids is intimately related with either, interparticle interactions or structural disorder. Here, we focus our attention in localization phenomenon induced by disorder in lattices in two dimensions. The typical model for a real electronic material that predicts localisation for finite disorder strength (Anderson transition)[1], states that at zero temperature () this occurs for dimensions [2]. In other words, a metallic system has always localised states for any disorder strength in . Several refinements to those predictions have been formulated for isotropic media[3-8] having a random distribution of scatterers. However, since Anderson localization is indeed a universal phenomenon, it than can be found also in anisotropic media, and in different situations than those belonging to the context of electronic transport. For instance, propagation of light in solids[9,10] in either anisotropic or orientationally ordered scatterers have been explored. In the condensed matter scenario, it is difficult to isolate the effects of the underlaying lattice structure against those associated to interparticle interactions. In contrast, ultracold neutral gases are perhaps the most versatile and simple setting to explore new possibilities in the study of transport properties, as compared to the analogous real materials with electrically charged particles[11-14]. Neutral ultracold gases confined in an optical lattice allow to have each aspect related to either dynamical or stationary properties, to be externally controlled with high precision[15]. Mott insulating transitions addressed theoretically in boson[16-18], and fermion[17,19,20] cases, have been also experimentally tested in both, Fermi[21] and Bose[22] samples. Even more, Bose glass phases[23,24] and Anderson localization and have been realized in bosonic[25-31] and fermionic systems[32,33]. More recently, the interplay of disorder and interactions has been tested in Bose[34,35] and Fermi[36] gases confined optical lattices. In particular, in that experiment, the difficulty to reach equilibrium configurations caused by disorder and the interplay with strong two body on-site interactions were explored. From the theoretical perspective, there has been an extensive amount of work addressing the study of transport in disordered optical lattices[15,37-40]. Additionally, proposals to investigate fractal structure induced by disorder[41] and measurement of system replicas simultaneously[42] have been put forward. The interplay of disorder with spin-orbit coupling[43], spin degrees of freedom[44,45], supersolidity[46], and Dirac fermions[47] are being explored. New developments to diagnose critical behaviour in the Hubbard model using neural networks have been considered[48], as well as, engineering of bound states via disorder[49]. Some standard schemes used to study the referred phenomenology are approximations to the extended and hybrid Bose-Hubbard model using Gützwiller ansatz[50], Stochastic Mean-Field Theory[15,51,52], Self-energy Functional Theory[53] and the use of Monte-Carlo simulations[54,55]. Numerical exact diagonalization simulations in the full quantum limit are difficult to perform since the size of Hilbert space grows exponentially with the number of lattice sites. Moreover, Renormalization Group schemes are computationally very expensive for systems in dimensions larger than one. Additionally, in the weakly interacting regime Monte-Carlo simulations show strong finite size effects[56]. This makes difficult the characterisation of the states in the simulations and the identification of the corresponding quantum many-body phases of matter. Therefore, other alternatives to address the problem in the lattice with many sites are desirable. A suitable approximation scheme is given by using the Gross-Pitaevskii (GP) equation to effectively describe the ground state of weakly interacting bosonic systems at low temperatures[57,58]. The limitations of this approach are that fluctuations in the number of particles are always poissonian in GP calculations, unable to represent MI states or ensembles of them while depletion is neglected. However, this mean field treatment allows to investigate either stationary or dynamical transport properties in the superfluid (SF) regime (compressible phases) away from insulating states driven by interactions (incompressible phases). The aim of our study is to give an effective picture of the role of geometry in the transport properties of the quantum degenerate gas in the limit of weak interactions in an optical lattice. To this end, we employ an approximate treatment using a mean-field description, GP picture. In this work, we restrict our attention to the characterisation of disorder induced phenomenology and the interplay with the lattice topology via numerical experiments.

In order to identify different signatures of localisation, we describe the weakly interacting bosonic gas in the limit of low temperatures (), where the system is accurately described by a macroscopic wave function describing the condensate component of the system[58,59]. The time evolution of the system is governed by the GP equation under the confinement of the optical lattice. This equation allows us to model the stationary states of the system subject to different lattice geometries, disorder strength , and effective weak on-site interactions. The analysis of our simulations allows us to identify localisation features by varying the disorder strength () and the system size (). The number of minima induced by the optical lattice where the atoms localise having maximum density peaks defines precisely such a system size. In the absence of disorder (), the system is a quantum fluid with translational invariance. As the disorder amplitude is increased, the formation of regions where the density is negligible appear, that is, a continuous density of the condensate is replaced by islands across the 2D lattice. This generates disconnected islands that lead to the suppression of atomic transport across the whole lattice. We arrived to this conclusion after performing a statistical analysis of several quantities, numerically determined, for the condensate wavefunction confined in the 2D disordered lattices. In particular, we studied the average of peak heights as well as the density across the lattice. Fragmentation of the condensate in islands is qualitatively similar to depletion in the sense that a fraction of the condensate extinguishes in certain regions. However, since we are only considering the ground-state of the condensate, that is the wave function, we are unable to consider finite effects[52,60]. As found in the literature, the inhomogeneity introduced by a random potential generates excitations of low energy finite momentum modes[61-65], or causes that a fraction of the superfluid become a normal fluid in the 3D case[66]. It has been found in experiments with optical lattices that depletion affects marginally () the behaviour of the atomic density in low dimensions (1D or 2D) when optical lattice depths are lower than those considered in this work[67]. This supports the use of an effective description given solely in terms of the behaviour of the condensate component in the system as a first approximation. In ultracold gases with weak interactions, this suggests that disorder leads to the fragmentation of the SF in the system. This is the well known scenario that leads to the formation of a Bose-Glass for larger values of disorder strength, which is also a fragmented SF state[57,68]. However, it is important to stress that in our study we considered both, effective interactions and disorder strengths for which the system remains in the superfluid phase, away from the Bose glass and insulating phases.

This work is organised as follows. In the section model, we describe the system and parameters used in our simulations to analyse localisation effects. In section Disorder induced localised states in 2D lattices, we characterise the localisation or density fragmentation, by performing a robust numerical analysis of several observables. Then in section Ground state energy spectra of localised states we study the role of disorder in the 2D lattices in the momentum space. Finally in the last section, we summarise our findings.

Model: Weakly Interacting Ultracold Atoms in 2D Disordered Lattices

To analyse the influence that disorder and the lattice geometry have in producing localised states in ultracold systems, we consider a weakly interacting gas of Bose atoms confined in several two dimensional lattices, subjected to white noise disorder of variable strength. The system is described by the mean field Gross-Pitaevskii equation for the amplitude of the wave function ,where , and being the effective coupling interaction in 2D. [This coupling constant is related to the standard 3D coupling interaction , as , with being the usual coupling interaction written in terms of the wave scattering length , , and being the natural length of the ground state harmonic oscillator of a condensate originally confined in 3D and then squeezed along direction, thus resulting into a disk shape condensate confined in the plane][69-74]. The potential is given by a two dimensional lattice with triangular, square or honeycomb symmetries[75,76],with the lattice constant and , labels identify honeycomb (), square () and triangular lattices (), with being the coordination number (the number of nearest neighbours), , being , , and . The function , representing a non correlated disorder, takes random values varying in the interval at each point being scaled in the same energy units as . Thus, the total potential depth at each point (disordered contribution and optical lattice) is the result of adding/subtracting a random number to the amplitude of the potential at each point . To perform a reliable analysis of the physical quantities and have meaningful predictions, we average over an ensemble of realisations for each value of the disorder amplitude . In Fig. 1 we show a fragment of the energy landscape used in our study, illustrating a particular realisation of disorder in a square geometry. Similar to real experiments, multiple realisations will be considered, as one random realisation fails to represent the typical behaviour of multiple scattering events due to uncorrelated disorder. The way in which the disorder has been simulated, warrants that, although the lattice symmetry is altered, the underlying structure is preserved. It is important to note, that this way of setting disorder through the whole continuous coordinate space where the lattice is defined, resembles analog real systems in condensed matter where each particle “feels” a different energy landscape associated to the disorder in the otherwise perfect lattice. Additionally, this way of implementing disorder is similar to consider a random energy shift at each lattice site, as originally formulated[1]. The disordered potentials experimentally created are the result of either, the light arising from an optical speckle field produced when a laser beam passes through a diffusing plate[26,28,29,77], or the combination of both, a standing wave with a defined lattice geometry and the light arising from an optical speckle field[27,31].

Figure 1

Schematic representation of the two dimensional square disordered lattice potential.

Simulation considerations and experimental parameters

In our simulations, we consider typical parameters of ultracold 87Rb atoms[35] confined in an optical lattice potential in 2D, having depths of and a lattice spacing nm, where is the recoil energy [][78]. At these optical lattice depths the single band Bose Hubbard model[15,18,78] is an accurate description. Other natural scales for the energy are the hopping amplitude between nearest neighbours , and the on-site interaction strength in the optical lattice . [This onsite effective interaction maintains a proportionality relationship with the coupling interaction , given by , being the Wannier functions]. For the lattices here analyzed, the tunneling amplitude can be estimated analytically in the case of the square lattice[78], and numerically in the cases of triangular and honeycomb lattices[79-82]. As it is well known, the appropriate parameter to verify the range of validity of the regime away from strong interaction effects, is the quotient between and . As reported in the literature, for a filling factor of , the SF-MI transition occurs for , and , for square[22], triangular and hexagonal[83] lattices respectively. Therefore, stronger effective on-site interactions than those bounds in each case are needed to be able to represent sub-poissonian on-site number fluctuations, that the GP model is unable to reproduce. In our calculations we use . In this case, the system is very far away from the interaction driven insulating states. Therefore, in our simulations there is no competition with Mott states. Moreover, we neglect quantum depletion effects as we are in the regime[67]. Therefore, as we are deep in the SF regime in 2D, the GP theory is a pertinent description in the presence of disorder.

A quantity that allows us to establish the range where disorder amplitude can be varied before the system enters into the phase in which the density of the condensate becomes discontinuous across the lattice, is the correlation function at the nearest neighbor distance, . This quantity can be defined in analogy to the so called phase coherence that measures the correlation between values of the Bose condensate wave function at points separated one lattice constant . The phase coherence defined as[84] , provides reliable information of the degree of coherence of the condensate inside the optical lattice[84]. In our case, the correlation function between nearest neighbors, is defined aswhere is the coordination number and is a label that identifies the nearest neighbors sites with respect a site located in position. For the particular case of the square lattice we have , while being for the triangular lattice, and for the honeycomb geometry. The correlation gets lost when the density of the condensate is discontinuous, thus preventing the GP approach. Therefore, besides considering values of effective interaction for which the system is well below the MI phase, we shall also consider values of disorder strength that guaranty that the interacting Bose gas remains sufficiently smooth, in order to ensure that our description of the condensate in terms of solutions of a differential equation still makes sense. We should mention here that the phase coherence is a quantity that allows to distinguish when the condensate is not a Bose glass[51]. In Fig. 2 we show the behaviour of as a function of for the analysed lattices. The values of the correlation in this plot are normalised with respect to its value at zero disorder. Error bars in this figure are associated to the ensemble of realisations for each value of (measured in units of recoil energy). As can be appreciated from Fig. 2, for values of below 1, the system remains distributed across the lattice. Thus we can safely consider disorder strengths below one hundred percent of the recoil energy to be in the superfluid regime.

Figure 2

Correlation function between nearest neighbours g1 as a function of δ for the hexagonal (z = 3), square (z = 4) and triangular (z = 6) lattices. Error bars in this figure are associated to the ensemble of 50 realisations for each value of δ.

In typical ultracold experiments besides disorder in the lattice, the atoms move under the influence of a harmonic confinement, so that , with . Although this contribution can also be considered in our mean field analysis, and in fact we performed calculations to investigate the modifications that the inhomogeneity introduces, we note that real analog systems are not under the presence of a harmonic potential, but instead each particle is immersed in a very large (infinite) disordered lattice. Also, near the center of the trap in a typical ultracold atom experiment it is reasonable to neglect harmonic confinement for sufficiently large lattices. Thus, our analysis will focus on lattices without the harmonic confinement for simplicity.

Disorder induced localised states in 2D lattices

Using the effective model given by Eq. (1), we study the emergence of localisation for honeycomb, square and triangular lattices as a function of the disorder amplitude and the system size. In order to describe the ground state of the system, we numerically solve Eq. (1) using imaginary time evolution (), which is equivalent to finding the lowest energy solution or steady-state solution . [The size of the grid was chosen to ensure the norm conservation. The Split-Step Crank-Nicholson method in imaginary time was used to find the stationary state. To ensure the reliability of the stationary state found, we evolved such a state in real time and verified its invariance]. This solution represents the amplitude of the macroscopic wave function of the system at in the lowest energy state. We simulate the GP equation in continuous coordinate space for a given disorder amplitude and analyse the obtained density profiles . As explained previously, the limitations of our treatment are that it does not account for depletion of the condensate, non-poissonian density distributions and it is a zero temperature model. However, it allows to study the combined effects of geometry and disorder, which is the main purpose of this work. The identification of localisation is made from the stationary state by first quantifying the density across the lattice sites, as a function of disorder for a given lattice geometry. Then, to further comprehend the emergence of localisation, we study the behaviour of the stationary localised states as a function of system size. For our calculations, we consider lattices having number of sites used in typical experiments performed in 2D lattices. In order to have meaningful quantities as a function of the disorder amplitude, we consider sets of ~50 realisations of random numbers for each value and a given lattice size . The size , is identified as the number of minima of the potential without disorder. The initial state used in all of our calculations to reach the steady state is a constant distribution [The size of the grid was chosen to ensure the norm conservation. The Split-Step Crank-Nicholson method in imaginary time was used to find the stationary state. To ensure the reliability of the stationary state found, we evolved such a state in real time and verified its invariance].

The information obtained directly from our numerical calculations can be summarised as follows. 1. At zero disorder strength, gaussian peaks centred at lattice sites fill the entire space, being the amplitude of those peaks the same across the lattice, except at the edges where the boundary condition (end of the lattice) give rise to peaks with lower amplitudes, 2. When the magnitude of disorder is non zero, gaussian peaks are not distributed uniformly in the whole lattice, instead of that, several sites in the lattice show a diminished density, thus exhibiting peaks with lower amplitudes. We fixed an arbitrary criterium to count the non negligible peaks contributing to the superfluid density (for a given disorder strength, only peaks having amplitudes larger than 5 percent of the highest amplitude are considered for the statistics), 3. As expected, different realisations of disorder associated to a given value of the disorder amplitude, produce different distributions of peaks across the lattice. What it is important to stress is that different realisations of disorder have, on average, the same distribution of heights and density peaks, 4. As the size of the disorder amplitude increases, the spatially distributed random peaks become sparse, and therefore, in accordance with the condition of constant density, the heights of the peaks become taller, and 5. As the magnitude of disorder is increased, the zero disorder scenario is replaced by a fragmented density across the lattice. For illustration purposes in Fig. 3, we show the stationary density profiles for one of the realisations and different lattice geometries (). These profiles correspond to disorder amplitudes of , and of the value of at each point . The analysis of many realisations of these density distributions is the core of our work.

Figure 3

Fragments of the density profiles for honeycomb, square and triangular lattices with a disorder amplitude δ = 0.2, 0.4, 0.8 with lattice sizes Ω ~ 103. Parameters are: V0 = 12E, U = 0.01E, spacing is in units of the lattice constant a. The bright regions correspond to density maxima while dark regions absence of SF.

In the typical model of localisation induced by disorder without interactions, the envelope of the wave function () gets localised exponentially[1]. In 1D, in an isotropic medium with randomly distributed scatterers where is some arbitrary point in the space, is a localisation length[7]. Within the self consistent theory[8], the localisation length depends on disorder as in 2D and in 1D, where is the scattering mean free path and the wave momentum. In low dimensions (1D and 2D) in the absence of interactions, the localisation length diverges for weak disorder strength, which means that there is no finite value of disorder for which the system can be an ideal conductor. Almost ten years ago, experiments performed in atoms of 87Rb confined in one dimension, in a quasi-periodic lattice[25], and in a weak disordered optical potential[26], showed the typical exponentially localised density profile behaviour. In our problem, instead of studying such a long tail exponential behaviour, we focus our attention on the short distance density variations to characterise the localisation phenomenon. The effects of disorder are captured by performing a robust study of the density profiles. In particular, taking into account the information obtained from our numerical calculations, the quantification of disorder effects in our system will be established in terms of two observables, the peaks heights distribution and the peak fraction , as a function of disorder amplitude . Here we stress that the observables to be studied can be experimentally accessed in actual ultracold atom experiments.

To investigate the dependence of and on we proceed as follows. For fixed disorder strength and lattice size , we performed ~50 realisations of random numbers. From each density profile of the condensate, we analyse the amplitude of the peaks at each lattice site by selecting the non-negligible peaks, using the criteria described above. With this knowledge we determine both, the relative heights and the peak fraction for each realisation. Then, we take the average of over the ensemble. In Fig. 4, we summarise the results obtained for honeycomb, square and triangular lattices. There, we plot the average of the peaks heights (top) and the peak fraction (bottom) as a function of disorder in units of the recoil energy . is referenced with respect to the average heigh value without disorder. Blue, purple and black symbols correspond to honeycomb (), square () and triangular () lattices respectively. As can be appreciated from Fig. 4, square and triangular lattices suggest a smooth transition to localisation for where both, and start to change. In other words, the density collapses into a single density peak randomly located across the lattice. In contrast, localisation is not observed for the honeycomb lattice for the same range of disorder strength. Interestingly, the continuous lines in Fig. 4, fitting the data of each geometry, allow to conclude that there is not sensitivity of the localisation appearance on the lattices having coordinations and . Blue dashed lines in the top and bottom figures of Fig. 4 correspond to the honeycomb lattice, but with a different value of the potential depth with respect to that considered for the continuous lines. While continuous lines correspond to , dashed lines label . Therefore, from this figure one can conclude that the stronger tunneling energy of the honeycomb lattice against the lower for square and triangular lattices, for the same lattice amplitude, plays a determinant role in the occurrence of localisation.

Figure 4

Average of the peaks heights 〈h〉 (top) and the peak fraction p (bottom) as a function of disorder δ in units of E. 〈h〉 is normalised to the largest height amplitude. Blue, purple and black symbols correspond to honeycomb (z = 3), square (z = 4) and triangular (z = 6) lattices respectively. Parameters are: Ω ~ 103, V0 = 12E (solid lines) and V0 = 24E (dashed line for honeycomb), U = 0.01E, with ~50 disorder realisations at each data point, the bars denote the standard deviation for each ensemble of disorder amplitude δ. The honeycomb lattice with V0 = 24E (dashed blue) and the square lattice with V0 = 12E (solid purple) have equivalent tunneling amplitudes in the tight binding limit, see the main text for details. The continuous lines for 〈h〉 are given for triangular and square by: 1 + exp (b0 + b1δ + b2δ2), with parameters (triangular, square): b0,Latt = (−12.53, −7.84), b1,Latt = (31.42, 19.98), and b2,Latt = (−15.58, −8.61). For the honeycomb lattice, 〈h〉 ≈ 1 for with V0 = 12E; for V0 = 24E, b0, = −26.04, b1, = 54.58, and b2, = −25.68. The solid lines for p are given by: 1 − αδ2 for and A exp (β1δ + β2δ2 + β3δ3 + β4δ4) for δ > δ. The fitting parameters are for V0 = 12E (triangular, square, honeycomb): αLatt = (0.29, 0.022, 0.00051), δ ≈ (0.3, 0.3, 0.6), ALatt = (1.023, 0.98, 0.99), β1,Latt = (1.26, 0.0, 0.21), β2,Latt = (−4.17, 3.59, −0.34), β3,Latt = (−1.93, −11.52, 0.005), β4,Latt = 0.0. For the honeycomb lattice with V0 = 24E, α = 0.006, δ ≈ 0.4, A = 0.98, β1, = 0, β2, = 0.69, β3, = 0, β4, = −3.64.

We also investigate finite size effects considering lattices of size . In particular, we study how the density fragmentation of the condensate component occurs as disorder amplitude is increased and the system size is varied. In Fig. 5, we plot the fraction of peaks as a function of for several values of , blue, purple and black points correspond to , and respectively, panels from top to down correspond to honeycomb, square and triangular lattices respectively. Several remarks can be extracted from this figure. First, in our simulations we observe that for disorder amplitudes of the peak fraction never reaches a constant behaviour for triangular and square geometries. In the region we observe that large fluctuations in occur preventing stabilisation. Although this behaviour is similar that observed in typical phase transitions in which critical slowing down and enhanced fluctuations occurs[85], we should stress that the treatment used for our study is unable of describing phase transitions. For , (in Fig. 5 we just illustrate and ) the peak fraction saturates to a small but constant value for these geometries. That is, we see that the fraction of peaks becomes independent of the lattice size for for . This result confirms that lattices with coordinations and suffer a noticeable density fragmentation for disorder strengths above . Although, we observe that the system does not reach the limit where , as in the standard notion of a single localised peak, this result for our lattices in 2D is analogous to the 1D case in which, for non-trapped gases and finite repulsive interactions, the Bose gas populates a finite number of localised single-particle Lifshitz states[68]. In agreement with the results found for lattices of size , honeycomb lattices do not exhibit such fragmentation for the same disorder strengths and the same lattice amplitude. It is important to stress that this density fragmentation observed for the analysed lattices is away of the Bose-Glass phase[57,68,86,87] since, as described in the previous section, the correlation function at the nearest neighbour distance is well above of zero, thus indicating that the system is not an insulator.

Figure 5

Normalized fraction of peaks as a function of the lattice size Ω. The values of the disorder amplitude are δ = 20% (blue), 60% (purple), 80% (black) of the potential depth V0 = 12E. From top to bottom panels correspond to honeycomb, square and triangular lattices respectively. Interestingly, the honeycomb lattice (z = 3) is resilient towards disorder induced localisation, as compared to other geometries with more nearest neighbours (z = 4, 6). The fitted lines different from p ~ 1 are given by: p ~ B exp (−cΩ) + D. For honeycomb: δ = 0.8, (B, D, c) = (0.42, 0.8, 0.0017). For square: δ = 0.6, (0.51, 0.03, 0.00073); δ = 0.8, (0.16, 0.0072, 0.00097). For triangular: δ = 0.8, (0.72, 0.15, 0.0023); δ = 0.8, (0.26, 0.06, 0.0037).

Ground state energy spectra of localised states

Complementing our analysis, we investigate the energy spectrum of the stationary states from our simulations. To proceed, we consider the energy associated to the stationary GP equation as a function of the disorder strength and a fixed value of the mean field interaction . The energy functional is given in terms of the stationary state by,where the disorder is considered in the effective potential for honeycomb (), square () and triangular () lattices (see section Model). The integrand is the ground state energy density in continuous coordinate space. To find the ground state energy density in the reciprocal space, namely the ground state energy spectra, we first rewrite the kinetic term in a quadratic form,Then, it follows from the use Parseval’s theorem[88], the expression in momentum space,with , and unit vectors in 2D. As it turns out from Eq. (4), the ground state energy spectrum associated to a particular realisation of disorder encodes the structure of the ground state energy in momentum space.

Analogously to the previous analysis for characterisation of density profiles, the ground state spectra analysis was done considering a significant number of realisations for a given disorder amplitude. Specifically, the results presented here correspond to the average over ~50 realisations. For all of our numerical calculations, we use and work within the first Brillouin zone (FBZ). For square lattices, this region is bounded by , while for honeycomb and triangular lattices the FBZ is the area of the primitive cell in the reciprocal lattice space generated by vectors and . In Fig. 6, we show the energy spectra of the lattices considered in our study for (left and right columns respectively) of the potential depth , in the first quadrant of the FBZ. From top to down, we plot vs. for honeycomb, square and triangular lattices respectively. The inset of each panel contains a density plot of ground state spectra in the full FBZ, associated to both, (left column) and (right column), we should emphasize that left column corresponds to the averaged ground state spectra. We observe that the system with disorder presents a ground state spectra which is no longer composed of sharp energy peaks, originated from the lattice symmetry. Disorder broadens the peaks in the distribution. In other words, instead of the sharp ground state spectra at zero disorder, we have a dense continuous behaviour when the disorder amplitude is different from zero. Interestingly, even though disorder is spanned across the entire lattice, and that this corresponds to the average over many realisations, we do not observe any cancellation effect due to randomness. Regarding the influence of the disorder strength, we found that as the disorder amplitude is increased the magnitude of the background increases until it becomes comparable to the amplitude of the signal due to the lattice symmetry. At this point, it is when the disorder effectively destroys the lattice symmetry, rendering the geometry effects indistinguishable. Summarising, the influence of disorder on the ground state spectra is manifested in replacing well defined peaks as a function of momenta with a super-imposed dense distribution due to incommensurability of spatial frequencies.

Figure 6

Energy spectrum in reciprocal space for honeycomb (top row), square (middle row) and triangular (bottom row) lattices in arbitrary units. Left and right columns correspond to disorder amplitudes of δ = 0 and δ = 80% of the potential depth V0 = 12E respectively. Each of the energy spectra on the right column is the average over ~50 realizations for a given disorder amplitude. The inset in each panel shows the corresponding density plot of the energy spectrum in the Brillouin zone.

Summary of findings and outlook

We have studied the stationary states of ultracold weakly interacting bosonic atoms confined in two dimensional disordered optical lattices having different lattice geometries. In particular, we investigated the influence of energetic disorder in lattices with honeycomb (graphene like), square and triangular geometries, in producing localised states. The characterisation of localised states was performed by means of a systematic analysis through numerical simulations at mean-field level, using the Gross-Pitaevskii equation. The quantities analysed in our study, as a function of the disorder amplitude, were the average of the heights amplitudes and the fraction of peaks obtained directly from the density profiles, and the ground state spectra determined by Fourier transforming the energy associated to the ground state. For our numerical analysis, we considered lattices of size sites, constant effective on-site interaction strength , and variable disorder amplitude in the interval , being these two quantities such that the Bose gas is well inside in the superfluid region. We also performed a finite size analysis considering the effects of the system size in the range .

We found that localisation in two dimensional lattices occurs as a smooth crossover for square and triangular lattices. On the contrary, honeycomb lattices do not exhibit such localisation for the same disorder amplitudes and the same lattice amplitude . This behaviour can be attributed to the characteristic larger tunneling energy of the honeycomb lattice with respect to that for square and triangular geometries. From the analysis of large lattice sizes, we found that starting from a certain disorder amplitude, the system exhibits partial localisation manifested as disconnected islands of density. We arrived at this conclusion from the analysis of the density since it saturates to a constant value as a function of the disorder strength. The original prediction made by Anderson is not completely restored for our weakly interacting disordered lattices in the sense that arbitrary disorder values do not produce localisation effects. Being this result analogous to that occurring in 1D, in which a weakly interacting Bose gas populates a finite number of localised single-particle Lifshitz states[68]. However, localisation in few sites is found for large disorder strengths. Therefore, although localisation or density fragmentation occurs as a smooth crossover, we can identify regions with marginal localisation for , weak localisation for and strong localisation for . We recognised such strong localisation as the collapse of the condensate density into a single density peak randomly located across the lattice, the limit when . Thus, transport is strongly suppressed for any lattice geometry for sufficiently large disorder strength. Our analysis also allowed to conclude that the coordination number determines how fast the smooth transition occurs as the system size increases, being the triangular lattice the first in which transport in inhibited against square and honeycomb lattices. It is interesting to note that although the honeycomb lattice do not constitute a bonafide Bravais lattice in 2D since it is composed by the superposition of two sub-lattices, our results are in good agreement with the fact that this lattice is an intrinsic good conductor.

From the analysis of the ground state spectra, we conclude that the influence of disorder is to replace a sharply peaked distribution for a continuous one. The disordered system ground state spectra is composed by continuous energy levels originated from disorder and the sharp energy level contribution originated from the underlaying symmetry of the lattice. In essence, the spatial random inhomogeneous distribution translates into broadening of the characteristic peaks of perfect lattices in the ground state energy distribution in momentum space.

Summarising, in our work we have been able to quantify the crossover to localisation as a function of the disorder amplitude for three different geometries. A general statement arising from our investigation is that disordered honeycomb lattices are robust against disorder as compared to triangular and square geometries for the same parameter ranges. As mentioned above, associated to the larger tunneling energy of the honeycomb lattices against their respective in square and triangular lattices, the honeycomb structure would have better conductivity in the presence of disorder in 2D.

The present work provides some insight to understand localisation phenomena induced by disorder in two dimensional weakly interacting systems for lattices of different geometries having coordination and 6. It is interesting to note although Gross-Pitaevskii equation has an intrinsic mean-field nature, the quantities analysed to investigate localisation phenomena, namely the density distribution across the lattice and the ground state spectra, indeed exhibit signatures associated to the lattice geometry. Additionally, the information displayed in the ground state spectra in the reciprocal space can be compared with direct experimental measurements in ultracold systems via time of flight measurements[89]. Further investigations of transport phenomena can be addressed with the approach followed in the present work. For instance, it is possible to consider in a square lattice, different disorder amplitudes in and , that is, setting . In view of the results obtained here, this could trigger the formation of stripes where transport could be manipulated via disorder induced localisation. Moreover, dynamical effects in adiabatic and non-adiabatic transfer protocols for the engineering of quantum states could be analysed[90,91]. Beyond ultracold atoms, systems of ions[92], superconducting devices (Circuit-QED)[93], polariton systems[94,95] and the addition of high-Q cavities[96] offer opportunities to test our findings analogously via quantum simulation[97,98]. Moreover, the interplay of these effects in lattices with dissipation offer another venue of exploration[99,100] or measurement induced dynamics[101]. The findings of our work can aid in the development and characterise of transport properties in two dimensional systems for the development of new analog quantum technologies[102,103].