Stability of a flattened dipolar binary condensate: emergence of the spin roton

We develop theory for a two-component miscible dipolar condensate in a planar trap. Using numerical solutions and a variational theory we solve for the excitation spectrum and identify regimes where density- and spin-roton excitations are favored. We characterize the various instabilities that can emerge in this system over a wide parameter regime and present results for the stability phase diagram. Importantly this allows us to identify the parameter regimes where a novel roton-immiscibility transition can occur, driven by the softening of the spin roton excitation.

We develop theory for a two-component miscible dipolar condensate in a planar trap. Using numerical solutions and a variational theory we solve for the excitation spectrum and identify regimes where density-and spin-roton excitations are favored. We characterize the various instabilities that can emerge in this system over a wide parameter regime and present results for the stability phase diagram. Importantly this allows us to identify the parameter regimes where a novel roton-immiscibility transition can occur, driven by the softening of the spin roton excitation.

I. INTRODUCTION
Two-component Bose-Einstein condensates in which one or both of the components have a large magnetic moment, present a new class of superfluids for exploring the interplay of mixture physics and long-ranged dipole-dipole interactions. Fantastic recent experimental progress has led to a number of possible realizations. Most notably, the production of Bose-Einstein condensates of two different isotope mixtures of the highly magnetic atoms Er and Dy [1,2], and the demonstration of magnetic Feshbach resonances to control the short ranged intra-and inter-species interactions of these mixtures [3]. Also by suppressing dipolar relaxation it may be possible to realize multi-component superfluids using several different spin states of a single isotope [4]. Another possibility involves a mixture of a highly magnetic atom with a weakly-or non-magnetic atom (cf. related work on degenerate fermionic mixtures [5,6]). The case of a two-component (binary) magnetic superfluid has been the subject of theoretical proposals for a new class of quantum droplet states [7][8][9] and supersolid phases [10][11][12][13]. The interplay of immiscibility and the long-ranged interactions has been the subject of investigations exploring pattern formation and novel instabilities [14,15]. Also, a number of studies have considered properties of the condensate ground states in stationary [16] and in rotating frames [17][18][19][20] (also see [21]).
The binary magnetic superfluid has a large number of parameters, including three contact interaction strengths, and the moments and orientation of the magnetic dipoles underlying the long-ranged dipole-dipole interactions (DDIs). This presents a rich landscape for exploring new behavior. In recent work the miscibility and stability of a dipolar condensate mixture was studied in a harmonic trap [16]. This work found that criteria derived from a homogeneous theory provided a good quantitative prediction of the instabilities for various regimes, with the general observation that the miscible regime was reduced as the strength of the DDIs increased. However, in pancake geometry traps (with the magnetic dipoles polarized along the tightly confined direction), the miscible regime extended over a much wider range of interaction parameters in quantitative disagreement with the homogeneous predictions. These results motivate the need for theoretical understanding of the physics of such a mixture in a flattened or planar trap.
Here we present theory for a binary dipolar condensate confined in a planar trap, allowing for the dipole polarization of the components to vary, including the extreme cases of parallel and anti-parallel dipoles (see Fig. 1). Our theory is developed for the case of the miscible phase where the two components overlap and are uniform in the plane of the trap. We solve for the quasiparticle excitation spectrum to study the collective modes of the system, and explore the conditions under which roton-like excitations emerge from the momentum dependence of the DDIs in the presence of confinement. These rotons can be of density, (pseudo) spin density, or mixed character, as can be revealed from looking at the density or spin-density structure factors. We use the excitation spectrum to locate the instabilities of the system, identified by excitations becoming dynamically unstable. These instabilities can be categorised by various features, including whether the unstable excitation causes density or spin-density fluctuations, and whether it is a phonon (long wavelength) or roton (short wavelength). From this analysis we obtain stability phase diagrams, identifying where the various types of instabilities arise as a function of system parameters. We also develop various analytic estimates to identify the boundaries of these instabilities, thus providing a useful characterisation of the system over a wide parameter regime.
The outline of the paper is as follows. In Sec. II we present the main formalism of the paper for the uniform miscible ground states and excitations. We also specialize the formalism to the balanced case, where the areal density and intraspecies contact interaction of each component is taken to be the same. This simplifies the treatment of the system and affords a variational approximation. Results for the balanced system are presented in Sec. III. Here our focus is on the stability phase diagrams for uniform miscible ground states and identifying the character of the instabilities that dominate in various parameter regimes. This analysis is aided by criteria developed from the limiting behavior of the interactions. We also consider how the roton or spin-rotons manifest in the density or spin-density structure factors. In Sec. IV we generalize our consideration and present results for the unbalanced case, before concluding in Sec. V.

