Metastable Patterns in one- and two-component dipolar Bose-Einstein Condensates

In this paper we study metastable states in single- and two-component dipolar Bose-Einstein condensates. We show that this system supports a rich spectrum of symmetries that are remarkably stable despite not being ground states. In a parameter region where striped phases are ground states, we find such metastable states that are energetically favourable compared to triangular and honeycomb lattices. Among these metastable states we report a peculiar ring-lattice state, which is led by the competition between triangular and honeycomb symmetries and rarely seen in other systems. In the case of dipolar mixtures we show that via tuning the miscibility these states can be stabilized in a broader domain by utilising inter-species interactions.


I. INTRODUCTION
Ultracold quantum gases with long-range interactions display remarkably intriguing behaviour and give access to probing fundamental quantum behaviour [1,2].Dipolar Bose-Einstein condensates (BECs) permit a controlled access to such effects [3,4].These include the recent experimental advances in the observation of quantum droplets [5][6][7], supersolids [8][9][10][11][12][13][14][15][16] and their excitation spectra [17][18][19][20].The existence of these states of matter can be viewed as a macroscopic signature for quantum fluctuations in that they suppress dipolar collapse that would otherwise occur [21].With that dipolar quantum gases emerged as an ideal platform to observe fundamentally interesting and surprising physical effects.One of these effects is the emergence of a point in phase-space where the superfluid-supersolid phase-transition becomes second-order [22], which means around that point the physics becomes practically linear in the modulation amplitude and acquires a glass-like nature and rich variety of patterns [23][24][25].
Generally, considering more than one species of atoms in ultracold quantum gases [26,27] adds complexity and can promote a range of intricate phenomena, including collapse suppression even for short-ranged interactions [28,29] by quantum-fluctuations in mixtures as well as tunable miscibility [30,31] between the components.
The recent observation of dipolar mixtures [32][33][34] paves the way for a range of new perspectives and qualitatively new behaviour as it permits the combination of the intriguing behaviour of two-component physics with long-range interactions [35][36][37][38].
Metastability is ubiquitous in nature ranging from physics [39,40], chemistry [41] to material [42] sci-ence.Metastabilty can lead to a variety of rich phenomena [43][44][45][46][47][48] and permits gaining further insights into the overall physics of the systems.The transition between metastable states has attracted attention as well [49].Previous work have found that long-range interaction is a crucial ingredient for the appearance of metastability [44-48, 50, 51].For example, it has been reported that dipolar atoms loaded in optical lattices can host a number of metastable states in Mott-insulator regimes [45,46].However, it is still unclear whether a continuous diploar gas can support metastable phases that feature a significantly different symmetry to the groundstate as well.
Here, we explore such metastable states that feature multiple length scales in single-as well as two-component dipolar BECs.Close to the second-order point deformations of small density modulations barely lead to a change in energy due to the shallowness of the energylandscape.Therefore, one can imagine that states with different symmetry can be created via linear superposition, however, such superpositions can be expected to be unstable.However, the dynamics can become so slow that the density almost appears frozen.Further away from the second-order point, where the periodic density modulations become larger and interactions lead to a more pronounced energy landscape, the possibility of sufficiently deep local energy minima appears more reasonable.
To address that idea systematically, we start off with considering single-component systems and investigate the emergence of superlattices by superposing two different patterns in regions close to lines where they are energetically degenerate.To avoid finite-size effects we consider the thermodynamic limit.By thermodynamic limit we refer to the situation where the plane perpendicular to the dipolar polarization direction z, which corresponds to both trapping and polarisation axis, is infinitely extended.In this plane, the average two-dimensional (2D) density ρ 2D is fixed, and a 2D symmetry breaking occurs.
After studying the single-component system we consider dipolar mixtures.Here, we have the additional degrees of freedom due to the interaction between the two components.By changing the miscibility we can tune between a situation where both components reach their maximum density at the center of the trap and triplelayered density distributions, where one component is "sandwiched" between two layers of the other component.The interaction between the layers might render such metastable states unstable or possibly stabilize them.

