Fingering instabilities and pattern formation in a two-component dipolar Bose-Einstein condensate

TL;DR

This paper investigates fingering instabilities and pattern formation at the interface of a two-component dipolar Bose-Einstein condensate, revealing spontaneous symmetry breaking and complex pattern emergence due to dipole-dipole interactions.

Contribution

It demonstrates the spontaneous breaking of rotational symmetry and pattern formation in dipolar BECs, with analytical estimates of domain sizes and effects of various parameters.

Findings

Fingering instability causes symmetry breaking.

Emergence of frog-shaped and mushroom-shaped patterns.

Pattern characteristics depend on atom ratio, trap ratio, and polarization tilt.

Abstract

We study fingering instabilities and pattern formation at the interface of an oppositely polarized two-component Bose-Einstein condensate with strong dipole-dipole interactions in three dimensions. It is shown that the rotational symmetry is spontaneously broken by fingering instability when the dipole-dipole interactions are strengthened. Frog-shaped and mushroom-shaped patterns emerge during the dynamics due to the dipolar interactions. We also demonstrate the spontaneous density modulation and domain growth of a two-component dipolar BEC in the dynamics. Bogoliubov analyses in the two-dimensional approximation are performed, and the characteristic lengths of the domains are estimated analytically. Patterns resembling those in magnetic classical fluids are modulated when the number ratio of atoms, the trap ratio of the external potential, or tilted polarization with respect to the z direction is varied.

Figures

Equations

$$i\hslash \frac{\partial }{\partial t}{\Psi }_i\left(\mathbf{r}\right)$$

$$U_{ij}\left(\mathbf{r}\right)=\frac{{\mu }_0{\mu }_i{\mu }_j\left[1-3{\left(\hat{d}\cdot \hat{r}\right)}^2\right]}{4\pi r^3},$$

$$U_{ij}\left(\mathbf{r}\right)=\gamma \frac{{\mu }_0{\mu }_i{\mu }_j\left(1-3{\cos }^2\theta \right)}{4\pi r^3},$$

$$\big{|}\psi_{i}\big{|}^{2}=\int\big{|}\Psi_{i}\big{|}^{2}dz,$$

$${\Psi }_i\left(\mathbf{r},t\right)={\psi }_i\left(\mathbf{\rho },t\right)\frac{1}{{\left(\pi {\zeta }^2\right)}^{1/4}}{e}^{-z^2/\left(2{\zeta }^2\right)},$$

$$i\hslash \frac{\partial }{\partial t}{\psi }_i$$

$$U_{ij}^{2D}\left(\mathbf{\rho }-{\mathbf{\rho }}^{\prime }\right)=\gamma \frac{{\mu }_0{\mu }_i{\mu }_j}{3\sqrt{2\pi {\zeta }^2}}\int \frac{d\mathbf{k}}{{\left(2\pi \right)}^2}{e}^{-i\mathbf{k}\cdot \left(\mathbf{\rho }-{\mathbf{\rho }}^{\prime }\right)}h\left(q\right),$$