II. FORMALISM
Here we consider a planar binary Bose-Einstein condensate of magnetic atoms with axial harmonic confinement of angular frequency ω z and no confinement in the transverse plane (see Fig. 1). We develop our theory under the assumption that the two atomic components have the same mass M and experience that same confinement potential. This should be a good approximation for a mixture of any two bosonic dysprosium or erbium isotopes, where the mass varies by less than 5% and the ac-polarizability is similar (e.g. see [1]).

A. Meanfield theory
The meanfield theory for this system is provided by the twocomponent dipolar Gross-Pitaevskii equation (GPE) [22]. Focusing on a miscible phase we take the component wavefunctions to be translationally invariant in the transverse plane, i.e. Ψ i (x) = √ n i ψ i (z), where i = 1, 2 is the component label, and n i is the areal component density with dz|ψ i | 2 = 1.
The magnetic dipoles of each component are taken to be polarized along the axisê =ẑ cos α +x sin α, that is at an angle of 0 ≤ α < π 2 to the z-axis and lying in the xz-plane (see Fig. 1). The DDI potential is of the form where g dd ij = µ 0 µ m i µ m j /3 is the DDI coupling constant between atoms in components i and j, with µ m i being the magnetic moment of component i on theê-axis. Note that we also consider the case of anti-parallel dipoles (see Fig. 1) where µ m 1 > 0 and µ m 2 < 0, for which g dd 12 < 0 and g dd ii > 0.
With these choices the GPE for component i takes the form where µ i is the associated chemical potential. The GPE operator is given by where g s ij = 4πa s ij 2 /M is the s-wave coupling constant describing the short ranged interactions between components i and j, and we have introduced the interaction potential in k-space 1

B. Excitations
We can quantify the excitations of the condensate by solving for the Bogoliubov-de Gennes (BdG) quasi-particles. The BdG equations can be obtained by linearizing the timedependent GPE i Ψ j = L j Ψ j about a stationary solution as where ψ j is a solution of the time-independent GPE (2) and ϑ j (x, t) is taken to be a small fluctuation of the form Here ρ = (x, y) is the planar position vector, ν labels the excitation band, and k ρ = (k x , k y ) is the planar momentum vector. The expansion (7) introduces the quasi-particle modes u νkρ (z) = (u νkρ1 (z), u νkρ2 (z)) and v νkρ (z) = (v νkρ1 (z), v νkρ2 (z)), with excitation frequency ω νkρ , where c νkρ is the (assumed small) amplitude of the mode. The BdG equations that must be solved to determine the quasi-particles take the form where w νkρ = (u νkρ , v νkρ ) T and which contains two sub-matrix operators and Here kρ = 2 k 2 ρ /2M and where F z denotes the one-dimensional Fourier transform of z coordinate, and F −1 z denotes the inverse transform.

C. General numerical treatment
For the numerical results presented in this paper we follow the methods shown in [23] closely. Importantly, we make use of a truncated interaction potential [23,24] to reduce finite grid range effects in the evaluation of the DDIs. For other aspects related to solving for the ground state and excitations we use the techniques described in [16,25].

D. Spin and density structure factors
The dynamic structure factor is a measure of the density fluctuations in a system, and can be used to characterize its structure and collective excitations. For a two component system the dynamic structure factor can be generalized to characterize the density (+) or spin-density (−) fluctuations in the plane and is defined as (see Ref. [16]) where δn ± νkρ ≡ δn νkρ1 ± δn νkρ2 , with being the integral of the density fluctuation of component j over z. The density dynamic structure factor has been measured in dipolar quantum gases using Bragg spectroscopy (e.g. see [26,27]). The dynamic structure factor, and the static structure factor defined by are particularly sensitive probes to roton excitations which we consider in the results sections.

