Pattern Formation in Quantum Ferrofluids: from Supersolids to Superglasses

Pattern formation is a ubiquitous phenomenon observed in nonlinear and out-of-equilibrium systems. In equilibrium, quantum ferrofluids formed from ultracold atoms were recently shown to spontaneously develop coherent density patterns, manifesting a supersolid. We theoretically investigate the phase diagram of such quantum ferrofluids in oblate trap geometries and find an even wider range of exotic states of matter. Two-dimensional supersolid crystals formed from individual ferrofluid quantum droplets dominate the phase diagram at low densities. For higher densities we find honeycomb and labyrinthine states, as well as a pumpkin phase. We discuss scaling relations which allow us to find these phases for a wide variety of trap geometries, interaction strengths, and atom numbers. Our study illuminates the origin of the various possible patterns of quantum ferrofluids and shows that their occurrence is generic of strongly dipolar interacting systems stabilized by beyond mean-field effects.

Here, we study single-component quantum ferrofluids confined in cylindrically symmetric geometries, including beyond mean-field quantum fluctuations. We find that extending the geometry from 1D to 2D in a trapped system extends not only the crystal structure of the droplet phase to the second dimension, but also gives rise to new morphologies. We show that dipolar BECs have a remarkably rich phase diagram as we find quantum liquid states of matter, including supersolid honeycomb and superglass labyrinthine morphologies beyond the supersolid droplet morphology.
In Sec. II, we briefly review our methodology and give an overview of the interactions in quantum ferrofluids. We present the ground-state phase diagram of quantum ferrofluids for an oblate trap geometry in Sec. III and discuss the types of morphologies, their location in the phase diagram, and the origin of the pattern formation (morphogenesis). In Sec. IV, we show that dimensionless units reveal scaling properties of quantum ferrofluids in the presence of quantum fluctuations and discuss the geometry dependence of the patterns. The scaling relations generalize the phase diagram discussed in Sec. III to a wide range of trap geometries and allow to tune the strength of the stabilization mechanism of the morphologies. Furthermore we show that, by simply adjusting the trapping confinement, geometric transitions between BEC, honeycomb, labyrinthine, and droplet states are possible. The characteristic length scale of the patterns follows the same scaling with trapping geometry that is known from the roton momentum of dipolar BECs and extends it to new and unexpected regimes. Finally, we conclude in Sec. V and provide an outlook of our study.
We consider a cylindrically symmetric harmonic trap V ext (r) = M ω 2 r (x 2 + y 2 + λ 2 z 2 )/2 with aspect ratio λ = ω z /ω r and the mass M of the atomic species.
The contact and dipolar interaction strengths g s = 4πh 2 a s /M and g dd = 4πh 2 a dd /M . These quantities are determined by the scattering length a s and the dipolar length a dd = µ 0 µ 2 m M/12πh 2 with the magnetic moment µ m . The long-range and anisotropic dipolar interaction with the dipoles aligned by a magnetic field along theẑ-direction is given by U dd (r) = (3/4π)(1 − 3z 2 /r 2 )/r 3 [11].
In the following, we are most interested in the ground states of the system for parameters where structured forms of matter arise. To understand structure formation as a result of competing interactions [1], we consider the underlying energy contributions of a state described by the eGPE in the context of a density functional theory [15,[104][105][106][107]. The eGPE can be formulated as ih∂ t ψ = δE/δψ * [108], where the right hand side is the functional derivative of the energy functional with respect to ψ * . We find ground states by a direct minimization of Eq. (2) using conjugate gradient techniques [101,106,109,110]. We denote the density n = |ψ| 2 and the integrands of Eq. (2) as an energy density E. Equation (2) contains the repulsive contributions by the contact interaction E con ∝ g s n 2 and quantum fluctuations E qf ∝ g qf n 5/2 , which importantly have a distinct scaling with the density [13,54,57,100,111]. The dipolar interaction is longrange and anisotropic and can give an attractive contribution E dd < 0 for particles that arrange in a head-totail configuration. The competition between the attractive dipolar and repulsive contact interaction can therefore lead to mean-field instabilities that are stabilized by the stronger density scaling of the quantum fluctuations. Repulsive and attractive interactions at different length and density scales are the key components in Eq. (2) that lead to structure formation and are also present in other systems such as optically coupled cold atoms [30-32, 92, 112], nuclear matter [38,111], helium droplets [15,16], and colloidal systems [1,3,105,113,114]. In the context of cold atomic physics, strongly dipolar BECs represent a realistic system holding the potential for complex pattern formation in equilibrium [12,13,92,100].

III. PATTERNS IN QUANTUM FERROFLUIDS
Here we first discuss the various morphologies that can be found in the phase diagram of quantum ferrofluids in oblate traps. Second, we turn to the origin of the pattern formation, the morphogenesis.
We consider a strongly dipolar BEC of 162 Dy atoms (a dd 130 a 0 ) confined in a cylindrically symmetric oblate trap with trapping frequencies ω/2π = (125, 125, 250) Hz, aspect ratio λ = 2, and a magnetic field alongẑ. The phase diagram for the chosen parameters is connected by scaling relations to similar phase diagrams in other trap geometries or with other atomic species as we show in Sec. IV.
In order to gain insight into the pattern formation of quantum ferrofluids we map out the ground-state phase diagram in a wide range of interaction strengths and atom numbers around the instability boundary from a BEC to structured states of matter, as shown in Fig. 1. We search for the ground state at every scattering length and atom number by setting a random initial wavefunction [115], allowing us to avoid hysteresis effects when crossing phase boundaries in parameter space [74].
The boundary below which the BEC transitions to structured phases is described by a critical scattering (a) The left-hand side show the N -as phase diagram for trap frequencies ω/2π = (125, 125, 250) Hz. The right hand side shows 2D density cuts n(x, y, 0) at relevant points in the phase diagram, shown by the corresponding markers. The density distributions for a specific marker are ordered in atom number from left to right. The BEC at high atom numbers has a ring of depleted density near its boundary (circles) and forms honeycomb structures toward smaller scattering lengths (diamonds). The honeycomb structures persist to higher atom numbers (squares) and move outwards to the rim of the density distributions, while the core of the BEC spatially saturates in density. An example of the pumpkin state can be seen in Fig. 2. Dashed lines indicate crossovers between different regions. (b) The transition between droplets and honeycombs occurs via stripes (b2) that break up into droplets at small scattering lengths (b1). At high atom numbers and low scattering lengths (b3-b6), labyrinth structures form that are almost degenerate with many other morphologically different labyrinth structures. The supersolid droplets form density connections toward higher atom numbers and transition to labyrinthine structures. The field of view for the 2D-densities in (a) and (b) is 14 × 14 µm 2 .
length a s,c . The structured states have a reduced symmetry compared to the rotationally symmetric BEC state, as the continuous rotational symmetry is spontaneously broken for scattering lengths below a s,c . The spontaneous rotational symmetry breaking characterizes the appearance of supersolid or superglass phases, where crystalline or amorphous spatial structure coexists with superfluid flow [53]. We find that the BEC can transition to a variety of patterns, namely supersolid droplet (SSD), honeycomb and stripe or labyrinth phases [12,13,92,116]. The phase diagram is shown in Fig. 1 As shown in Fig. 1(a) (circles), the BEC states near a s,c develop a radial substructure such that they differ from a Thomas-Fermi density distribution. The BEC states in the range N 60-200 × 10 3 near a s,c show a ring of depleted density near their boundary in addition to the maximum density in the center of the trap ( Fig. 1(a), circles, left column). At intermediate atom numbers (N 200-400 × 10 3 ) a second minimum in the center of the trap can occur and toward higher atom numbers, the trap center is filled with atoms and only the depleted density ring near the boundary remains ( Fig. 1(a), circles, right column). A special case of the BEC shape occurs toward lower atom numbers (N < ∼ 50×10 3 ), where the maximum density in the center of the trap vanishes, leaving only the density ring away from the trap center. These states are known as biconcave or blood cell states [89,94,101,[117][118][119][120][121][122][123][124][125][126] due to the similarity to the shape of a red blood cell. Indirect experimental evidence of these shapes has recently been found [89] and a theoretical study explained their connection to supersolid droplets by investigating elementary excitations across the transition [91].
The honeycomb phase ( Fig. 1(a), diamonds and squares) forms for sufficiently high atom numbers with a s < a s,c , where density bridges connect the central maximum and the outer ring. When another density minimum is present in the center of the trap, multiple rings with connecting density bridges and honeycomb patterns with six, seven or more density minima form. These structures feature strong density connections, facilitating superfluid flow along the honeycomb pattern [12,74,79,92,127]. In combination with the crystalline structure that develops, these states form a supersolid phase [29,100]. Comparing the three-, four-, six-droplet states (stars) with the three-, four-, six-minima honeycomb states (diamonds) shown in Fig. 1(a) suggests that there is a symmetry between positive droplets and negative droplets on top of a background density distribution. In the infinite quasi-2D system [92], it was shown that this is indeed a symmetry where the honeycomb structure becomes energetically favorable over the hexagonal droplet crystal beyond a critical density. We find that a similar symmetry exists in the harmonically trapped finite size system we consider here ( Fig. 1(a), stars and diamonds). The region in which the change from droplet to honeycomb occurs is determined by an interplay between the overall density and the quantum fluctuation strength [92].
In a window of atom numbers where the BEC-SSD boundary changes to the BEC-honeycomb boundary, the transition below a s,c can occur via stripes ( Fig. 1(b), b 2 ) or honeycomb patterns deforming into stripes toward smaller a s . The emergence of the stripe phase between supersolid droplets and honeycomb phases has been observed with Quantum Monte Carlo simulations [12] and in a mean-field theory in a scenario where three-body interactions ∝ n 3 [100] take the stabilizing role instead of quantum fluctuations ∝ n 5/2 [54,57,59]. We have con- Morphogenesis. Shown are density distributions of the ground state with N = 1000 × 10 3 atoms. At such high atom numbers, the BEC features a density saturated core and develops a density depleted ring in its crust when the scattering length as is reduced to its first boundary (as = 88.7 a0). For smaller as, the outer high-density ring becomes unstable and breaks in, yielding the pumpkin-like state (as = 88.6 a0). As as is reduced further, the depleted density wanders closer to the core and is closed off by an outer shell of density, yielding a dense core surrounded by a honeycomb structure (as = 88.0 a0). When the depleted density has expanded yet closer to the core a new depleted density ring forms (as = 86.2 a0) which again forms honeycomb structures (as = 86.0 a0). As the stabilizing repulsion becomes insufficient to uphold the fine density bridges of the honeycomb pattern, some of the connections break up and yield a labyrinthine pattern (as = 84.0 a0). The top and bottom row show 3D density distributions and 2D density cuts n(x, y, 0), respectively.
firmed that toward larger aspect ratios, yielding larger samples (toward the thermodynamic limit), the intermediate stripe phase is enlarged in the phase diagram [115]. When a s is further reduced, these stripes break up their connections and reenter the supersolid droplet phase (b 1 ). However, toward higher N and smaller a s , these stripes can curve and form overlap with neighboring stripes, representing a small region in the larger labyrinthine phase ( Fig. 1 This labyrinthine phase consists of elongated and curved density stripes. The amorphous spatial structure together with the strong density connections, supporting superfluid flow along the labyrinthine stripes, classify the labyrinth as a superglass. In the labyrinthine regime ( Fig. 1(b), b 3 -b 6 ) we cannot unequivocally determine the true ground state by a random initial wavefunction or by choosing a previously found low-energy state, since we find for fixed N and a s many morphologically distinct labyrinthine patterns that are almost degenerate [5,6,8,93,94,96,100], with total energy differences of a few single Hz per atom. However, we find the labyrinth states to be robust against small perturbations [6,8,93,96], be it in changes of scattering length or trap deformations.
With these observations about the morphologies, we now turn to the important change occurring in the phase diagram of Fig. 1(a), namely that the critical scattering length a s,c changes from rising to falling with increasing atom number. Qualitatively, the shape of the phase boundaries in Fig. 1(a) can be understood by noting that Eq. (2) contains the three distinct scalings ∝ n (single-particle), ∝ n 2 (mean-field), and ∝ n 5/2 (quantum fluctuations) [13,57,59]. While the phase diagram for low atom numbers is dominated by stabilization due to quantum pressure (kinetic energy) [108,128], the interplay between mean-field interactions and quantum fluctuations determines where a s,c rises quickly with atom number ( Fig. 1(a)). For a high density, the stabilizing quantum fluctuations dominate and allow for a smaller contact repulsion with the same effective stabilization, hence the phase boundaries (including a s,c ) decrease with atom number [92,115]. This change of a s,c coincides with a peak density saturation in the ground state distributions as the honeycomb and labyrinthine phases appear for a s < a s,c [115]. A saturating density is a defining feature of self-bound and isolated quantum droplets [13,[57][58][59][60][61][62][66][67][68][69], which develop a flat-top (spatially saturated) density distribution toward high atom numbers. The saturation signals an increasingly quantum liquidlike behavior and reduced compressibility compared to the BEC state, like for a liquid compared to a gas. Similarly for the honeycomb and labyrinthine phases, the observation of a saturating density leads to an intuitive understanding of the morphogenesis.
The effect of a saturated density in the ground state distributions for the morphogensis is best understood by following a BEC state at a high atom number through the various instability boundaries toward smaller a s , as shown in Fig. 2. Toward the atom number shown in Fig. 2, the BEC close to a s,c grows and develops a shellstructure reminiscent of ultra-dense neutron stars [37][38][39]129]. In the study of neutron stars, the occurrence of stable and nonuniform states of matter below the saturation density in the crust of the stars is known as nuclear pasta [37,39]. Analogously as seen in Fig. 2, the dense "core" of the BEC is saturated and the density varies spatially mostly in the "crust" of the BEC. Quantum fluctuations stabilize the core and prevent crystallization by an increasing density. Instead the system minimizes its energy by depleting density toward smaller a s . The first stage of this behavior is presented by the depleted density ring occuring in the crust of the BEC due to the inward pressure provided by the external harmonic trap (Fig. 2, a s = 88.7 a 0 ). The atom number determines how close to the boundary of the BEC this depletion occurs. Toward higher N , the core region of the BEC grows and the depleted density ring shifts outwards. While the BEC-honeycomb transition is crossed toward smaller a s up to around N 700 × 10 3 (cf. Fig. 1(a)), for N > ∼ 700 × 10 3 the depleted ring is located so close to the boundary (cf. Fig. 2), that an instability similar to the fingering instability known from classical ferrofluids occurs at a s,c [2, 3, 5-7, 9, 10]. The BEC at these high atom numbers passes through an intermediate state when a s is reduced, which we call the pumpkin state (Fig. 2, a s = 88.6 a 0 ). Toward smaller a s the repulsive contact interaction and quantum fluctuations become weaker and destabilize the core region as transitions through honeycomb to labyrinthine states occur, as detailed in Fig. 2, by a cascade of depleted density rings that form holes and wander closer to the core region.
One can connect the decrease in a s,c and the associated morphologies for a s < a s,c to the infinite system case [92]. In the infinite system, the decrease happens roughly above a critical density at which the three phases of BEC, droplet and honeycomb are connected by a second-order phase transition [130]. Generically below or above this critical density, the BEC is connected by a first-order transition to the honeycomb or droplet patterns in the infinite 2D system [92]. Consistent with the observations in the infinite system, here in the finite size system we find that the transition from BEC to the stripe states around the turning point of a s,c occurs more smoothly with no clear jump in peak density between N 60×10 3 and N 110 × 10 3 compared to the transition from BEC to the supersolid droplet or honeycomb phase at lower or higher atom numbers, respectively.
The morphogenesis of supersolid droplets for a s < a s,c at low atom numbers (cf. Fig. 1(a)) is a special case as the system can minimize its energy by locally increasing density with the crystallization of supersolid droplets, which are not density-saturated. Studied in detail recently [89,91], their morphogenesis is explained by the softening of elementary excitations called angular roton modes near a s a s,c , which provide an angular instability and split the rotationally symmetric BEC structure into droplets.

IV. SCALING PROPERTIES OF QUANTUM FERROFLUIDS
The pattern formation studied above is by no means the outcome of fine-tuning of system parameters. Indeed, here we show that they are generic features of a phase diagram for dipolar quantum gases that can be discussed using dimensionless parameters and scaling laws.
We note that the ground state solution of Eq. (2) is uniquely specified by the external potential parameters {ω i } and the interaction parameters (a s , a dd , N ). In our present case, the external potential parameters correspond to the trap frequencies of the harmonic confinement, but may be left general in case of other external potentials.
Since Q only explicitly depends on C, N , and on the ratio D/C through dd , a generalization of the phase diagram ( Fig. 1(a)) to different atomic species is straightforward. For a fixed trap geometry, we base the length unit on the dipolar length x s = 4πa dd and obtain (C, D) = ( −1 dd N, N ) [68]. Therefore Q is only a function of dd and N . Consequently in a fixed trap, the only parameters determining the type of morphology are the atom number N and the relative dipolar strength dd and, for the trap discussed in Sec. III, the phase diagram generalizes to different atomic species by replacing the a s -axis with −1 dd for any given a dd . For a fixed atomic species in varying cylindrically symmetric traps, choosing x s = h/M ω r (therefore ω 0 = ω r ) is useful as this choice leaves only the aspect ratio λ = ω z /ω 0 as an independent parameter for the external trapping potential V ext (r) = (x 2 + y 2 + λ 2 z 2 )/2. In this formulation, Eqs. (4)-(5) reveal that the contact and dipolar interaction strengths C ∝ D ∝ N √ ω 0 follow the same scaling with atom number and trap frequency. Therefore E nl (C, D, 0) is scale invariant when N √ ω 0 is kept constant [132] and quantum ferrofluids in Shows how the broken scale invariance can be used to relate similarities between different parameter regimes as an example for a trap frequency modified by a factor of two. The transition from BEC to honeycomb for a scaled atom number and correspondingly scaled trapping frequency ( √ 2N, ω0/2) occurs at higher scattering lengths due to the effective reduction in Q, which is compensated by larger as. The dimensionless densityñ = nx 3 s /N at z = 0 is shown in (a) and (b) to compare the ground state density in both trap geometries.
the absence of quantum fluctuations obey an important scaling property. Once a solution for a certain (C, D) is known, an entire family of solutions with higher atom numbers and smaller trapping frequencies or vice versa has been found [101,123,124,132]. In the presence of quantum fluctuations (Q > 0), the scale invariance is broken due to the explicit atom number dependence of Q ∝ C 5/2 /N . Therefore the strength of the stabilizing quantum fluctuations can effectively be tuned along the contours N √ ω 0 = const. Such scaling properties have also proven useful for BECs interacting with an induced gravity-like interaction [133,134] and one-dimensional systems [135,136], where they enabled the reduction of the parameter space dimension by one. In our case, the scaling behavior of Q along the contours N √ ω 0 = const.
allows to tune the strength of the stabilization mechanism of the structured quantum ferrofluid states of matter, as we show in the following. In Fig. 3(a), we illustrate the utility of tuning the quantum fluctuations in a quantum ferrofluid for the example of a honeycomb state. We take a four-minimum honeycomb ground state (cf. Fig. 1(a)) and vary the parameter Q by a few percent to understand the effect of this scaling on the ground states. We see that changes in Q and a s are similar since both provide a repulsive and stabilizing effect, only with a different density scaling. To this end one may note that Cn 2 + Qn 5/2 = C(r)n 2 acts as an effective contact interaction, with a spatially dependent scattering length whose spatial dependence is given by C(r) = C + Qn(r) 1/2 . Figure 3(b) shows how this scaling can be realized by reducing the trapping frequencies by a factor of two while keeping the aspect ratio λ = 2 constant. Due to the reduction of the stabilizing quantum fluctuations in lower confinements the BEC-honeycomb transition has shifted to higher a s . Therefore at the same scattering length as in the higher confinement, the state in the lower confinement is already in the droplet regime with a s = 89.5 a 0 . Toward this scattering length, the ground state in lower confinement has transitioned from the honeycomb phase through a stripe phase and finally to the supersolid droplet regime. Intuitively, Q (and similarly a s ) controls the tendency of the density in the ground state to bond with nearby density structures. Therefore the reduced Q leads to structures that bond less, the droplet state being the result of a labyrinthine state losing its tendency to bond. Generally, similar (C, D, Q) provide an efficient way to locate similar phases in the parameter space of the energy functional parametrized by the physical quantities (a s , a dd , N ). In particular the scaling between atom number and trap frequencies suggests that the quantum liquid states of matter shown in Figs. 1-2 might be observable in more tightly confined traps at experimentally accessible atom numbers [91,115,137], provided that loss mechanisms are negligible and a high optical resolution is available to resolve these fine structures. With higher trap frequencies (smaller x s ) one may trade off the benefit of well-separated structures (larger x s ) for similar ones with enhanced quantum fluctuations at smaller atom numbers. While here we only considered trap aspect ratios of two, these arguments are also valid for different cylindrically symmetric traps [115,138], as we show in the following.
An interesting property of quantum ferrofluids derives from the anisotropy of the dipolar interaction, which is their geometry dependent stability [122,[138][139][140][141]. The tunability of the trapping frequencies allows to investigate this geometry dependent stability continuously both theoretically and experimentally. Above, we showed that an overall scaling of trapping frequencies can be absorbed into the dimensionless interaction strengths. In the following we investigate how the morphologies are influenced by the only independent geometric parameter in the system -the aspect ratio λ = ω z /ω r . There is a difference between changing the aspect ratio by modifying ω z with constant ω r and vice versa since the magnetic field alongẑ breaks the symmetry between the radial and axial directions. Two cases arise, namely either a change in vertical confinement or radial confinement, as we show in Fig. 4 and Fig. 5, respectively. Figure 4(a) shows that the ring-state in a nearly spherical trap transitions to the BEC purely by a geometric change of the trapping confinement. The state transitions through the labyrinthine phase, an increasingly macroscopically developed honeycomb phase and finally a pumpkin state. The patterns become finer as the vertical confinement increases (Fig. 4(a)). Analogous to the situation in classical ferrofluids confined be- The insets in the lower right corners show the spatial power spectrum (PS) Sn(kx, ky) in arbitrary units. The crystallinity can be seen from the diffuseness of the PS along the ring with radius |k| = k * . Labyrinthine states have a powdered (diffuse) PS at k * , reflecting the amorphous or glassy density distribution [116,142]. Toward honeycomb states, the PS concentrates in a triangular pattern indicating the increasing crystallinity. The pumpkin state (λ 3.8) PS shows more angular peaks at k * corresponding to its higher discrete rotational symmetry. (b) The characteristic momentum at radial wavevector |k| = k * scales as k * ∝ 1/lz ∝ √ ωz and defines the characteristic spacing of the morphologies 2π/k * ∝ lz, where lz = h/M ωz is the harmonic oscillator length along the magnetic field direction (vertical direction). A least-squares fit to k * /2π = c/lz as a function of vertical confinement yields c = 0.206(2) (lr = h/M ωr 0.71 µm). Doubling N or changing as by 1 a0 yields a similar behavior with a deviation of c by less than 2% (see main text).
tween two plates [1,[6][7][8], the higher vertical confinement frustrates the morphologies more strongly and leads to their thinning. The spatial power spectrum (PS) S n (k x , k y ) = |F[n(x, y, 0)](k x , k y )| 2 , shown in the insets of Fig. 4(a), reveals information about how many length scales are involved in the morphologies, the crystallinity, and the spacing (fineness) of the structures. We have denoted F[g](k) = g(r)e ik·r d 2 r as the Fourier transform of a function g. Since the states have no modulation along z the PS of the cut suffices to analyze the structures. The PS is concentrated radially around a single characteristic momentum |k| = k * . This single radial concentration shows that there is only a single characteristic length scale in the morphologies, corresponding to 2π/k * . The spacing (fineness) of the structures can be seen in the absolute value of k * as a function of vertical confinement. Figure 4(b) reveals that the spacing scales as 2π/k * ∝ l z , where l z = h/M ω z is the harmonic oscillator length along the magnetic field direction.
This scaling behavior is known from the roton momentum k rot , defining the characteristic momentum at which the dispersion relation of a dipolar BEC shows a distinct roton minimum [55,56,143]. The collective excitations associated to this minimum, the roton modes, are precursors to a structural phase transition when the roton minimum softens near zero excitation energy. Representing the dominant fluctuations driving this transition [81,89,91], the roton modes carry their length scale, the roton wavelength λ rot = 2π/k rot , over into the newly emerging ground state structure and provide its characteristic structural length scale 2π/k * . The fact that the characteristic length scale across a structural phase transition can be interpreted to originate from softening or energetically low-lying excitations on the highersymmetry-side of the transition is a generic result of linear stability analysis in nonlinearly interacting systems, such as classical ferrofluids [5,6] or nonlinear optics [26,27,[29][30][31] and is therefore general beyond the situation in quantum ferrofluids [105,107]. In the supersolid droplet regime, this behavior has been thoroughly studied recently [78,79,81,82,91]. Figure 4 shows that this scaling behavior persists from the BEC state to the honeycomb phase, throughout the multistable labyrinthine phase to the ring state.
Relating the domain spacing to the roton momentum suggests that the coefficient c = 0.206(2) for the characteristic momentum k * /2π = c/l z mostly depends on chemical potential and maximum density in the system [55,56,143,144]. As the density is saturated for the labyrinthine and honeycomb phases, the chemical potential varies weakly with atom number in these regimes. Therefore, c varies weakly with atom number and yields a robust characterization of the fineness of the structures for a given interaction strength and trap geometry. We have repeated the analysis shown in Fig. 4 with a different scattering length a s = 84 a 0 and atom numbers N = {700, 1000} × 10 3 and find that c varies by less than than 2% at these different parameters. Toward lower atom numbers, the peak density and chemical potential become more sensitive to interaction parameters and trapping frequencies and c is generally a function of these parameters. However, the scaling k * ∝ 1/l z remains. Figure 5 shows the behavior of the morphologies with decreasing radial confinement with a fixed vertical confinement ω z /2π = 250 Hz for the same a s and N as in Fig. 4. Instead of a transition from labyrinthine phase through honeycomb and pumpkin states toward the BEC (Fig. 4), one finds in Fig. 5 that the labyrinthine phase loses its density connections and transitions into the crystalline droplet phase. The PS (insets in Fig. 5) shows that the characteristic momentum k * does not change during the transition. These observations can be understood as follows.  Fig. 4. Across the labyrinthine to supersolid droplet transition, the characteristic momentum k * /2π 0.43 µm −1 defining the length scale of the phases stays roughly constant up to λ = 5. Toward higher λ, k * weakly increases. The most significant change is that the weight of the PS (insets) concentrates into a triangular pattern, presenting the emerging multiple Brillouin zones of the macroscopic crystalline pattern with a lattice constant 2π/k * formed by the droplets seen in position space. As the labyrinthine patterns lose some density connections, they transition first to slightly noncylindrical droplets (λ 3.1) for which the PS still is slightly diffuse on the ring with |k| = k * toward the pristine crystal at larger aspect ratios (λ 5). Reducing the radial confinement or increasing the vertical confinement (see Fig. 4) both increase the aspect ratio, but the effective change in the morphologies is drastically different between the two cases.
Equations (4)-(6) with x s = h/M ω r show that a decreasing radial confinement leads to a reduction in the dimensionless interaction strengths similar to a decreasing atom number [115]. In the phase diagram ( Fig. 1) this decrease corresponds to a crossing of the labyrinthine-SSD boundary at constant a s , explaining the labyrinthine to supersolid droplet transition seen in Fig. 5. Since decreasing ω r additionally leads to an increase of the natural length scale x s ∝ 1/ √ ω r and aspect ratio √ λ ∝ x s , the transition is not exactly equivalent to a change in atom number but corresponds to a trajectory through four-dimensional parameter space (C, D, Q, λ) [115]. As the spacing of the structures at constant (C, D, Q) decreases as 2π/k * ∝ 1/ √ λ (Fig. 4), but for the case of decreasing ω r the natural length scale x s ∝ √ λ expands at the same rate, these two effects roughly balance and lead to a constant k * .
Finally, we note that a change in the aspect ratio combined with a change in atom number according to λ → ∞, N → ∞, n 0 = const. corresponds to systems approaching the thermodynamic limit [55,56,143,[145][146][147][148]. Accordingly one expects quantum ferrofluids to form more macroscopic structures toward larger aspect ratios. Repeating the calculation for the phase diagram toward larger aspect ratios, we indeed find that the structures become more macroscopic and that the morphologies discussed in Sec. III prevail [115].

V. CONCLUSION AND OUTLOOK
In conclusion we identify new quantum liquid forms of matter in quantum ferrofluids beyond the supersolid droplet regime. We have shown a general phase diagram of quantum ferrofluids in an oblate trap, which features supersolid droplets at low densities and labyrinthine, honeycomb, and pumpkin states toward higher densities. The emergence of these morphologies can be traced back to the increasingly dominant role of quantum fluctuations toward higher densities, providing the underlying stabilizing mechanism. The strength of this stabilization can be tuned by adjusting the overall trapping confinement. Due to the anisotropy of the dipolar interaction, the morphologies can be transformed into one another by a simple adjustment of the trap aspect ratio. Squeezing the quantum ferrofluid morphologies along the magnetic field direction reveals that the characteristic length scale of the morphologies follows the same scaling behavior as the roton wavelength known from ordinary BEC states.
The labyrinthine states hint at a large degeneracy of the ground state within the framework of an effective mean-field description. This calls for a more elaborate theory beyond the effective description in this labyrinth phase, which however is beyond the scope of the current work. In particular, an interesting possibility is that the various labyrinthine morphologies we find to be degenerate in our effective description might actually be selected upon by quantum fluctuations [149].
Another direction worth investigating is to obtain further insight into the dominant collective excitations giving rise to the honeycomb and labyrinthine morphologies. A linear stability analysis similar to studies on the BEC to supersolid droplet transition [78,79,82,91] may allow identification of modes characteristic of the supersolid or superglass nature of these patterns.
We anticipate that an extension of our study to molecules [137,150] with tunable electric dipole moments could reveal further interesting phases in regimes where strong correlations and the granular nature of matter play an important role [105,114,151].
Note added. Upon submission of the present work, we became aware of a related and very recent preprint [152].

A. Phase diagram
We find ground states using conjugate gradient techniques [101,106,109,110]. When searching the ground state from a random initial wavefunction for the phase diagram shown in Fig. 1, we use gradient noise [153]. We typically find a faster convergence with gradient noise compared to white noise or gaussian states. The mean-field dipolar potential is effectively calculated using Fourier transforms, where we use a spherical cutoff for the dipolar potential. The cutoff radius is set to the size of the simulation space such that there is no spurious interaction between periodic images [101,123,154].
In order to understand the behavior of the morphologies towards the thermodynamic limit, we recalculate the phase diagram (see Fig. 1(a)) for an aspect ratio of λ = 3 by keeping ω z /2π = 250 Hz constant and reducing the radial trapping frequencies to ω r /2π = 83.3 Hz. The droplet, labyrinth, honeycomb and pumpkin phases can be found in the new phase diagram as well and the relative location of the phase boundaries are similar to Fig. 1(a). Examples of the droplet, stripe, honeycomb and pumpkin states for an aspect ratio of λ = 3 are shown in Fig. S1. Compared to smaller aspect ratios, the morphologies have expanded and become more macroscopic, as expected. More droplets, stripes, honeycomb minima and fingers of the pumpkin state form. The boundary described by a s,c has shifted to higher atom numbers and scattering lengths, and the rate is smaller with which a s,c decreases toward higher atom numbers above the critical atom number N 200 × 10 3 . The shift of the boundaries in lower radial confinements can be intuitively understood by considering that quantum fluctuations are reduced (see Sec. IV), and therefore higher scattering lengths and atom numbers are required to reach similar patterns. The fact that in different trap geometries the overall structure of the phase diagram is similar, in particular that the superglass and supersolid states of matter prevail, shows that no fine-tuning of the trap geometry or atom numbers is necessary to observe these structures. The morphologies are not fine-tuned states but rather phases of matter in the complex phase diagram of quantum ferrofluids. Furthermore, the scaling relation provided in Sec. IV provide an intuitive understanding of the changes induced on the boundaries between these phases by an overall change in trapping frequencies or the atom number.
For completeness of the discussion regarding the peak density in the main text, we show the peak density n 0 along vertical and horizontal cuts of the phase diagram in the main text ( Fig. 1(a)) in Fig. S2. Figure S2(a) shows that the peak density has a jump at small atom numbers, when the supersolid droplet regime is entered and that the discontinuity becomes smaller toward higher atom numbers. Above N 120 × 10 3 , where the honeycomb phase separates the BEC phase from the labyrinth, the critical scattering length a s,c decreases with increasing atom number. In the labyrinthine phase, there are fluctuations in peak density as the scatttering length is varied since labyrinths with different forms can have slightly different peak densities. However, toward smaller scattering lengths they follow the same general functional form n 0 (a s ) regardless of the atom number, which indicates that the peak density is saturated in the labyrinthine phase. The saturation can also be seen in Fig. S2(b) as a function of atom number for fixed scattering lengths. In the BEC regime (a s = {90, 90.5} a 0 ) the peak density rises relatively quickly up to N 120 × 10 3 , where the behavior of the critical scattering length a s,c changes, as described in the main text. Slightly below this atom number, the stripe phase appears as an intermediate region between the honeycomb and the supersolid droplet (SSD) phase (see Fig. 1). Above N 120 × 10 3 in the BEC regime, the peak density grows significantly slower compared to the initial increase and only weakly depends on the scattering length (cf. Fig. S2(a)). Spatially, the core region of the BEC close to the honeycomb transition is roughly density saturated (see Fig. 2) and grows slowly when the atom number is increased as shown in Fig. S2(b). One can see from the peak density with a s = 89 a 0 , where the honeycomb phase is entered and exited as a function of atom number, that in the honeycomb phase the peak density is saturated and when it is exited, follows the same behavior of the BEC. At smaller scattering lengths, one can see that the peak density still grows in the droplet regime when the atom number is increased, but is saturated in the stripe and labyrinthine phases.
A change of the peak density behavior can be seen in the BEC phase close to the instability boundary where the atom number is high enough to support a density maximum in the center of the trap, surrounded by a ring of depleted density near the boundary ( Fig. 1(a), circles). Increasing the atom number from there on mainly leads to an overall growth of the BEC structure while maintaining the depleted density ring near its boundary. The overall peak density in the BEC still grows for higher atom numbers, but at a smaller rate [115].
Towards smaller scattering lengths in the honeycomb  and labyrinthine phases, the depleted density needs to redistribute itself among the remaining density connections, leading to a moderate increase of density in these structures. This process leads to the transition from honeycomb to labyrinthine states, as some density connections weaken sufficiently toward low scattering lengths to break up. The states in the honeycomb and labyrinthine phases do not increase their peak density toward higher atom numbers but only grow in size and change their morphology. Along the density lines in the honeycomb and labyrinthine phases, the density remains spatially almost flat. Honeycomb and labyrinthine phases with their macroscopically saturated density distribution realize a quantum liquid that is even further extended in space than the previously studied isolated and self-bound quantum droplets [13,47,54,[57][58][59][60][61][62][63][65][66][67][68][69]. Despite their saturated density n 0 , these phases are still ultradilute compared to strongly interacting or ordinary liquids [63,65] as the gas parameter in the entire phase diagram of Fig. 1(a) stays below n 0 a 3 s < ∼ 3 × 10 −4 . There is an analogy to the situation of elongated supersolids [115], where a decrease in critical scattering length beyond a critical density has also been noticed when quantum fluctuations are included in the description [75,80,84]. In elongated geometries, the transition from BEC to droplets is smooth in an intermediate density regime [80] and the decrease in critical scattering length is only observed when the LHY correction is included [84]. Beyond the intermediate atom number regime, it is seen that the density modulation forms first around the outer boundary and moves inwards for decreasing a s [75]. Analogously in the oblate trap, we find that the density modulation at low atom numbers occurs in the center of the trap (blood cell and subsequent droplet formation) and at higher atom numbers the density minimum of the blood cell moves outwards and bridges form between the central maximum and the outer ring, yielding the honeycomb phase, as shown in Fig. 1. The smooth transition in an intermediate atom number regime can also be seen in the round trap as we discussed in the main text, as well as the infinite system [92].
We point out a similarity to self-assembling collodial systems with competing interactions [105,114]. In these systems, the phase diagram has a similar basic structure as shown in Fig. 1 [105], and minimum-energy configurations can show the formation of shell structures [114], in which multiple depleted density regions form in a honeycomb pattern.
We conclude with final remarks regarding the phase diagram on the transition between the honeycomb and labyrinthine phases.
When setting a previously found ground state across the honeycomb-labyrinthine transition, we find that the discrete rotational symmetry of the honeycomb state persists to smaller scattering lengths compared to the state one finds when searching from a random initial state. In some cases, honeycomb structures can lose their outer connections toward smaller scattering length which yields droplets that surround circularly symmetric density rings. When comparing these states with stripe and labyrinth states that we obtain by searching for the ground state from a random initial wavefunction, we find that the ring states are typically a few Hz up to ten Hz higher in total energy per atom, which indicates that these are metastable states originating from a hysteresis. Nonetheless we find that these states may be relevant for future work, as we performed real-time simulations where we slowly ramp the scattering length across the transition and find that these metastable ring states can be long-lived (we have evolved these states up to 120 ms after the ramp is complete and find them to be stable) and might therefore be experimentally observable.

B. Reduced units and scaling properties
As described in the main text, we define the dimensionless variablest = tω 0 ,r = r/x s ,ψ = ψ x 3 s /N to nondimensionalize the energy functional [68,91,92,123,124,131]. Here, ω −1 0 and x s can at first be taken as arbitrary quantities with units of time and length, respectively. For Schrödinger-like equations it is convenient to define the energy and time units as =h 2 /M x 2 s and ω −1 0 = M x 2 s /h based on the unit of length x s . A significant consequence is that the contact and dipolar interaction terms (Eqs. (4)-(5)) for a given atomic species are only ever modified by the product N √ ω 0 .
This is a result of the contact and dipolar interaction both being quadratic in ψ and since the threedimensional convolution (∝ d 3 r) of the dipolar interaction (U dd ∝ 1/r 3 ) stays invariant when scaling space by a factor of x s . With other nonlocal interactions following different power law behaviors, for instance in systems with an induced gravity-like interaction (∝ (1/r) * |ψ| 2 ), different scaling properties with the atom number can