A. Single-component BECs
A single-component dipolar Bose-Einstein condensate at zero-temperature composed of N dipolar Bosonic atoms of mass m including quantum fluctuations can be described via Here, ρ(r) ≡ |ψ(r)| 2 represents the condensate density, we assume trapping along the polarization direction only, in this case the z-direction, and ω z is the respective frequency of the harmonic trap.Furthermore, g = as 3a dd is the fraction of the s-wave scattering length a s to the dipolar length a dd , the latter of which characterizing the dipole-dipole interaction strength.V (r) = 1 4πr 3 (1 − 3z 2 /r 2 ) is the usual dipole-dipole interaction.The Lee-Huang-Yang (LHY) correction due to quantum fluctuations µ LHY is given by [52][53][54][55][56].

B. Two-component BECs
To generalize the single-component description to model two-component dipolar condensates we assume the local density approximation and employ the model recently introduced in [35,36]: where α, β = 1, 2, ρ α (r) ≡ |ψ α (r)| 2 , ω z is the frequency of the harmonic trap along the polarization direction, is the usual dipoledipole interaction with a αβ s being the s-wave scattering length and a αα dd the typical dipolar length characterizing the dipole-dipole interaction strength.And the LHY correction µ (α) LHY is given by with For simplicity, we assume that the atomic masses of the two components are equal (i.e., m 1 = m 2 = m), which is justified for the typical Dy-Dy [35] as well as reasonable for Dy-Er [36] mixtures, and the above equations have been nondimensionalized through scaling spatial coordinates and time by l = 12πa 11 dd and ml 2 /ℏ, respectively.Hereafter, our discussion will focus on the Dy-Er mixture, i.e., a 11 dd = 132a 0 and a 22 dd = 65.5a0 with a 0 being the Bohr radius.

III. RESULTS
In this section we present that both in single-as well as in two-component dipolar BECs we can find metastable states that can even feature two length scales, despite the fact that there is only one roton minimum in the dispersion relation.For that matter, in the first subsection III A we illustrate in a single component BEC that new states with multiple length scales can be thought of as a certain superposition of other metastable states.Then, in III B we show that similar arguments also apply for two-component systems.

A. Single-Component BECs
Let us start with the single component system and present the groundstate phase-diagram that has already been presented for a finite density, i.e. trapped in all three spatial directions, in [23,24] and in the thermodynamic limit including the stripe phase in [25,57].
The groundstate phase-diagram is shown in Fig. 1.We find that all phases emerge from the point where the superfluid-supersolid phase-transition becomes secondorder and coexistence terminates.The region where the stripe-phase is the ground state can be separated into two areas, depending on whether the honeycomb (downhexagons) or the triangular (hexagonal) states are energetically favorable with respect to each other.This transition is indicated by the white dashed line in Fig. 1.Whereas in this area both the honeycomb as well as the triangular state are metastable, they feature remarkable stability and robustness.Therefore, it appears reasonable to ask whether this system robustly supports also more complex metastable states, such as states with multiple length-scales, two examples of which we present in the following.
A natural approach for finding metastable states that feature more involved geometries is to inspect regions in the phase-diagram Fig. 1 close to where the meta-stable states (i.e.honeycomb/triangular) feature equal energy, that is close to the earlier mentioned dashed line in Fig. 1.That is due to the fact that at these points the system does not favour either of them and, therefore, one might be tempted to expect that in this region they can be admixed in some way.This will be used to explain the emergence of surprising metastable phases in the next subsection.