$$\displaystyle\psi_{i}\left(\bm{\rho},t\right)=e^{-iM_{i}t/\hbar}\Big{[}\phi_{i}\left(\rho\right)$$

$$\hslash \Omega u_i$$

$$-\hslash \Omega v_i$$

$$E_{dd}$$

$$y=\eta \left(x,t\right).$$

$$f\left(y\right)=\left\{\begin{array}[]{ll}0&\mbox{if $y\rightarrow-\infty$};\\ 1&\mbox{if $y\rightarrow+\infty$},\end{array}\right.$$

$$n_1\left(\mathbf{\rho },t\right)$$

$$n_2\left(\mathbf{\rho },t\right)$$

$$n_1\left(\mathbf{\rho },t\right)$$

$$n_2\left(\mathbf{\rho },t\right)$$

$$\frac{\delta E_{dd}}{\delta \eta \left(x,t\right)}$$

$$g\left(k_y\right)=\int dy{f}^{\prime }\left(y\right)e^{-i{k}_{y}y}.$$

$$f^{\prime }\left(y\right)=\frac{e^{-y^2/l_y^2}}{\sqrt{\pi }{l}_y},$$

$$g\left(k_y\right)=e^{-k_y^2{l}_y^2/4},$$

$$\frac{\delta E_{dd}}{\delta \eta \left(x,t\right)}$$

$$\equiv$$

$$\varpi =\sqrt{\frac{F\left(\kappa \right)\kappa +\sigma {\kappa }^3}{m\left(n_1^0+n_2^0\right)}},$$

$$\psi_{i}=e^{-iMt/\hbar}\Big{[}\sqrt{n/2}+u_{i}e^{i\bm{k}\cdot\bm{\rho}}e^{-i\Omega t}+v_{i}^{*}e^{-i\bm{k}\cdot\bm{\rho}}e^{i\Omega^{*}t}\Big{]},$$

$$\hslash \Omega$$

$${\psi }_1$$

$${\psi }_2$$

$$\frac{E}{N}=\frac{{\hslash }^2{K}^2}{8m}+\frac{g_{dd}n}{4}h\left(K\zeta /\sqrt{2}\right),$$

Full text

I Introduction

Fingering instabilities (e.g., the Saffman-Taylor instability and Rayleigh-Taylor instability) lead to complicated pattern formations of an interface instarmp1986 ; instarmp2009 , and ubiquitously occur at the interface between two fluids of different densities and viscosities, such as water suspended atop oil in gravity keiser2017 , Hele-Shaw cells between two plates hscell2012 ; hscell2015 ; hscell2017 , biological systems bio2010 , mushroom clouds from volcanic eruptions and nuclear explosions in the atmosphere, plasma fusion plasma2000-1 ; plasma2000-2 , supernova explosions supernova2002 ; supernova2013 , and the Crab Nebula crab2008 . In addition, domains of magnetic fluids are also known to undergo fingering instability and form complex labyrinthine patterns magfluid1975 ; magfluid1994 . In resemblance to the classical fluids, a system of two superfluids can also exhibit a Rayleigh-Taylor interfacial instability, for instance, a system of a two-component Bose-Einstein condensate (BEC) fibec1 ; fibec2 ; fibec3 ; fibec4 ; fibec5 . A two-component BEC of an atomic gas with strong dipole-dipole interactions (DDIs) exhibits similarities with those in classical magnetic fluids, such as hexagonal, solitonlike, and labyrinthine pattern formations saito2009 .

Recently, a dipolar BEC with strong DDIs was shown to exhibit self-organized crystallization, which is observed in the quench dynamics kadau2016 , and leads to novel self-bound droplet states of superfluids schmitt2016 ; barbut2016 . These emerging phenomena have inspired many theoretical studies to clarify the implicit new physics behind stabilization against collapse and droplet pattern formation xi2016 ; bisset2016_1 ; blakie2016 ; saito2016 ; wachtler2016_1 ; wachtler2016_2 ; bisset2016_2 ; baillie2016 .

In this paper, we consider a two-component BEC with strong DDIs in which the dipole moments of the two components are polarized oppositely. By modulating the strength of DDI, we demonstrate that spontaneous symmetry breaking and fingering instabilities occur at the interface between the two components, which leads to frog-shaped and mushroom-shaped pattern formations. The Bogoliubov dispersion relation corresponding to the fingering instability is analyzed in a two-dimensional approximation. Spontaneous density modulation and domain growth are also shown in the dynamics, along with analytical estimations of the characteristic lengths of the dipolar domains. We observe the formation of droplet patterns when the population imbalance in two components becomes large. Furthermore, a labyrinthine pattern grows as the trap ratio increases and a stripe phase occurs as the tilted angle increases.

This paper is organized as follows. In Sec. II, we formulate the problem. In Sec. III, we demonstrate the symmetry-breaking fingering instability (Sec. III.1) and domain dynamics of the system (Sec. III.2), where a Bogoliubov analysis and an estimation of the characteristic lengths of the dipolar domains are performed. In Sec. IV, we show stationary pattern formation as the number ratio of atoms (Sec. IV.1), trap ratio (Sec. IV.2), and tilted polarization (Sec. IV.3) are varied. In Sec. V, we provide the conclusions of our study.

II Formulation of the problem

A two-component BEC with DDI at zero temperature is described by the macroscopic wave functions $\Psi_{1}\left(\bm{r},t\right)$ and $\Psi_{2}\left(\bm{r},t\right)$ in the mean-field regime, which obey the nonlocal Gross-Pitaevskii (GP) equations given by

[TABLE]

where the atomic mass $m$ is assumed to be the same for both components. The two-body coupling constants are $g_{ij}=4\pi\hbar^{2}a_{ij}/m$, where $a_{ij}$ is the $s$-wave scattering length between components $i$ and $j$. The macroscopic wave function $\Psi_{i}$ is normalized as $\int|\Psi_{i}|^{2}d\bm{r}=N_{i}$, where $N_{i}$ is the number of atoms in component $i$. The DDI has the form

[TABLE]

where $\mu_{0}$ is the magnetic permeability of vacuum, $\mu_{i}$ is the magnetic dipole moment of component $i$, $\hat{d}=\cos\alpha\hat{z}+\sin\alpha\hat{x}$ is the direction of the dipole polarization with $\alpha$ being the angle between $\hat{d}$ and the $z$ axis, $\bm{r}$ is the vector between the two dipoles, and $\hat{r}=\bm{r}/|\bm{r}|$. In our calculations, the atoms are trapped in the same axisymmetric harmonic potential $V=m[\omega_{\perp}^{2}(x^{2}+y^{2})+\omega_{z}^{2}z^{2}]/2$, where $\omega_{\perp}$ and $\omega_{z}$ are the radial and axial trap frequencies, respectively.

In our numerical simulations, we solve the three-dimensional nonlocal GP equations in Eq. (1) using the pseudospectral method with a fast Fourier transform, where the dipolar interaction terms in Eq. (2) are calculated using the convolution theorem. We take 52Cr as the atoms in the two-component system, with $\mu_{1}=6\mu_{B}$ and $\mu_{2}=-6\mu_{B}$, where $\mu_{B}$ is the Bohr magneton. Experimentally, such a two-component dipolar BEC system can be implemented by using ${}^{7}S_{3}$ $m_{J}=-3$ and $+3$ states of 52Cr crbec2005 ; crbec2006 , and the $s$-wave scattering lengths $a_{ij}$ can be modulated by using Feshbach resonance lahaye2007 .

III Instabilities and dynamics

III.1 Fingering instabilities

We now investigate the fingering instabilities at the interface between two components with different strengths of dipolar interactions with $\alpha=0$, i.e., there is no tilted polarization in this case. To investigate the dependence of the DDI strength on the dynamics, we introduce the strength coefficient $\gamma$ of the DDI as

[TABLE]

where $\cos\theta=\hat{d}\cdot\hat{r}=\hat{z}\cdot\hat{r}$. Experimentally, the coefficient $\gamma$ can be tuned by fast rotation of the magnetic field gamma . We initially prepare stationary states without DDI, which become metastable states with unstable interfaces in the presence of DDI. The DDI is then suddenly introduced at $t=0$, followed by time evolution. The parameters are taken to be $a_{11}=100a_{B}$, $a_{22}=1.1a_{11}$, $a_{12}=1.2a_{11}$, $\left(\omega_{\perp},\omega_{z}\right)=2\pi\times\left(100,800\right)$ Hz, and $N_{1}=1\times 10^{6}$, $N_{2}=2\times 10^{6}$. These scattering lengths satisfy the phase-separation condition, $g_{11}g_{22}<g_{12}^{2}$, and the ground state without DDI has a circular interface between the two components.

The column density profiles obtained by the integration of the densities along the $z$ direction are

[TABLE]

where $i=1$, $2$. These column densities are shown in Figs. 1 and 2 for $\gamma=0.18$ and $1$, respectively. In the absence of the DDI, the immiscible states of the two components exhibit rotational symmetry, as expected, which are shown in Fig. 1(a) and Fig. 2(a). After the DDI is induced, the rotational symmetry of the circular interface is spontaneously broken, leading to a fingering instability at the interface between the two components, which is similar to the interface behavior of classical dipolar fluids magfluid1994 . We find from Figs. 1 and 2 that the system exhibits spontaneous rotational symmetry breaking with a smaller length scale for a larger DDI. For $\gamma=0.18$, a frog-shaped pattern is formed in component 1 at $t=100$ ms, as shown in Fig. 1(g). For $\gamma=1$, a mushroom-shaped pattern first appears at the interface, as shown in Fig. 2(c) and Fig. 2(d), and then the droplet pattern is formed at a later time, as shown in Fig. 2(f) and Fig. 2(g). It is shown that with stronger dipolar interactions, the interface of a two-component BEC becomes more unstable.

For comparison, we then investigate the case in the absence of DDI, in which the scattering lengths are suddenly changed from immiscible to miscible, $g_{11}g_{22}>g_{12}^{2}$ at $t=0$. As shown in Fig. 3, no fingering instability occurs at the interface without DDI (i.e., $\gamma=0$ together with the same values of the other parameters), and the initially phase-separated components mix in the dynamics, which demonstrates that the immiscible-miscible transition has no contribution to the fingering instability at the interface. These results reveal that fingering instability never appears due to the immiscible-miscible transition.

For a deeper understanding of the instability, we perform a Bogoliubov analysis for the system. To simplify the calculation, the dimensionality of the system is reduced to two dimensions (2D). We assume that the $z$ dependence of the wave function is decomposed as

[TABLE]

where $\bm{\rho}=(x,y)$ and $\zeta$ is the width of the density profile in the $z$ direction. The value of $\zeta$ is determined by fitting the Gaussian to the numerically obtained wave function. Equation (1) is then rewritten as

[TABLE]

where $V(\rho)=m\omega_{\perp}^{2}(x^{2}+y^{2})/2$, $a_{z}=\sqrt{\hbar/(m\omega_{z})}$ and $g_{ij}^{2D}=g_{ij}/\small\sqrt{2\pi\zeta^{2}}$. The 2D dipolar interaction reads tiltmean

[TABLE]

where the dipolar kernel $h(q)=2-3\small\sqrt{\pi}qe^{q^{2}}\mbox{erfc}(q)$ with $q=k\zeta/\small\sqrt{2}$ for polarized spins in the $z$ direction.

We define $\phi_{i}(\rho)$ as the stationary state with rotational symmetry, and then weak perturbations of the wave functions can be expressed as,

[TABLE]

where $M_{i}$ is the chemical potential of component $i$, the amplitudes $u_{i}(\rho)$ and $v_{i}(\rho)$ are treated as small, $\varphi=\tan^{-1}(y/x)$, and $\Omega$ is the frequency of the oscillation. The quantum number $L$ characterizes the angular momentum of the modes.

Substituting Eq. (8) into Eq. (6), and collecting the first-order terms proportional separately to $e^{\pm i\Omega t}$, we then obtain the following Bogoliubov de Gennes (BdG) equations,

[TABLE]

where $H_{0}=-\hbar^{2}\nabla^{2}/(2m)+V(\rho)$, and $i=1$, $2$.

We numerically diagonalize Eqs. (9) and (10) to study the instability of the system. If $\Omega$ has a positive imaginary part, the system exhibits a dynamical instability, and the unstable modes will grow exponentially to the evolution limitation arising from the nonlinearities. The imaginary part of the Bogoliubov excitation frequency ${\rm Im}[\Omega]$ as a function of dipolar strength coefficient $\gamma$ is shown in Fig. 4 with modes $L=2-10$, among which the modes with the largest ${\rm Im}[\Omega]$ dominate the instability dynamics. For the $L=0$ and $1$ modes (monopole and dipole modes), ${\rm Im}[\Omega]=0$ and these modes are stable. Each mode in $L=2-10$ has its critical strength of dipolar interaction, above which ${\rm Im}[\Omega]>0$ and monotonically rises. When $\gamma=0.18$ (the vertical black dashed line in Fig. 4) corresponding to Fig. 1, the interface is unstable with modes $L=3-6$, in which the modes of $L=4$ and $5$ with the relatively largest ${\rm Im}[\Omega]$ are predominant in the instability dynamics. We thus see that the results shown in Fig. 1 are in agreement with the Bogoliubov analysis.

We then analyze the fingering instability at the interface induced by the DDI. The DDI energy of the system in the 2D approximation is given by

[TABLE]

where g_{dd}=\gamma\mu_{0}(6\mu_{B})^{2}/\big{(}3\small\sqrt{2\pi\zeta^{2}}\big{)}, and $n_{i}=|\psi_{i}|^{2}$ is the 2D density for component $i$. For simplicity, we assume that components $1$ and $2$ have uniform distributions in $y\lesssim 0$ and $y\gtrsim 0$, respectively. The interface is located at $y\simeq 0$, and the position of the interface is defined by

[TABLE]

A smoothed step function is introduced as

[TABLE]

and $f^{\prime}(y)>0$ is localized at $y\simeq 0$. Using Eqs. (12) and (13), the densities of the two components can be expressed as

[TABLE]

We assume that $\eta$ is small and the densities can be expanded as

[TABLE]

where we assume $f^{\prime}\left(-y\right)=f^{\prime}\left(y\right)$. By neglecting the terms that are independent of $\eta$, we obtain the pressure difference at the interface as

[TABLE]

where

[TABLE]

The function $f^{\prime}\left(y\right)$ is assumed to be Gaussian,

[TABLE]

which gives

[TABLE]

where $l_{y}=\sqrt{\hbar/(m\omega_{y})}$. Assuming that $\eta\left(x,t\right)\propto\sin\left(\kappa x-\varpi t\right)$, we obtain

[TABLE]

which corresponds to the pressure acting on the interface. In the Rayleigh-Taylor instability at the interface between heavy and light fluids in the gravity, the pressure acting on the interface is given by $g(\rho_{\rm heavy}-\rho_{\rm light})\eta$. Thus, $F(\kappa)$ in Eq. (23) plays the role of the difference in the gravitational force on the unit volume of heavy and light fluids in the classical Rayleigh-Taylor instability. According to the Rayleigh-Taylor dispersion relation of the interface mode instabook , we then obtain

[TABLE]

where $\sigma$ is the 2D interface tension coefficient. When $F\left(\kappa\right)\kappa+\sigma\kappa^{3}<0$, the interface mode with wave number $\kappa$ is dynamically unstable. When the healing lengths in the two components are almost the same, $\xi_{1}\simeq\xi_{2}\equiv\xi$, the interface tension coefficient in Ref. tension is approximated as $\sigma\simeq g_{11}^{2D}n_{1}^{0}\xi\sqrt{K-1}$, where $K=a_{12}/\sqrt{a_{11}a_{22}}$. Using the parameters in Fig. 1 and fitting the numerically obtained interface profile to Eq. (20), we obtain $\sigma\simeq 0.015\hbar\omega_{\perp}/a_{\perp}$, and the ratio is $\zeta/a_{\perp}\simeq 1.3$. Then the most unstable wavelength is estimated to be $2\pi/\kappa\simeq 5.5$ $\mu{\rm m}$ for $\gamma=0.18$ and $2\pi/\kappa\simeq 2.2$ $\mu{\rm m}$ for $\gamma=1$, which are in qualitative agreement with those in Figs. 1 and 2. For a more accurate estimation of the unstable wave number, we must consider the circular geometry of the system. Moreover, the interface tension is modified by the DDI, which should be taken into account.

III.2 Domain dynamics

We now investigate the domain dynamics for another initial condition. The initial state is the stationary state $\Psi_{0}$ of a single-component dipolar BEC with the parameters taken as $\gamma=1$, $N=4\times 10^{6}$, $a=100a_{B}$, $\alpha=0$, and $\left(\omega_{\perp},\omega_{z}\right)=2\pi\times\left(100,800\right)$ Hz. Then the stationary state is split into two components with opposite polarization as $\Psi_{1}=\Psi_{2}=\Psi_{0}/\sqrt{2}$. During the real-time evolution of the two components, we take $\gamma=1$, $a_{11}=a_{22}=a_{12}=100a_{B}$, $\left(\omega_{\perp},\omega_{z}\right)=2\pi\times\left(100,800\right)$ Hz, and $N_{1}=N_{2}=2\times 10^{6}$. In this case, the $s$-wave interactions have no contribution to the phase separation of the system because of the equal values of $s$-wave scattering lengths.

The results of the dynamic density modulation in a two-component dipolar BEC are shown in Fig. 5. At the initial time, the phase separation is induced by the DDI, which leads to a complicated domain structure with a large kinetic energy [Fig. 5(b)]. The dipolar domains of the system then grow with time [Figs. 5(b) – 5(e)], and the domain structure is gradually rearranged at later times [Figs. 5(e) – 5(g)]. The domain size in the dynamics is shown as a function of time in Fig. 6. These results demonstrate that the growth of the domain size is similar to the recent results regarding the hydrodynamics in the magnetic domains of a spin-1 ferromagnetic BEC domain2011 .

To obtain a deeper understanding of the domain dynamics, we now perform an analytical estimation of the characteristic lengths both at the beginning and at longer times of the dipolar domain formation process in the 2D approximation. At the beginning of the formation of the dipolar domains, the characteristic lengths are estimated by using a Bogoliubov analysis of a uniform system. We define $g_{11}^{2D}=g_{22}^{2D}=g$, and write the wave functions as

[TABLE]

where $n$ is the 2D total density and $M=n(g+g_{12}^{2D})/2$. Substituting Eq. (25) into Eq. (6) and diagonalizing the Bogoliubov matrix, we obtain

[TABLE]

where $\epsilon_{k}=\hbar^{2}k^{2}/(2m)$. Equations (26) and (27) correspond to the eigenfrequencies of the density and spin waves, respectively. If $\Omega$ is complex, the corresponding mode grows exponentially, and the system is dynamically unstable. Therefore, Equation (27) indicates that the phase separation can occur due to the DDI even when the two components are miscible ($g>g_{12}$) using only $s$-wave scattering. Substituting $g_{12}^{2D}=g$ and the parameters used in Fig. 5 into Eq. (27), the most unstable mode, i.e., the wavelength that maximizes the imaginary part of Eq. (27), is $2\pi/k_{u}=1.05a_{\perp}$ with $a_{\perp}=\sqrt{\hbar/(m\omega_{\perp})}$, which corresponds to the characteristic wave number for the phase separation at the beginning of the time evolution, as shown in Fig. 5(b).

At longer times, for the case with a striped labyrinthine pattern in Fig. 5(g), we apply the trial wave functions of the following form,

[TABLE]

The total energy then becomes

[TABLE]

where the $s$-wave interaction energy is independent of $K$ and neglected. Substituting the parameters used in Fig. 5 into Eq. (30), the characteristic domain size at later times , i.e., the wavelength that minimizes Eq. (30), is $2\pi/K\simeq 2.4a_{\perp}$. Figure 6 shows the time evolution of the average domain size in the dynamics of Fig. 5. The domain size is obtained from the average radius of the excitation ring in the Fourier-space density profiles. The initial and final domain sizes in Fig. 6 are in qualitative agreement with $2\pi/k_{u}$ and $2\pi/K$ estimated above.

IV Domain patterns in stationary states

IV.1 Number ratio

Using imaginary-time propagation, we investigate the stationary states with different numbers of atoms $N_{2}$ in component 2 with fixed $N_{1}$ in component 1. We take $a_{11}=a_{22}=a_{12}=100a_{B}$, $\alpha=0$, and $\left(\omega_{\perp},\omega_{z}\right)=2\pi\times\left(100,800\right)$ Hz. The number of atoms in component 1 is fixed to $N_{1}=2\times 10^{6}$. In all the following simulations, the dipolar strength coefficient $\gamma=1$ is maintained. The $s$-wave scattering lengths herein and in the following are set to fulfill the condition of the miscible states. The ratio between the numbers of atoms $N_{1}/N_{2}$ is changed from $1$ to $10$.

The column density profiles are plotted in Fig. 7. Phase separation shown in Fig. 7(a) arises from the opposite polarization of dipole moments in the two components, and rotational symmetry breaking (labyrinthine pattern) appears spontaneously due to the strong DDI. Moreover, a droplet pattern occurs when $N_{1}/N_{2}>1$, where the droplets tend to form a triangular lattice. The crystallization of the dipolar droplets is similar to that observed in Ref. kadau2016 . The number of droplets decreases as $N_{2}$ is decreased. For component 2 (bottom row in Fig. 7), which has a smaller number of atoms, the filament structure splits into droplets and tends to assemble toward the center. Our results are similar to those in Ref.ryan2012 , in which only one component is dipolar and the system is strictly two dimensional.

IV.2 Trap ratio

We then investigate the trap geometry effect on the pattern formation in the two-component dipolar BEC. The number of atoms in the two components is the same, $N_{1}=N_{2}=2\times 10^{6}$, $\alpha=0$, and $a_{11}=a_{22}=a_{12}=100a_{B}$ are taken. The radial trap frequency is fixed as $\omega_{\perp}=2\pi\times 100$ Hz, while the trap ratio defined by $\lambda=\omega_{z}/\omega_{\perp}$ is changed from $1$ to $16$.

The column density profiles with $\lambda=1$, $4$, $8$, and $16$ are plotted in Fig. 8. As the trap ratio increases with the fixed radial trap potential, the axial trap potential is enhanced, and the atoms in the trap tend to be expanded toward the outside along the radial directions, since the density increases with an increase in the axial confinement. Thus, a more complicated labyrinthine pattern is formed for larger $\lambda$ with the radial trap frequency being fixed. The pattern in Fig. 8(c) is the same as that in Fig. 7(a) due to the same parameter settings. In the imaginary-time evolution, small noise is introduced in the initial state, and the complicated labyrinthine structures obtained after convergence depend on the initial noise. Therefore, a lot of similar structures are almost degenerate near the ground state, and each of them is obtained depending on the initial noise.

IV.3 Tilted polarization

The effect of a tilted polarization on the 2D dipolar Bose gas has recently been studied by mean-field theory tiltmean and dynamic many-body theory tiltpimc . We now investigate the effect of tilted polarization on the interface of a three-dimensional two-component dipolar BEC system. The parameters are taken to be $a_{11}=a_{22}=a_{12}=100a_{B}$, $\left(\omega_{\perp},\omega_{z}\right)=2\pi\times\left(100,800\right)$ Hz, and $N_{1}=N_{2}=2\times 10^{6}$. The angle $\alpha$ of tilted polarization with respect to the $z$ direction is increased from $\pi/12$ to $\pi/3$.

We plot the column density profiles in Fig. 9. The results with $\alpha=\pi/12$, $\pi/6$, and $\pi/3$ are shown. As the tilted angle $\alpha$ increases, the filaments of the labyrinthine structure in Fig. 7(a) are gradually straightened, and the labyrinthine pattern is reorganized into a stripe phase, as shown in Fig. 9(c), which is similar to the quantum phase transition in the 2D case in Ref. tiltpimc . It is shown that the tilted polarization overcomes the effect of labyrinthine instability, and develops an ordered stripe pattern replacing the disordered labyrinthine pattern.

We also modify the number ratio of atoms $N_{1}/N_{2}=6$ in the two components. The column density profiles are plotted in Fig. 10 with $\alpha=\pi/12$, $\pi/6$ and $\pi/3$. The droplet pattern induced by the large number ratio is also reshaped into a stripe pattern as the tilted polarization angle increases. It is shown that a crystallized droplet pattern occurs when $\alpha=\pi/6$ and $N_{1}/N_{2}=6$. Comparing Fig. 10(c) with Fig. 9(c), we find that the number of stripes is decreased as the number of atoms in component 2 is reduced. This is because there are insufficient atoms in component 2 to form many stripes.

V Conclusions

In conclusion, we have investigated the pattern formation of fingering instabilities and domain dynamics in an oppositely polarized two-component BEC with strong DDIs in three dimensions. We have demonstrated that dynamical fingering instabilities occur at the interface between two components of dipolar BECs, which breaks the rotational symmetry of the interface. Frog-shaped and mushroom-shaped patterns have been shown with different strengths of the DDIs. A Bogoliubov analysis gives a qualitative understanding of the interfacial instabilities of the two dipolar BECs, and a dispersion relation similar to that in classical fluids is obtained. Spontaneous density modulation and dipolar domain growth in the dynamics have also been demonstrated, in which we have analyzed the characteristic sizes of the dipolar domains corresponding to different patterns at the initial and later times in the evolution. We have also investigated the parameter dependence of the ground states, and found that the droplet patterns are formed due to the population imbalance in the two components. Labyrinthine patterns grow as the trap ratio increases, and a striped phase appears as the angle of tilted polarization increases.

The findings in this and other recent works have revealed that including anisotropy and long-range dipole-dipole interactions results in rich physics of multicomponent BECs, with similarities to conventional viscous classical fluids. These findings shed light on the nature of dipolar BECs, which should be within reach of current experiments. For example, the system studied in this paper can be implemented by using the ${}^{7}S_{3}$ $m_{j}=-3$ and $+3$ states of 52Cr crbec2005 ; crbec2006 , the $s$-wave scattering lengths can be modulated by Feshbach resonance lahaye2007 , and the strength of the DDI can be tuned by fast rotation of the magnetic field gamma . In this paper, we ignore the dipolar relaxation and spinor dynamics. To observe the phenomena presented in this paper, the application of different dipolar species twospecies may provide a more achievable approach. The experimental resolution of in situ imaging is typically $\sim\mu{\rm m}$, so the patterns of the BECs shown in the present paper can be measured. Recently, it was numerically shown that superfluid flow over a rough surface of a wire exhibits the properties of a boundary layer similar to those in classical fluids stagg2017 , which also establishes further connections between superfluids and classical fluids. Such investigations are likely to be of relevance when examining superfluid flow in practical applications such as atomtronics and quantum metrology seaman2007 ; byrnes2015 ; moxley2016 ; lamico2017 .

Acknowledgements.

K.-T.X. would like to thank E. Braaten for his kind support and comments. H.S. acknowledges support by JSPS KAKENHI Grants No. JP17K05595, No. JP16K05505, No. JP17K05596, and No. JP25103007. T.B. acknowledges support by the Shanghai Research Challenge Fund; New York University Global Seed Grants for Collaborative Research; National Natural Science Foundation of China (61571301); the Thousand Talents Program for Distinguished Young Scholars (D1210036A); the NSFC Research Fund for International Young Scientists (11650110425); NYU-ECNU Institute of Physics at NYU Shanghai; the Science and Technology Commission of Shanghai Municipality (17ZR1443600); and the China Science and Technology Exchange Center (NGA-16-001).