E. Balanced regime
In this work we initially focus on results for the balanced regime where the component densities and intra-species interactions are identical, i.e. n = n 1 = n 2 , g s 11 = g s 22 , and |µ m 1 | = |µ m 2 |. This latter condition means that g dd 11 = g dd 22 = |g dd 12 |, but allows for two possible orientation classes: parallel dipoles (i.e. µ m 1 = µ m 2 ) with g dd 12 = g dd ii , which we denote as ↑↑; and anti-parallel dipoles (i.e. µ m 1 = −µ m 2 ) with g dd 12 = −g dd ii , which we denote as ↑↓.
The balanced configuration allows us to concentrate on how the excitations of the system depend on the effect of the interspecies contact interaction g s 12 , the relative strength of the DDIs to the contact interactions, and the parallel versus antiparallel orientation of the dipoles.
Under the balanced conditions, when the system is miscible, the ground state profile is the same for both components, we set ψ j → ψ with µ j → µ. The GPE (2) reduces to where g p 0+ denotes the effective condensate interaction with p = {↑↑, ↑↓} denoting the two dipole orientations: Notably, the effect of the DDIs on the ground state cancel for the anti-parallel case. For the balanced case the excitations can be further chosen to have the components fluctuating in-phase (λ = +1) or outof-phase (λ = −1), i.e.
where σ x = [ 0 1 1 0 ] is the x Pauli matrix, i.e. u ±νkρ = (u ±νkρ (z), ±u ±νkρ (z)), etc. The BdG equations are then block diagonal and can be decoupled into the two uncoupled 2-by-2 systems: where where we definẽ noting that in the balanced caseŨ λ (k) =Ũ ii (k) + λŨ 12 (k), and we discuss the extension to the unbalanced case in Sec. IV C.
The excitations with λ = +1 describe (total) density fluctuations whereas those with λ = −1 describe (pseudo) spin fluctuations (also see discussion of the dynamic structure factors in Sec. II D). For this reason we will refer to the relevant excitations as density and spin branches, andŨ + (Ũ − ) as the density (spin) interaction.

F. Properties ofŨ λ
Many of the excitation properties we examine can be understood from the limiting behavior ofŨ λ (k). In this subsection we collect the key results that we use later.
The long wavelength behavior is characterized by This quantity is independent of k z [see Eqs. (18) and (19)], and g p 0+ corresponds to the condensate interaction. The full set of values for g p 0λ , including those for λ = −1, are given in Table I.
It is convenient to expressŨ λ in the form 2 With this identification we have introduced (also listed in Table I). In cases where g p ∆λ = 0 (i.e. the density interaction for the anti-parallel case, and the spin interaction for the parallel case), the interactions are momentum independent, and thus contact like. For cases where g p ∆λ = 0 we have to carefully look at the finite momentum excitations for the development of roton-like excitations.
It is also useful to consider the interactions for shortwavelength in-plane excitations. In practice for our choice of tilting dipoles in the xz-plane the most stringent conditions for stability (i.e. least repulsive interactions) occurs along the k y -axis, and thus we define so that g p ∆λ = (g p 0λ − g p yλ )/3, with values for the particular cases evaluated in Table I.

G. Variational treatment
A variational Gaussian treatment can be applied to the balanced case. In this approach each component wavefunction has the form where σ is the dimensionless variational axial width parameter and l z = /M ω z . The total energy of the condensate is given by where γ ≡ dz l z |ψ σ | 4 = 1/ √ 2πσ. We can use the variational approximation to describe the lowest in-phase and the lowest out-of-phase excitation bands (i.e. ν = 0 and λ = ±1) by employing the "same-shape approximation" [23,28,29], i.e. taking that We can then integrate out the z coordinates to obtain where [cf. Eq. (23)] with G 0 (q) = √ πqe q 2 erfc(q).

A. Instabilities for the parallel case
First we consider the instabilities that can occur for the case of a balanced binary dipolar fluid with parallel dipoles. We show three example stability phase diagrams for this system in Figs. 2(a1)-(a3). The various instabilities are identified from the type of excitation that becomes dynamically unstable (i.e. the excitation develops an imaginary energy).