Rings and ring-droplets
Let us come back to the numerical results shown in Fig. 1 and discuss the yellow region of Fig. 1, where isolated ring-like densities emerge [cf.Fig. 2(g)].It appears that the yellow region shrinks upon approaching the second-order point, but does not converge to the secondorder point.This can be expected, as only groundstates can converge to the second-order point, whereas supporting meta-stable states requires sufficiently large amplitude modulations, as the latter permit the formation of a sufficiently deep local energy minima.Therefore, the nu-merical results matches the qualitative expectation: Initially the "center-of-mass" of the ring-region appears to nicely follow the dashed transition-line upon decreasing density, yet at a certain point the domain has to stay sufficiently far away from the second-order point and deviates from that trend.
Whereas the ring state is metastable with respect to the stripe phase, it features a lower energy than both the triangular as well as the honeycomb density distribution for the solid yellow lines.Upon deviating from the dashed white line, the ring state becomes energetically unfavourable compared to both triangular and honeycomb states, yet continues to exist.To distinguish that case we draw a dashed yellow line rather than a solid yellow line in Fig. 1 closest to the second-order point.
After this mostly qualitative discussion, let us now inspect the energies of the states involved in more detail.The energy per particle E is given by the functional where N is the total particle number and the wave function ψ(r) has been normalized to 1.
The situation is plotted in Fig. 2 for a fixed density ρ 2D = 312.5.Fig. 2(a) shows the comparison of the energies of different states supported by the system in the region where the striped density distribution is the ground state as function of the a s /a dd .To add further visual clarity, Fig. 2(b) shows the energy difference between the respective states to the straight dashed line as depicted in Fig. 2(a).
We find that for small values of a s /a dd the ground state is given by a triangular distribution of density droplets Fig. 2(c).Upon increasing a s /a dd sufficiently the groundstate corresponds to a stripe phase shown in Fig. 2(d) and finally, for large values of a s /a dd is a honeycomb as depicted in Fig. 2(e).
Let us now discuss the metastable states of Fig. 2(b) shown in Fig. 2(f,g) and start with the domain where the triangular or hexagonal lattice is the ground state.One might be tempted to expect that the stripe state ought to be the first metastable state, as it becomes the groundstate for larger a s /a dd .However, for small a s /a dd the first metastable state we find corresponds to a droplet lattice that features two length scales [see Fig. 2(f)].Upon further increasing a s /a dd we find that the stripe phase becomes the first metastable state before becoming the groundstate.In the region where the stripe phase is the groundstate we can again distuinguish two regions, one where the triangular lattice is the first metastable state (small a s /a dd ) and one that emerges upon increasing a s /a dd .This state is curious, as it features a density distribution that is ring-like [see Fig.Let us now return to the discussion of Fig. 1, but from the perspective of the length scales involved rather than the earlier presented energetic arguments.The deviation of the "center-of-mass" of the ring-state region from the white dashed line can be understood as follows.At the second-order point co-existence of all phases with a single, fixed wave-vector converge.This wave vector k can be found by the Bogoliubov excitation of the unmodulated state.In order to obtain a density distribution that features multiple length scales, the dynamics needs to be "sufficiently nonlinear" in the modulation amplitude, i.e., sufficiently far away from the second-order point at which there can only be one length scale λ = 2π/k.It is the nonlinearity that gives rise to a second length scale.This is again consistent with the fact that these ring states typically feature a large contrast and we were not able to find ring states with small modulation amplitude.
Rings on a triangular lattice are dramatically different from the earlier mentioned more common patterns, as they can be thought of as a mix of both triangular and honeycomb lattice.To understand how they emerge and to furthermore identify the two earlier mentioned lengthscales, let us superpose two states such that metastable ring states close to the white dashed line in Fig. 1 in the following fashion: Here, we assume 0 < A ≪ 1 as an ansatz for a weakly modulated condensate with triangular symmetry ρ T or honeycomb symmetry ρ H , respectively.A denotes the small amplitude of the density modulation.The three wave vectors form an equilateral triangle in the transverse plane with k 1 + k 2 + k 3 = 0 and |k j | = k.
To see how we can obtain a ring state with these states, let us superpose a triangular droplet lattice ρ T with a honeycomb lattice ρ H in the following manner: Here, R(ϕ) denotes a rotation perpendicular to the polarization direction by an angle of ϕ and T represents a translational shift operation.The idea of this ansatz is, basically, that the ring state looks similar to a honeycomb lattice with certain connections [see the yellow dashed circles in Fig. 3] being removed.To remove these connections as well as the background, we subtracted a triangular lattice.To further clarify the situation, the superposition process is shown in Fig. 3.This ansatz appears to capture the essence of the ring state, and we can directly read off that there are two sets of wavevectors k: a (cos(2π/n), sin(2π/n)) and k ′ n = 2π a (cos(2π/n + π/6), sin(2π/n + π/6)) with n = 1 . . .6 [also see Fig. 4(c)].Here, a denotes the lattice constant.
Surprisingly, even upon employing variational analysis, we find that such ring states can indeed be energetically preferable compared to both honeycomb as well as triangular lattice state, however, the stripe phase remains the ground state.To make sure that this result is not an artificial feature due to the variational ansatz, we show numerical results to confirm both existence as well as stability of such ring states.
To validate whether this ansatz is consistent with the states we found numerically earlier in Fig. 2, we also show the Fourier transform of the ring state as well as the triangular and honeycomb states in Fig. 4. Evidently, the symmetry is correct and the Fourier transform shows what has been expected from this linear expansion [58].Interestingly, this state resembles the so-called fairy circles that can be found in waterscarce areas [59][60][61][62], as it shares the property of multiple inherent lengthscales (ring-size and ring-to-ring distance).Moreover, we would like to point out this ring state solely results from the strong nonlinear effect of dipolar condensates.This is in sharp contrast to similar phenomena observed in spin-orbit coupled BECs as well as unbalanced binary atomic mixtures, where the ring-like lattices are led by either gauge fields or external trapping with a ring geometry [63][64][65][66].

B. Two-component BECs
It is not priori clear whether our findings in the single-component case can be simply transferred to twocomponent systems.The biggest difference is the possibility to form miscible and immiscible or layered structures, and due to the cross-component interaction these layers interact.This interaction can be thought of as either being inhibiting or catalysing the formation of metastable states with multiple length scales.
As a starting point to explore this more complicated system it is reasonable to first analytically write down an expression that captures the immiscibility in a two component system neglecting the modulations.In other words, we first consider unmodulated states and only focus on the new feature of immiscibility in this simplified system.After that we aim at transcribing the findings from single-component to the two-component case.

Miscibility of two-component unmodulated states
In this section we seek to find analytical approximations for the unmodulated density distributions including that captures their miscibility.As mentioned before, twocomponent systems can have different degrees of miscibility depending on the parameter regime [30,31,35,36,67].
Here, we consider an unmodulated state for both components and approximate the density distribution (i.e., ρ 1 and ρ 2 ) of that state with a Thomas-Fermi profile in the trapping direction z.After a bit of algebra, that allows to find with For η = 1 represents the point where the transition from the miscible to the immiscible regime occurs.In the miscible regime (η < 1), both ρ 1 and ρ 2 reach maximum value at z = 0 as can be seen from Fig. 5(a).In contrast, for η > 1 the density of ρ 1 becomes a parabola close to z = 0 [see Eq. ( 9) and Fig. 5(b,c)].Hence, ρ 1 no longer has its maximum at z = 0, but it acquires two maxima at the point where ρ 2 (z) becomes zero, while ρ 2 retains the normal Thomas-Fermi profile.In contrast, both ρ 1 and ρ 2 deviates from the normal Thomas-Fermi distribution when η > µ 1 /µ 2 , where the strong contact inter-species repulsion eventually fully separates the two components as shown in Fig. 5(d,e)].In that case, ρ 1 vanishes in a finite region around z = 0.In other words the binary condensates enter the immiscible regime.
In the strong immiscible regime the density profiles of the two species can be expressed as below where the chemical potentials µ 1,2 take the same form as in the former case.The critical points σ z1,z2,z3 can be determined by the constraint ρ 1,2 (z)dz = ρ 2D .For the case of intermediate cross interaction [cf.Fig. 6(d)], there remains overlap of the two species in the region σ z3 ≤ |z| ≤ σ z2 .However, if the inter-species contact interaction exceeds the critical strength (a 12 s ) c = (a 11 s + 2a 11 dd )(a 22 s + 2a 22 dd ) − 2 a 11 dd a 22 dd , the overlap between the two components disappears completely [see Fig. 6(e)] and it enters the completely immiscible regime.