Density instabilities and the density roton
The phonon (long-wavelength) interaction for this system, g ↑↑ 0+ , must be non-negative for the relevant modes of the system to be stable. We identify g ↑↑ 0+ = 0 as the phonon instability boundary (also see discussion of component stability below).
Because g ↑↑ ∆+ = 0 the density interaction is momentum dependent for the parallel system. If g ↑↑ y+ < 0 then the high-k y interaction is attractive and this can allow a roton excitation to form [e.g. see Fig. 2(b)]. For g ↑↑ y+ sufficiently attractive the roton excitation goes soft, marking the onset of a roton instability. The condition g ↑↑ y+ = 0 thus provides a lower bound for where a roton can form. In practice g ↑↑ y+ has to be sufficiently attractive for the interactions to overcome the kinetic energy of the excitation, allowing a local minimum to form in the dispersion. Thus, the precise conditions for where the rotons form have to be determined from a calculation of the excitation spectrum. Such a spectrum is shown in Fig. 2(b) for a case close to where the roton goes to zero energy. The wavevector of the roton k roton is identified as the location of the local minimum in the lowest density-branch excitation (i.e. ω kρ0+ ). We numerically determine the roton boundary by searching for the interaction parameters where the roton energy is zero. The blue and orange lines in Figs. 2(a1)-(a3) show the roton instability boundary from the BdG solution and variational theory, respectively. Figure 2(c) shows the roton wavevector along the roton instability boundary and Fig. 2(d) gives the interaction parameter coordinates of the roton instability boundary. The groundstate depends only on g ↑↑ 0+ , and the matrix H + (22) depends only on g ↑↑ 0+ and g ↑↑ y+ , so the density roton boundary g ↑↑ y+ , can be written as a universal function of only g ↑↑ 0+ . We also note that comparison of numerical calculations of the BdG equations to the variational results in these figures shows good qualitative agreement.
The density modes of the parallel balanced system map onto those of an equivalent scalar dipolar gas (in the same planar geometry) with the identification of g s ii + g s 12 and 2g dd ii as the effective contact and dipolar coupling contacts of the scalar system, respectively. The basic density instabilities of this system are thus identical to those of the scalar system explored in Ref. [23]. Importantly, the roton stability boundary is universal in the units of Fig. 2(d) (c.f Fig. 5(a) of Ref. [23]), and is the same curve appropriately mapped with the choice of α and differing g s 12 values in Figs. 2(a1)-(a3). Of course the additional features of the binary dipolar gas contained in the spin structure and the individual components (discussed below), add additional structure to the stability phase diagram.

Spin instabilities
The spin instabilities can be analyzed in a similar way to the density instabilities, but using the relevant results from theŨ − interaction. However, as g ↑↑ ∆− = 0, there is no momentum dependent behavior in the spin interaction. Thus we only need to consider the long wavelength behavior characterized by g ↑↑ 0− . The condition g ↑↑ 0− = 0 marks the boundary to a long wavelength instability. Because this occurs in the spin interaction, the instability marks the development of immiscibility (phase separation) of the two components, and we refer to this instability as phonon immiscibility.

Component stability
For the uniform binary system to be stable, we also require that each individual component is stable. This requires that the coefficient of |ψ i | 2 in Eq. (3) is non-negative. This coefficient corresponds to the long-wavelength limit g ii . Thus we identify the condition for the component-i stability boundary as g ii = 0. This marks the onset of a phonon instability in this component. For the results shown Figs. 2(a1)-(a3) we see that the density and single component stability both contribute to determining the region in which the system has a phonon instability.
The density instabilities (of the phonon or roton form) and the single component instabilities, indicate a mechanical instability of the system where the total density or the density in a component can collapse and form density spikes. In this regime the quantum fluctuation effects are necessary to stabilize the collapse (e.g. see [7][8][9]).

B. Instabilities for the anti-parallel case
We now consider the instabilities that can occur in a balanced binary dipolar fluid where the dipole moments of the two components are anti-parallel. We show three example stability phase diagrams for this system in Figs. 3(a1)-(a3).

Density instabilities
The phonon interaction for this system is g ↑↓ 0+ and we identify g ↑↓ 0+ = 0 as a phonon instability boundary. Because g ↑↓ ∆+ = 0 the density interaction is momentum independent and a (density) roton cannot occur.

Spin instabilities and the spin-roton
The condition g ↑↓ 0− = 0 marks the long-wavelength phonon immiscibility transition. Because g ↑↓ ∆− = 0 the spin inter- action is momentum dependent for the anti-parallel system. If g ↑↓ y− < 0 then a spin-roton excitation can form [e.g. see Fig. 3(b)]. The spin-roton going to zero energy marks the onset of a roton immiscibility transition, i.e. where phase separation is driven by an unstable mode of a finite wavevector k roton . The condition g ↑↓ y− = 0 thus provides a lower bound for where spin-roton can form. Figure 3(c) shows the roton wavevector at the roton instability boundary and Fig. 3(d) gives the interaction parameter coordinates of the roton immiscibility boundary.
We note that for the anti-parallel case, H − depends on g ↑↓ 0+ , g ↑↓ 0− , and g ↑↓ y− , and the spin roton boundary g ↑↓ y− is a function of both g ↑↓ 0± . So, the roton immiscibility boundary does not display the universality found for the density roton in the par-allel case [cf. Fig. 2(d)]. The critical g ↑↓ y− line terminates when g ↑↓ 0+ = 0 and a phonon mode becomes unstable.