Excitation spectrum of two-component unmodulated states
For convenience we report here an (approximated) excitation spectrum of two-component unmodulated states as their expressions are useful for the next subsection.The (approximated) excitation spectrum can be easily obtained from linearisation of the governing equation of motion Eq. ( 3).For that matter, we use the expression obtained before Eq. ( 9) for ρ 1 (z), ρ 2 (z) and integrate out the z-direction.With that we obtain the following expression that only depends on k ⊥ as follows: Here ρ 2D α = ρ α (z)dz (α = 1, 2) denotes the average 2D density of the condensate in the plane perpendicular to the polarization direction and corresponds to the Fourier transform of ρ α and g αβ = contact interaction.To obtain the above simple analytical excitation spectrum, we have neglected the LHY correction terms.This is not at all justified and serves for a qualitative discussion only and is certainly quantitatively wrong.The blue solid line in figure Fig. 6 shows the position where he roton minimum of the unmodulated state's excitation spectrum touches zero as function of the scattering length a 12 s .Remarkably, this line displays distinct behaviors in the weak and strong cross contact interaction regime.In the case of small scattering length a 12 s , the roton instability can be promoted either by the intercomponent interaction a 12 s or by the intra-component interaction a ii s under a general constraint.In contrast, in the strong cross interaction regime, this critical line no longer depends on a 12  s and is solely determined by the intra-component interaction a ii dd .This can be understood as follows.As can be seen from Fig. 5(e), there is no overlap between the distributions of the two components at large a 12 s .Therefore, g 12 becomes zero and the excitation spectrum ω depends on the intracomponent contact interaction only [see Eq. ( 12)].In contrast, in case of small a 12 s , the two components are miscible and their overlap is determined by their cross interaction.This is why the excitation spectrum exhibits the two different scalings.s /a 11 dd = a 22 s /a 22 dd .Again, the ring state remains metastable with respect to the stripe phase, however, in the region indicated by the dashed yellow line it is energetically preferred compared to the triangular and honeycomb states.The overall trend of all lines can be estimated via the approximated excitation spectrum Eq. ( 12), which amounts to the blue solid line.The two horizontal lines correspond to the values of a ii s /a ii dd (i = 1, 2) for which the miscibility ρ 2 (z=0)−ρ 1 (z=0) ρ 2 (z=0)+ρ 1 (z=0) is shown in Fig. 5.We find that the ring state persists for a large range of miscibility.Furthermore, for the values at hand, it appears that the miscibility does not affect the existence regions of the states, as despite of a dramatic change in the former, the overall trend guessed from the (approximated) excitation spectrum remains.The domain of triangular states is separated into two regions by the white solid line.Left to this line the groundstate is the usual triangular lattice, whereas right to this line the groundstate features the triangular superlattice that will be discussed further in subsection III B 5.
The qualitative behavior is confirmed by numerical simulations [see the black solid line in Fig. 6].The significant difference between the full numerical results (black line) and the approximated excitation-spectrum (blue line) is clearly visible, and it is clear that for quantita-tively correct results the LHY correction has to be taken into account.Since the LHY correction acts like contact repulsive interactions which tends to stabilize the unmodulated state, the roton-instability critical line would be shifted towards smaller a ij s by quantum fluctuations as the numerical results show.

Rings and ring-droplets in two-component systems
We already established that we can find ring states in the single component system, however, it is unclear whether we can find their analogue in two-component systems.That is due to the fact that the different layers interact and thereby might prevent the formation of ringstates.
To explore that, consider Fig. 7(a), which displays energy differences to the energy of the stripe phase.The supplementary figures (b,c,d,e) display the densities of the different states involved.We find again that close to where the triangular and honeycomb state become energetically degenerate and the striped state is the ground state, a metastable ring state emerges.Close to the aforementioned energy degeneration, this ring state features a lower energy than both the triangular state as well as the honeycomb state, akin to what we already found in the single-component case.= 0.9.The ground state for most of the considered region remains a stripe phase.In analogy to the the single component case shown in Fig. 1, we find that in a two-component BEC ring states are possible as well.As in the single-component BEC, these ring states occur in a region close to where honeycomb and triangular lattices become energetically comparable.Now that we established that ring states are also possible in dipolar BEC mixtures, let us address whether their existence is promoted or suppressed when tuning the miscibility.As we have seen before in Eq. (10), one can alter the miscibility by tuning the value of a 12 s .The result is provided in Fig. 6, where we show the phase diagram for a fixed density and varying cross-contact interaction and contact interaction in the upper panel (a).Here, we focus on the case of balanced intra-component interactions, i.e., a 11 s /a 11 dd = a 22 s /a 22 dd .We see that for this set of parameters the regions of the different phases are practically unaffected despite tuning the miscibility dramatically from 0.3 to 1 [see Fig. 6].In other words, the ring state, rather than ceasing to exist, even retains its domain of existence.This shows that the ring states are robust, as the cross-interaction can be considered as a perturbation to the single-component physics that deforms the energy landscape significantly.
In fact, the opposite is true for the displayed case: The cross-interaction with the other component actually stabilizes the ring state as compare to the single-component case beyond the otherwise critical line of a s /a dd shown in Fig. 1 of around a s /a dd ≈ 0.79, where in the single component case all modulated states cease to exist in favor of the unmodulated state.Therefore, in a way one might say that these somewhat peculiar states actually emerge in a broad parameter region in both single-as well as two-component dipolar BECs.Moreover, the cross interactions in two-component systems can stabilize the existence of these exotic states.8.The superfluid fraction of different stable states vary with the interspecies contact interaction while the intraspecies interaction is fixed at a 11 s /a 11 dd = a 22 s /a 22 dd = 0.9, where the full (dashed) lines indicate the superfluid fraction of the first (second) component.Here the red, black, dark yellow, and blue lines represent the corresponding superfluid fraction f s 1 (f s 2 ) of the triangular, stripe, ring, and honeycomb states, respectively.The regime of the stripe ground state is denoted by the light gray color, while the dark gray shadow illustrates that the ring state is energetically favorable than the triangular and honeycomb states in this region.
Thus far we did not consider the superfluid properties the different patterns we We will do that now using Legget's bound [68,69].
where α = 1, 2, L x is the size of the numerical box along x direction, and we take the minimum with respect to all possible directions defined by the angle θ with x = x cos θ − y sin θ and ȳ = x sin θ + y cos θ [22].This is a useful quantity to estimate the superfluid fraction of the condensate [70,71] and has recently been refined [72].Fig. 8 presents the superfluid fraction of the triangular (red lines), stripe (black lines), ring (dark yellow lines) as well as honeycomb (blue lines) states in the vicinity of the critical region from modulated states to a flat state.Here, the interspecies interaction is fixed at a 11 s /a 11 dd = a 22 s /a 22 dd = 0.9 (i.e., along the black dashed line in Fig. 6).We note that the first component (solid lines) possesses a large superfluid fraction compared to the superfluid fraction of the second component (dashed lines) which is practically negligible for small values of a 12 s .However, for sufficiently large values of a 12 s , such that the honeycomb can form in the second component (blue dashed line), also the second component is able to increase its superfluid fraction significantly due to the overlapping wavefunction, i.e. atoms can flow along the bridges of the honeycomb [22].
To gain further insight into why the two components behave so distinctively different, let's consider the singlecomponent case.For the second component this individual species has a much lower dipolar length (i.e. a 22 dd = 65.5a0 vs. a 11 dd = 132a 0 ), it would require a much larger average 2D density to reach the second-order point where superfluidity is large and from where new phases with higher superfluid fraction emerge [22].Thus, for the case of balanced dipolar mixtures considered here, the second component remains almost insulating.Nevertheless, as stated before the situation changes when second component enters the honeycomb regime.