Component stability
The single component stability condition for the antiparallel system is g ii ≥ 0, i.e. the same as for the parallel balanced system, as was given in Sec. III A 3. FIG. 4. The density and spin-density dynamic and static structure factors for a balanced system with parallel (a1)-(a3) and anti-parallel (b1)-(b3) dipole orientations. Results in (a1)-(a3) are for the parallel case considered in Fig. 2(b) and (b1)-(b3) are for the antiparallel case in Fig. 3(b). White solid (dashed) lines show the inphase/density (out-of-phase/spin) excitation spectrum for reference. The dynamic structure factors are frequency broadened by a frequency width of ωB = 0.1ωz.
In Fig. 4 we show the dynamic and static structure factors for the parallel and anti-parallel balanced cases.
The parallel case is shown in Figs. 4(a1)-(a3) and corresponds to the spectrum shown in Fig. 2(b), where we observed a density roton feature. The roton portion of the lowest λ = +1 excitation branch makes a strongly weighted contribution to the density dynamic structure factor [ Fig. 4(a1)] and appears as a prominent peak in the density static structure factor [ Fig. 4(a3)].
The anti-parallel case [corresponding to the spectrum in Fig. 3(b)] exhibits a spin-roton. This feature makes a strongly weighted contribution to the spin-density dynamic structure factor [ Fig. 4(b2)] and appears as a peak in the spin-density static structure factor [ Fig. 4(b3)].
In Fig. 5 we present results for the dynamic and static structure factors for two unbalanced cases. In the first case . These results can be considered as a progression from the results of Figs. 4(a1)-(a3), in which the DDIs of the second component are reduced, but the other parameters are adjusted to be close to the roton boundary. As these cases are unbalanced, the theory of Sec. II E is inapplicable and the excitations do not separate into in-phase/density and out-of-phase/spin classes, so the theory of Secs. II A and II B is required.
The insets to Figs. 5(a3) and (b3) show that the axial condensate density profiles for these cases differ, as g 11 > g 22 , so the variational theory, which assumes that both components are the same, is no longer applicable.
The rotons observed are now of mixed density and spin character, and they contribute weight to both the density and spin-density dynamic structure factors. In practice, we identify the character of the roton by the dominant peak of the static structure factors. For example, the case in Fig. 5(a3) is a density-dominant roton (i.e., S + has a higher peak than S − ), whereas the case in Fig. 5(b3) is a spin-dominant roton.