Deformed triangular states at strong cross interaction
As discussed in Sec.III B 1, the two-component dipolar mixtures feature a transition from a miscible to an immiscible distribution that can be controlled by the cross contact interaction a 12 s .Thus far, the discussion of miscibility has been limited to polarisation direction as we considered states that are unmodulated in the transverse plane (perpendicular to the polarisation direction).
What is missing is whether there is a parameter region of modulated states where we can see clear signatures of immiscibility in the transverse plane.
As a proof-of-principle that such regions exist and to illustrate the effect of "transverse immiscibility", we go into the region of deeply crystallized triangular states that appears for large values of a 12 s (cf.Fig. 6, right to the white line).Fig. 9 shows an example of deformed triangular states at strong cross contact interaction a 12 s / a 11 dd a 22 dd = 0.86.From the 3D slice plot in the panels (a) and (b), one can discern that the first component is not fully separated from the second component in z (not completely immiscible).Instead, new "droplets" turn up at the center of three neighbouring regular droplets as shown in Fig. 9(a) and (c).Furthermore, these new "droplets" spread across z = 0 plane and bridge between the upper and lower layers of the first component.The reason that the triangular state is composed of well separated droplets in the deep modulation regime.Upon increasing the cross contact interaction, the first component is pushed towards large z in order to minimize the total energy by reducing the overlap between the two components.If a 12 s is too large, part of the atoms prefer to stay around z = 0 with low trapping potential energy, and eventually give rise to these new links between the two dominant layers.In this situation the second component can no longer maintain its regular triangular distribution and resorts to a deformed triangular profile as displayed in Fig. 9(d).As mentioned, the phase-boundary between regular and deformed triangular states is denoted by the white line in Fig. 6.In contrast, the fully immiscible honeycomb or flat states do not have links in between the outside layers, since the second component possesses a significant superfluid background which prevents the formation of such layer links due to the strong interspecies contact repulsion.