B. Stability phase diagram
A stability phase diagram for a general unbalanced case is shown in Fig. 6(a). We have used g s 12 / g s 11 g s 22 as the vertical axis of the phase diagram, which characterizes the stability for a two-component non-dipolar condensate 3 . Along the horizontal axis we vary the ratio of the dipole components, with the left and right limits being the anti-parallel and parallel cases.
Similar to the procedure for the balanced case, we identify the instabilities from numerical calculations of the excitation spectrum. A significant difference, as discussed in the last subsection, is that for the unbalanced case we need to use the structure factors to identify whether any dynamically unstable modes have dominant density (i.e. unstable regions) or spin character (i.e. immiscible regions). The boundaries are found using a square tracing algorithm [30], and we refine results using bisection. We note that the ground state in the stable miscible regime is always symmetric along z for the parameters we consider in this phase diagram (cf to Fig. 2 of Ref. [16], also see Ref. [31,32]).
In the stable region blue and red shading is used to quantify the dominant roton peak in the S + or S − structure factor, represented by log 10 S peak ≡ log 10 max kρ,± S ± (k ρ ). We apply a negative sign if the peak of S − is greater than the peak of S + , and if there is no peak we set the overall result to zero. Thus the upper red region indicates a spin-dominant roton, whereas the lower blue region indicates a density-dominant roton. These regions occur close to the boundaries to the portions of the phase diagram that are roton immiscible and roton unstable, respectively.
The phonon unstable region occurs in the lower left part of the phase diagram for an anti-parallel mixture. The phonon immiscible region occurs in the upper right part for a parallelmixture. The criteria for these phonon regions is similar to the non-dipolar case 4 , demonstrating that systems with parallel dipoles, the transition to immiscibility is largely unaffected by dipolar physics, while for anti-parallel dipoles the mechanical stability is largely unaffected (also see [16]). 3 A homogeneous binary non-dipolar condensate is unstable to collapse for g s 12 < − g s 11 g s 22 , while for g s 12 > g s 11 g s 22 it is immiscible. 4  In Fig. 6(b) we show the same stability phase diagram using g 12 / √ g 11 g 22 as the vertical axis, emphasising the g 0λ = 0 approximations to the immiscible (upper) and unstable (lower) phonon boundaries.
C. Approximate stability criteria For the balanced system, the interaction potentialŨ λ introduced in Eq. (24) appears in the BdG equations (23), and the k ρ → 0 and k y → ∞ limits ofŨ λ characterize the phonon instability and immiscibility boundaries, and provide lower bounds for the onset of the roton boundaries. In the unbalanced case, the lack of symmetry between the two components, ψ 1 = ψ 2 , means that there is no associated decoupling into in-phase and out-of-phase modes. Nevertheless we have found empirically that theŨ λ -potential given by Eq. (24) is useful for the unbalanced case, where we still interpret λ = 1 and −1 as indicating the effective density and spin interactions, respectively. Following the analysis of the balanced case we can find values for the effective low momentum and high momentum interactions as g 0λ = lim kρ→0Ũλ (k) and g yλ = lim ky→∞Ũλ (k), and these quantities are given in Table I. Results using these limits for estimating the boundaries are shown in Figs. 6(a) and (b). Generally we have found that the analytical results provide close lower bounds when ψ 1 ≈ ψ 2 . For example, g 0+ = 0 accurately predicts the phonon instability boundary, as in this region ψ 1 ≈ ψ 2 [ Fig. 6(c)]. Conversely, although g 0− is qualitatively correct, it is less accurate since ψ 1 ≈ ψ 2 . In particular there is a dip in the component with the highest interaction coupling constant, ψ 1 [see Fig. 6(e)] 5 . For a non-dipolar system, within the Thomas-Fermi approximation, this dip would appear when g 22 < g 12 [see [33], above the g 22 = g 12 curve in Fig. 6(a)] 6 .
As the analytical results are based on interactions only, the agreement between numerical and analytical results is best when interactions are dominant, i.e. ng ij ω z 7 .
We have compared our planar results to those for a system in a three-dimensional harmonic trap. We have chosen ng s ii / ω z l z = 5 in Fig. 6 to approximately match the value of this quantity in Fig. 2(c) of Ref. [16] at µ m 2 = −µ m 1 , and g s 12 /g s ii = 40/140 (using the peak areal density from the trapped system for n). There is good qualitative agreement between these two figures for the predictions of the stable region and nature of the instability, noting that quantitative agreement is not expected as the areal density of the trapped system varies over the phase diagram and across the harmonic confinement.

V. CONCLUSIONS AND OUTLOOK
In this paper we have developed theory for a twocomponent miscible dipolar condensate in a planar trap. Our theory allows us to quantify the various instabilities arising because of phonon or roton modes becoming dynamically unstable. Notably we develop a number of analytic results that characterize the boundaries or provide lower bounds for these instabilities to develop. The roton unstable and roton immiscible regimes emerge as the new features from the DDIs under confinement. We have also explored how the roton excitations driving these instabilities are revealed through the density-or spin-density fluctuations, as characterized by the dynamic and static structure factors.
In this work we do not explore the nature of the new ground states that would emerge from these instabilities. Where the instability is of a density type (i.e. phonon or roton unstable) the system is expected to mechanically collapse or implode. In this case quantum fluctuation effects can stabilize the collapse leading to a stable higher density quantum droplet array or a density supersolid state forming. Where the instability is of the spin-density type (i.e. phonon or roton immiscible), we instead expect phase separation of the two components to occur. For the case of an unstable spin-roton the ground state in this regime may be an immiscible pattern characterized by a finite microscopic length-scale (e.g., striped pattern or two-dimensional crystal). Aspects of roton-immiscibility dynamics and ground states have been explored in a dipolarnon-dipolar mixture in planar systems [10,14]. This is also a regime where a domain-supersolid can develop in which the two components form a series of alternating domains, producing an immiscible double supersolid (e.g. see [12,13]). Interestingly, this can occur at lower densities where the quantum fluctuation effects are not significant and loss rates will be much lower. Our results show that the anti-parallel system might be favourable for studying the roton-immiscibility transition and that it can occur in a broad ranges of cases for a suitable choice of interaction parameters.
Here we have focused on the planar system, however our approach can be extended to the infinite tubular regime that provides a model of the cigar shape traps that have been extensively used in recent dipolar Bose-Einstein condensate experiments exploring rotons and supersolids (e.g. see [34][35][36][37]).