IV. CONCLUSIONS
In this paper we established that dipolar BECs support unusual metastable robust states featuring multiple length scales.These states have the shape of a ringlike density distribution whose azimuthal density modulation can be modulated via tuning the scattering length.They appear in a domain where the stripe phase is the ground state and in a region around the line where the (metastable) triangular and honeycomb lattice become energetically degenerate.
Moreover, these states can be stabilized in a much broader regime in the binary dipolar BECs where only unmodulated flat state exist in the single-component counterparts.In sharp contrast to other ring-like phenomena induced by gauge fields and ring-shape confinement in spin-orbit coupled or unbalanced quantum gas mixtures, the ring-lattice state reported here is solely induced by the strong nonlinear effects.
Although such ring states do not emerge as ground states, they underpin the variety of stable self-organized structures in long-range interacting systems, as already featured in, e.g., [45,46,[73][74][75].Therefore, dipolar BECs represent a promising platform to explore metastablestate quantum phase transition as well [39,76].
In addition to the various metastable states emerging close to where the triangular and honeycomb states become energetically degenerate, we also discovered a deformed triangular superlattice groundstate in the deeply modulated and immiscible regime.Usually, immiscibility is discussed as a phenomenon occurring in the polarization direction (e.g., [35,36]) rather than the plane transverse to the polarization.In this case the emergence of the triangular superlattices with periodic density-bridges displayed in Fig. 9 is a clear new feature due to "transverse immiscibility" of two-component BECs that is different to what has been discussed and has no analogue with a single component dipolar BEC.
As an outlook, we think that further understanding of the topology of this phase-diagram can be found by studying the bifurcation diagram [77].Its dependency on temperature seems to represent an interesting endeavour as well [78][79][80].Furthermore, in this work we restricted our consideration to a small subset of parameters, and extending that to e.g.unequal intra-component interactions and unequal masses to explore the rich spectrum of metastable states remains to be done.
2(g)] and is rarely seen in other pattern-forming systems.The state shown in Fig.2(f) resembles the state Fig.2(g), as it carries the same underlying symmetry apart from an additional azimuthal modulation along the ring.

FIG. 2
FIG. 2. (a) depicts the energy landscape of the ground as well as metastable states at the density ρ2D = 312.5.The energy differences of each state with respect to the dashed line in (a) is plotted in (b) to represent the energies of different states more clearly.The density profiles at z = 0 of the triangular, stripe, honeycomb, ring-droplet, and ring states are shown in the subplots (c-g), respectively.Here the stripe state presented in (d), the ring state shown in (g) and the honeycomb state exhibited in (e) correspond to the phases at as/a dd = 0.76, 0.75 and 0.78, respectively, while the remaining states are at as/a dd = 0.74.

FIG. 3 .
FIG. 3. (a) shows a honeycomb lattice ρ H that has been shifted to locate its minimum on top of the maximum of the triangular state and rotated by π/6, (b) a triangular lattice ρ T , (c) the superposition of the two lattices by subtracting (b) from (a).

FIG. 4 .
FIG. 4. Fourier transforms of the (a) triangular, (b) honeycomb, and (c) ring states.This shows that ring state emerges as a mixture between a honeycomb lattice with a triangular lattice.

FIG. 6 .
FIG. 6.The phase diagram of two-component dipolar condensates with fixed density of ρ 2D 1 = ρ 2D 2 = 625 and balanced intra-component interactions a 11s /a 11 dd = a 22 s /a 22 dd .Again, the ring state remains metastable with respect to the stripe phase, however, in the region indicated by the dashed yellow line it is energetically preferred compared to the triangular and honeycomb states.The overall trend of all lines can be estimated via the approximated excitation spectrum Eq. (12), which amounts to the blue solid line.The two horizontal lines correspond to the values of a ii s /a ii dd (i = 1, 2) for which the miscibility ρ 2 (z=0)−ρ 1 (z=0) ρ 2 (z=0)+ρ 1 (z=0) is shown in Fig.5.We find that the ring state persists for a large range of miscibility.Furthermore, for the values at hand, it appears that the miscibility does not affect the existence regions of the states, as despite of a dramatic change in the former, the overall trend guessed from the (approximated) excitation spectrum remains.The domain of triangular states is separated into two regions by the white solid line.Left to this line the groundstate is the usual triangular lattice, whereas right to this line the groundstate features the triangular superlattice that will be discussed further in subsection III B 5.

FIG. 7 .
FIG. 7. The energy for four different stable states is shown in (a) for a 11 s a 11 dd FIG.8.The superfluid fraction of different stable states vary with the interspecies contact interaction while the intraspecies interaction is fixed at a 11 s /a 11 dd = a 22 s /a 22 dd = 0.9, where the full (dashed) lines indicate the superfluid fraction of the first (second) component.Here the red, black, dark yellow, and blue lines represent the corresponding superfluid fraction f s 1 (f s 2 ) of the triangular, stripe, ring, and honeycomb states, respectively.The regime of the stripe ground state is denoted by the light gray color, while the dark gray shadow illustrates that the ring state is energetically favorable than the triangular and honeycomb states in this region.

FIG. 9 .
FIG.9.The density profiles of the deformed triangular states at strong inter-species contact interaction, i.e., the right side of the white line in Fig.6(a), where the interactions are fixed at a ii s /a ii dd = 0.74 and a 12 s / a 11 dd a 22 dd = 0.86, while the density is the same as Fig.6(a).Panels (a) and (b) present the 3D density profiles of ρ1 and ρ2, respectively, via the slice plots.To clearly see the geometry of such states, the 2D distributions ρ1(z = zc) and ρ2(z = 0) are shown in panels (c) and (d) as well.Here zc denotes the position where ρ1 reaches its maximum in the z direction.
This work was supported by the National Nature Science Foundation of China (Grant No.: 12104359), National Key Research and Development Program of China (Grant No.: 2021YFA1401700), Shaanxi Academy of Fundamental Sciences (Mathematics, Physics) (Grant No.: 22JSY036), and the Danish National Research Foundation through the Center of Excellence "CCQ" (Grant No.: DNRF156).Y.C.Z.acknowledges the support of Xi'an Jiaotong University through the "Young Top Talents Support Plan" and Basic Research Funding as well as the High-performance Computing Platform of Xi'an Jiaotong University for the computing facilities.F.M. acknowledges the funding from the Ministerio de Economía y Competitividad (PID2021-128910NB-100).
2D for the two components.ρ 2 is the component that remains localized around z = 0 and ρ 1 gets pushed out in the immiscible regime.That is reflected by the fact that ρ 1 is a piece-wisely defined