Transverse instability and disintegration of domain wall of relative phase in coherently coupled two-component Bose-Einstein condensates

TL;DR

This paper analyzes the stability and disintegration of domain walls in two-component Bose-Einstein condensates with Rabi coupling, providing an analytical phase diagram and numerical simulations of instability evolution.

Contribution

It presents the first analytical stability phase diagram for domain walls in coupled BECs and demonstrates nonlinear evolution of transverse instabilities through simulations.

Findings

Identified the stability region in the Rabi frequency and coupling constant plane.

Showed the domain wall's transverse instability outside the stable region.

Numerically demonstrated the nonlinear evolution of the instability.

Abstract

We study transverse instability and disintegration dynamics of a domain wall of a relative phase in two-component Bose-Einstein condensates with a coherent Rabi coupling. We obtain analytically the stability phase diagram of the stationary solution of the domain wall for the one-dimensional coupled Gross-Pitaevskii equations in the plane of the Rabi frequency and the intercomponent coupling constant. Outside the stable region, the domain wall is dynamically unstable for the transverse modulation along the direction perpendicular to the phase kink. The nonlinear evolution associated with the instability is demonstrated through numerical simulations for both the domain wall without edges and that with edges formed by the quantized vortices.

Equations

$$\displaystyle E[\Psi_{1},\Psi_{2}]=\int d\bm{r}\biggl{[}\sum_{j=1,2}\biggl{(}\Psi_{j}^{\ast}h\Psi_{j}+\frac{g_{j}}{2}|\Psi_{j}|^{4}\biggr{)}$$

$$\displaystyle+g_{12}|\Psi_{1}|^{2}|\Psi_{2}|^{2}-\frac{\hbar\Omega_{\rm R}}{2}(\Psi_{1}^{\ast}\Psi_{2}+\Psi_{1}\Psi_{2}^{\ast})\biggr{]}.$$

$$i\hslash \frac{\partial {\Psi }_1}{\partial t}$$

$$i\hslash \frac{\partial {\Psi }_2}{\partial t}$$

$$n_0=\frac{\mu +\hslash {\Omega }_R/2}{g+g_{12}}.$$

$$i\frac{\partial {\psi }_1}{\partial t}=-{\nabla }^2{\psi }_1-\tilde{\mu }{\psi }_1+|{\psi }_1{|}^2{\psi }_1+\gamma |{\psi }_2{|}^2{\psi }_1-{\omega }_{\text{R}}{\psi }_2,$$

$$i\frac{\partial {\psi }_2}{\partial t}=-{\nabla }^2{\psi }_2-\tilde{\mu }{\psi }_2+|{\psi }_2{|}^2{\psi }_2+\gamma |{\psi }_1{|}^2{\psi }_2-{\omega }_{\text{R}}{\psi }_1.$$

$${\omega }_{\text{R}}=\frac{\hslash {\Omega }_{\text{R}}}{2g{n}_0},\gamma =\frac{g_{12}}{g},\tilde{\mu }=1+\gamma -{\omega }_{\text{R}}.$$

$$\frac{{\partial }^2\theta }{\partial x^2}=2{\omega }_{\text{R}}\sin \theta$$

$$\theta (x)=4\arctan e^{\sqrt{2{\omega }_{\text{R}}}x}.$$

$$-\frac{{\partial }^2{\psi }_1}{\partial x^2}-\tilde{\mu }{\psi }_1+(1+\gamma )|{\psi }_1{|}^2{\psi }_1-{\omega }_{\text{R}}{\psi }_1^{\ast }=0,$$

$${\psi }_1$$

$$A$$

$${\omega }_{\text{R}}^c=\frac{1}{3}{n}_{\text{min}}({\omega }_{\text{R}}^c)(1-\gamma ).$$

$${\psi }_j={\psi }_j^0+\left[u_j(x)e^{iky-i\omega t}-v_j^{\ast }(x)e^{-iky+i\omega t}\right].$$

$$\displaystyle\hat{\cal H}{\bf u}=\hbar\omega{\bf u},\quad\quad\hat{\cal H}=\left(\begin{array}[]{cccc}\hat{h}_{1}&-\left(\psi_{1}^{0}\right)^{2}&\gamma\psi_{1}^{0}\psi_{2}^{0\ast}-\omega_{\text{R}}&-\gamma\psi_{1}^{0}\psi_{2}^{0}\\ \left(\psi_{1}^{0\ast}\right)^{2}&-\hat{h}_{1}&\gamma\psi_{1}^{0\ast}\psi_{2}^{0\ast}&-\left(\gamma\psi_{1}^{0\ast}\psi_{2}^{0}-\omega_{\text{R}}\right)\\ \gamma\psi_{1}^{0\ast}\psi_{2}^{0}-\omega_{\text{R}}&-\gamma\psi_{1}^{0}\psi_{2}^{0}&\hat{h}_{2}&-\left(\psi_{2}^{0}\right)^{2}\\ \gamma\psi_{1}^{0\ast}\psi_{2}^{0\ast}&-\left(\gamma\psi_{1}^{0}\psi_{2}^{0\ast}-\omega_{\text{R}}\right)&\left(\psi_{2}^{0\ast}\right)^{2}&-\hat{h}_{2}\\ \end{array}\right),$$

$$V_{\text{ext}}=V_0\left[2+\tanh (ar-R)-\tanh (ar+R)\right]$$

$$i\frac{\partial {\psi }_j}{\partial t}=-\frac{{\partial }^2{\psi }_j}{\partial x^2}-\tilde{\mu }{\psi }_j+|{\psi }_j{|}^2{\psi }_j+\gamma |{\psi }_{\bar{j}}{|}^2{\psi }_j-{\omega }_{\text{R}}{\psi }_{\bar{j}}+iV\frac{\partial {\psi }_j}{\partial x}.$$

Full text

I Introduction

Solitons in multidimensional systems are generally unstable, known as the transverse instability Kivshar and Luther-Davies (1998); Carr and Brand (2008), where the solitonic structure is dynamically unstable against the symmetry-breaking modulation along the unbounded spatial dimension. Bose-Einstein condensates (BECs) are well-suited system to study such an unstable properties of solitons. This is because: (i) The properties of BECs are described by the Gross-Pitaevskii (GP) (nonlinear Schöredinger) equation, which allows various solitary wave solutions. (ii) The system consists of dilute gases and is almost far from dissipation so that nonequilibrium processes of unstable dynamics can be directly observed from the well-defined initial states. (iii) The topological solitons can be created by well-developed phase-engineering techniques. Given all three of these advantages, dynamical instabilities of solitons have indeed been observed in dilute BECs. For example, the dark solitons in two-dimensional (2D) BECs are broken up to an array of vortex pairs due to the transverse instability, called “snake instability” Muryshev et al. (1999); Brand and Reinhardt (2002). The decay of a planer dark soliton into a vortex ring was observed experimentally in a 3D BEC Feder et al. (2000); Anderson et al. (2001).

For two-component BECs described by the two-component order parameters, the structures of solitons are richer than those in single-component BECs, e.g., dark solitons Öhberg and Santos (2001), dark-bright solitons Busch and Anglin (2001); Becker et al. (2008), domain walls Coen and Haelterman (2001), and magnetic solitons Qu et al. (2016). For a typical two-component mixture of BECs, the U(1) phases in each component are independent variables, because the two components are coupled only through the density-density coupling. However, when a Rabi coupling is applied between the internal states of two-component BECs, the symmetry associated with one of the two U(1) phases is broken, and the relative phase between the two components makes sense. Then, the new type of soliton can exist as “a domain wall of the relative phase” which is obtained as a kink solution of the sine-Gordon equation for the two-component BECs Son and Stephanov (2002). This domain wall can exist as a bound string between two vortices in two-component BECs, which results in a molecule of the vortices Kasamatsu et al. (2004).

Recently, real time dynamics of the aforementioned vortex molecule in the Rabi-coupled condensates have been studied by some authors Tylutki et al. (2016); Calderaro et al. (2017); Eto and Nitta (2018). The authors in Refs.Tylutki et al. (2016); Eto and Nitta (2018) have especially stressed that this system can be used as the analogous simulations of the “confinement” in quark–anti-quark phenomena. Tylutki et al. studied the precession dynamics of the vortex molecule as a function of the Rabi coupling and the molecular distance Tylutki et al. (2016). In a certain regime, the domain wall was found to disintegrate into some parts, which is analogous to the string breaking phenomena in quantum chromodynamics. This observation was further confirmed by the simulations by Eto and Nitta Eto and Nitta (2018).

Here, we reveal the physical origin of this disintegration as the instability associated with the transverse displacement of the domain wall, i.e, the snake instability. We obtain an exact solution of the domain wall for the 1D coupled GP equations for two-component BECs, making a phase diagram of the stationary solutions in the plane of the intercomponent coupling strength and the Rabi frequency. The energetic stability of the solution reproduces the previous work by Usui and Takeuchi Usui and Takeuchi (2015). In the unstable regime, the dynamical instability takes place for the transverse modulation along the direction perpendicular to the sine-Gordon kink, which is analyzed by the Bogoliubov-de Genne (BdG) equation. The nonlinear dynamics associated with this instability are shown by direct numerical simulations of the 2D GP equation for the domain wall without edges and with edges formed by quantized vortices. The disintegration observed in previous works Tylutki et al. (2016); Eto and Nitta (2018) can be explained as the manifestation of the snake instability.

The paper is organized as follows. After introducing the formulation of the problem in Sec. II, we first study the stability of the domain wall of the relative phase based on the exact solution of the GP equation in Sec. III.1. Next, the transverse instability of the domain wall is discussed through the BdG analysis in Sec. III.2. In Sec. IV we demonstrate the dynamics due to the transverse instability using the direct numerical simulations of the 2D GP equations. Section V is devoted to the conclusion.

II formulation

We consider two-component BECs of ultracold atoms having the same mass $m$ and residing in two different hyperfine states. The BECs are described by the condensate wave functions $\Psi_{j}$ ($j=1,2$). The equilibrium state of the system is obtained by minimizing the Gross-Pitaevskii (GP) energy functional

[TABLE]

Here, $h=-\hbar^{2}\nabla^{2}/(2m)+V_{\text{ext}}$ is the single-particle Hamiltonian with the trapping potential $V_{\text{ext}}$. The coupling constants $g_{j}$ ($j=1,2$) and $g_{12}$ represent the strength of intra- and the intercomponent interactions, respectively, described as $g_{j}=4\pi\hbar^{2}a_{j}/m$ and $g_{12}=4\pi\hbar^{2}a_{12}/m$ with the $s$-wave scattering lengths $a_{i}$ and $a_{12}$ between the corresponding atoms. We assume that the intracomponent coupling constants satisfy $g_{1}=g_{2}=g>0$ for simplicity. The last term in Eq. (1) describes a coherent Rabi coupling induced by an external electromagnetic field, which allows atoms to transfer their internal states coherently Hall et al. (1998); Matthews et al. (1999); $\Omega_{\rm R}$ stands for the Rabi frequency. Then, the total number $N=N_{1}+N_{2}=\int d\bm{r}(|\Psi_{1}|^{2}+|\Psi_{2}|^{2})$ is a constant of motion. In the following analysis except the numerical simulations in Sec. IV.2, we assume the homogeneous system by setting $V_{\text{ext}}=0$. For the coupling constants, we confine ourselves to the range $-1<g_{12}/g<1$ corresponding to a miscible regime without the Rabi coupling, otherwise the components phase separate for $g_{12}>g$ Timmermans (1998); Ao and Chui (1998) or undergo mean-field collapse for $g_{12}<-g$ Adhikari (2001). The properties of a domain wall of the condensate densities for $g_{12}>g$ in the presence of the Rabi-coupling have been studied in detail in Ref. Dror et al. (2011).

The time-dependent GP equations are given by the variational procedure $i\hbar\partial\Psi_{j}/\partial t=\delta E/\delta\Psi_{j}^{\ast}$ as

[TABLE]

First, we consider the ground state in the homogeneous system by ignoring the time and spatial derivative terms in Eqs. (2) and (3). Because of the conservation of the total particle number, the chemical potential as a Lagrange multiplier for both components should be common as $\mu_{1}=\mu_{2}=\mu$, the stationary wave function being written as $\Psi_{j}=\sqrt{n_{j}}e^{i\theta_{j}-i\mu t/\hbar}$. Then, the relative phase $\theta=\theta_{1}-\theta_{2}$ should be vanished, because the energy of the Rabi-coupling term $-\hbar\Omega_{\text{R}}\sqrt{n_{1}n_{2}}\cos\theta$ is minimized at $\theta=0$. As a result, the densities in the miscible regime satisfy $n_{1}=n_{2}\equiv n_{0}$ with

[TABLE]

The miscible regime takes place when the intercompnent coupling constant satisfies $g_{12}<g+\hbar\Omega_{R}/(n_{1}+n_{2})$ Abad and Recati (2013), otherwise the equilibrium state involves spontaneous density imbalance.

We scale the wave function as $\Psi_{j}=\sqrt{n_{0}}\psi_{j}$, and introduce the length, time, and energy scale as $\xi=\hbar/\sqrt{2mgn_{0}}$, $\hbar/(gn_{0})$, and $gn_{0}$, respectively. Then, we get the dimensionless GP equation,

[TABLE]

Here, all the variables are dimensionless and the coefficients are given as

[TABLE]

III Transverse instability of domain wall of relative phase

In the seminal paper Son and Stephanov (2002), Son and Stephanov showed that the GP equation for the Rabi-coupled two-component BECs can be reduced to the sine-Gordon equation when the gradient of the density is neglected. Considering the 1D system along the $x$-axis and substituting the expression $\psi_{j}(x)=e^{i\theta_{j}(x)}$ into Eqs. (5) and (6), we get the sine-Gordon equation

[TABLE]

with the relative phase $\theta\equiv\theta_{1}-\theta_{2}$ and the constant $C$ of integral. The stationary solution for the boundary conditions $\theta=\pi\pm\pi$ and $\partial_{x}\theta=0$ for $x\to\pm\infty$ is given by

[TABLE]

This solution is known as the sine-Gordon soliton Son and Stephanov (2002). The stability criterion was obtained as $\omega_{\text{R}}^{c}=(1-\gamma)/3$ by neglecting the density gradient, which is valid for $\gamma\sim 1$. Son and Stephanov proposed that, when vortices are present in two-component BECs, the Rabi coupling can bind a pair of vortices in the different components via the sine-Gordon domain wall. After that the binding of vortices and composite structures of the resulting vortex molecules have been studied in various situations Kasamatsu et al. (2005a, b); Eto and Nitta (2012); Cipriani and Nitta (2013); Dantas et al. (2015); Aftalion and Mason (2016); Uranga and Lamacraft (2018); Tylutki et al. (2016); Eto and Nitta (2018); Calderaro et al. (2017); Shinn and Fischer (2018).

In this section, we consider the stability of the domain wall of the relative phase based on the GP equations (5) and (6), where we take account of the contribution of the density gradient. We first consider the energetic stability of the domain wall. Although this problem has been considered by several papers Son and Stephanov (2002); Usui and Takeuchi (2015); Qu et al. (2017), we employ the exact solution of the domain wall and extend the phase diagram to the region of the negative intercomponent coupling $\gamma<0$. Next, we consider the transverse instability by extending the 1D domain wall to the additional spatial dimension. The stability can be studied by the BdG analysis, where the signal of the instability is shown by the appearance of the imaginary excitation frequency of the Bogoliubov modes.

III.1 Exact solution and energetic instability

We consider the energetic stability of the domain wall of the relative phase by solving the 1D version of the GP equations (5) and (6). We seek the stationary solutions which satisfy $|\psi_{j}|=1$ at infinity and the phases change continuously as $0\to+(-)\pi$ for $\psi_{1}$ ($\psi_{2}$) from $x=-\infty$ to $x=+\infty$. For the sine-Gordon soliton and our symmetric parameters, the solution should satisfy the relation $\psi_{1}=\psi_{2}^{\ast}$. This restriction reduces the equations to

[TABLE]

which allows us to get an expression of the exact solution of Eqs. (5) and (6) as

[TABLE]

The typical profile of the solution is shown in Fig. 1(a) and (b). The relative phase changes $2\pi$ around the origin over the length scale $\sim(2\omega_{\text{R}})^{-1/2}$. The density profile is written as $n_{1}=n_{2}=\tanh^{2}(\sqrt{2\omega_{\text{R}}}x)+A^{2}\>\text{sech}^{2}(\sqrt{2\omega_{\text{R}}}x)$, being uniform only for $\omega_{\text{R}}=0$. With increasing $\omega_{\text{R}}$, the spatial profile of $\theta$ approaches to a step function and the density depression becomes deeper. The solution of Eq. (11) is effective below the upper bound of the Rabi frequency $\omega_{\text{R}}^{u}=(1+\gamma)/4$, at which the solution coincides with the form of the dark soliton with exactly zero density at the center. Physically, this boundary is interpreted as that the length scale $(2\omega_{\text{R}})^{-1/2}$ is equal to the healing length given by the effective coupling constant $1+\gamma$ [see Eq. (10)].

Figure 1(c) shows the phase diagram representing the energetic stability of the single domain wall in the $\gamma$-$\Omega_{\text{R}}$ plane. Below the line $\omega_{\text{R}}=\omega_{\text{R}}^{u}$, we have the solution of the domain wall of Eq. (11). Here, we check numerically the energetic stability of the solution Eq. (11) through the imaginary time propagation of the GP equations (5) and (6), which is plotted in Fig. 1(c). When the solution is unstable, the density difference grows and unwinds the 2$\pi$ difference of the relative phase to zero, leading to the uniform solution Usui and Takeuchi (2015). The stable range of $\omega_{\text{R}}$ becomes narrower with increasing $|\gamma|$ and vanishes at $\gamma=\pm 1$.

This result agrees with the previous works Son and Stephanov (2002); Usui and Takeuchi (2015). Near $\gamma\sim 1$, the plots are fitted $\omega^{c}_{\text{R}}=(1-\gamma)/3$ obtained in Son and Stephanov (2002). Usui and Takeuchi extended the analysis by considering the density depression through the numerical and variational analyses Usui and Takeuchi (2015). They proposed that the instability is associated with the Landau instability of the local counterflow across the domain wall Son and Stephanov (2002) and the critical velocity was estimated by using the local density at the depression. Since the instability can occur when the wavelength of the unstable excitation is smaller than the length scale of the domain wall, they obtained the expression for the stability as

[TABLE]

Here, $n_{\text{min}}(\omega_{\text{R}}^{c})$ represents the density at the density minimum for $\omega_{\text{R}}=\omega_{\text{R}}^{c}$, being written as $n_{\text{min}}=(1+\gamma-4\omega_{\text{R}}^{c})/(1+\gamma)$ by using Eq. (12). This stability criterion is also shown in the dashed curve in Fig. 1, which is in good agreement with the numerical result.

III.2 Transverse instability of the domain wall

Next, we further study the stability of the sine-Gordon domain wall by the BdG analysis. Especially, we include the fluctuation along the direction ($y$-axis) perpendicular to the sine-Gordon soliton in the $x$-axis to study the transverse instability.

In the standard BdG analysis, the wave function is expanded around the stationary solution $\psi_{j}^{0}$ as

[TABLE]

Here, the fluctuation along the $y$-direction is included by the plane wave $\propto e^{iky}$ Takeuchi et al. (2010); Suzuki et al. (2010). Substituting this expression into Eqs. (5) and (6), we get the eigenvalue equation. The eigenfrequency $\omega$ is calculated by solving the BdG equation

[TABLE]

where ${\bf u}=(u_{1},v_{1},u_{2},v_{2})^{\text{T}}$ and $\hat{h}_{j}=-\partial_{x}^{2}+k^{2}-\tilde{\mu}+2|\psi_{j}^{0}|^{2}+\gamma|\psi_{\bar{j}}^{0}|^{2}$, where $\bar{j}=1(2)$ for $j=2(1)$. Employing the domain wall solution for $\psi_{j}^{0}$, we numerically solve Eq. (19) with the finite system size $-60\leq x\leq 60$ and the 600 grid points. To this end, starting from Eq. (11) which is a solution in an infinite system, we make imaginary time propagation of the GP equation to get the proper solution $\psi_{j}^{0}$ for the finite-size system.

The eigenvalues of the BdG equation can be used to clarify the stability properties of the stationary solutions. Let us first consider the situation $k=0$, where the only 1D perturbation is present. In the stable region in Fig. 1, there are positive eigenfrequencies and one zero-energy mode when we take only eigenmodes with positive norm $\sum_{j}\int dx(|u_{j}|^{2}-|v_{j}|^{2})>0$. When the solution enters the unstable region in Fig. 1, there appears the negative eigenvalue, a signature of the energetic Landau instability. The magnitude of the negative eigenvalue is very small in the most of unstable region, which implies that the Landau instability is very weak. The imaginary time propagation confirms this feature, because the decay of the initial domain wall solution to the uniform one needs very long time. Thus, we expect that this instability is not significant in cold gas experiments at ultralow temperatures.

The more significant instability is the dynamical one. When we include the eigenmodes with finite $k$, there appear imaginary components in the excitation frequencies. Figure 2(a) shows the imaginary part of the excitation frequency as a function of $k$ for $\gamma=0$ and several values of $\omega_{\text{R}}$. The imaginary part appears for the finite range $0<k<k_{\text{max}}$ and takes a maximum value at a certain wave number $k=k_{0}$. The unstable range of $k$ is extended and the maximum of $\text{Im}[\omega]$ is enhanced with increasing $\omega_{\text{R}}$. Thus, the instability becomes stronger with increasing $\omega_{\text{R}}$. In Fig. 2(b), we also plot $k_{0}$ and $k_{\text{max}}$ as a function of $\omega_{\text{R}}$ for $\gamma=0,\>\pm 0.5$. Extrapolation of the curves to $k_{0}=k_{\text{max}}=0$ gives the critical values $\omega_{\text{R}}^{c}$ of the dynamical stability for each $\gamma$, which agrees with the numerically obtained $\omega_{\text{R}}^{c}$ in Fig. 1. We see that the parameter region exhibiting the transverse instability coincides with the unstable region of Fig. 1. For $\omega_{\text{R}}>\omega_{\text{R}}^{u}$ the solution reduces to the dark soliton, where the transverse instability is expected according to the previous literature Muryshev et al. (1999); Brand and Reinhardt (2002). As seen in Fig. 2(b), the plots are connected to the values of $k_{\text{max}}$ for the transverse instability of the dark soliton at $\omega_{\text{R}}=\omega_{\text{R}}^{u}=(1+\gamma)/4$; we have $k_{\text{max}}=\sqrt{(1+\gamma)/2}$ for a dark soliton (see Eq. (7.3) in Ref. Kivshar and Luther-Davies (1998)) following our notation of Eq.(10).

IV Disintegration dynamics of a domain wall of relative phase

In this section, we perform numerical simulations of the GP equations to demonstrate the disintegration of the domain wall of the relative phase through the transverse instability. We consider two cases: (i) an extended domain wall without the edges in a uniform system, and (ii) a domain wall with a finite length in a cylindrical trap. In the latter case, the edges of the wall correspond to the vortices with the same circulation in each component. The disintegration of the domain wall has been found in the numerical simulations in Refs. Tylutki et al. (2016); Eto and Nitta (2018), when the separation of the vortices and the Rabi frequency are large. The transverse instability of a moving domain wall without edges in trapped BECs has been reported by C. Qu et al. Qu et al. (2017).

Typical structures of the domain wall corresponding to the situations (i) and (ii) are shown in Fig. 3. In the contour plots of $\theta$, we show it within the range $-\pi\leq\theta\leq\pi$ for clarity instead of $0\leq\theta\leq 2\pi$, where we have $2\pi$ phase jump at the center of the wall. Then, the wall can be visualized as the localized pattern as shown in Figs. 3(a) and (c), as done in the previous papers Kasamatsu et al. (2004); Tylutki et al. (2016); Eto and Nitta (2018). Figure 3 (a) shows the domain wall extended along the $y$-direction, which gives the initial state of the time development in Fig. 4(b).

This situation (ii) is deeply connected with the experimental observation. A straightfoward way to prepare the domain wall of the relative phase in two-component BECs is that vortices are prepared in each component and then the internal states are coupled by rf fields to induce the Rabi coupling. Then, the vortices are connected with the domain wall of the relative phase. Figures 3 (c) and (d) show the typical structure of the vortex molecule. The $\psi_{1}$-component has a vortex at $(x,y)=(0,22.5)$ with the positive unit winding and the $\psi_{2}$-component has a vortex at $(x,y)=(0,-22.5)$ with the same positive unit winding. When the Rabi coupling is applied, the relative phase between two components tends to becomes zero over the entire space. However, the phase kinks due to the vortices leave a localized structure in the relative phase between the two vortices, the domain wall being naturally formed.

IV.1 domain wall in a uniform system

We first demonstrate the nonlinear dynamics associated with the transverse instability through the simulation of the 2D GP equations (5) and (6) in a uniform system. We first prepare the initial state as $\psi_{j}^{\text{ini}}(x,y)=\psi_{j}(x)$ with the solution of Eq. (11), and calculate the time development using the Crank-Nicholson method. The system size is $[-50,50]$ in the $x$-$y$ plane with numerical grids is 1000 $\times$ 1000. We take the Neumann and periodic boundary condition for $x$- and $y$- direction, respectively. To initiate the dynamical instability, we add a small random noise $\sim 10^{-6}$ to the initial wave functions. We confirm that the domain wall in the stable region in Fig. 1 is certainly stable in the real time development.

Figure 4 shows the typical snapshots of the disintegration dynamics of the domain wall of the relative phase in the unstable regime for $\gamma=0,\>\pm 0.5$. The snapshots show that the small transverse modulation of the wall is amplified after some time and leads to the disintegration of the domain wall. According to the BdG analysis, the growing time scale and the wave length of the unstable excitations are estimated as $\tau\sim 2\pi/\text{Im}[\omega(k_{0})]$ and $\lambda\sim 2\pi/k_{0}$, respectively; for the parameters in Fig. 4, we have $(\tau,\lambda)=(114.9,\>32.2)$ for (a), $(\tau,\lambda)=(247.2,\>64.1)$ for (b), and $(\tau,\lambda)=(149.5,\>46.5)$ for (c). The simulations results are in fairly agreement with these estimations; the quantitative deviations may be due to the finite size effect in the simulations.

After the disintegration, the dynamics of the walls exhibits different behaviors depending on $\gamma$. There, the domain wall with the finite size involves the vortices at the edges. For $\gamma=-0.5$ in Fig. 4(a), the walls rapidly shrink to zero size due to the combined attractive force by the intercomponent coupling and the Rabi coupling between the two components. In other words, the vortices at the wall edges attract each other due the attractive vortex-vortex interaction Eto et al. (2011) as well as the string tension caused by the Rabi coupling. For $\gamma=0$ [Fig. 4(b)], some fragmented walls also tend to shrink due to the attraction by the Rabi coupling, but the others keep their separation. Such vortex molecules undergo center-of-mass motions, going to outside. Eventually, all the walls shrink to zero size at the later stage. Contrary to these, the subsequent dynamics for $\gamma=0.5$ is different. Although the initial wall disintegrates into small pieces, they keep the disintegration-merge cycle for a while, and after that some walls go to outsides. In this case, the repulsive intercomponent interaction $\gamma>0$ prevents the shrink of the wall. In equilibrium, the balance of the repulsive vortex-vortex interaction Eto et al. (2011) and the attractive force by the Rabi coupling realizes the stable vortex molecule Kasamatsu et al. (2004, 2005b); Eto and Nitta (2012); Dantas et al. (2015); Shinn and Fischer (2018).

IV.2 domain wall connecting vortices in a cylindrical trap

We next consider the dynamics of the domain wall having initially a finite length. This situation has been demonstrated in Refs. Tylutki et al. (2016); Eto and Nitta (2018). The authors observed that the domain wall rotates due to the attractive tension between the two vortices by the Rabi coupling. However, when the separation between vortices or the Rabi coupling becomes large, the wall tends to disintegrate. Although, in general, the BdG analysis is not allowed to apply such a non-stationary configuration, we can understand qualitatively that their numerical observations are certainly due to the transverse instability as discussed below.

In simulations, we numerically solve the 2D GP equation with the cylindrical trap

[TABLE]

with the radial coordinate $r$, the potential depth $V_{0}=10$ and the radius $R=50$ of the cylinder. The sharpness of the wall boundary is represented by the parameter $a=1$. We use the trap of Eq. (20) because there is additional contribution to the vortex dynamics from the density inhomogeneity when we employ the harmonic trap Tylutki et al. (2016); Calderaro et al. (2017).

We simulate the time development with the following procedure. First, the initial state without vorticity is prepared for given $\gamma$ and $\omega_{\text{R}}$ through the imaginary time evolution of the GP equation. Next, we imprint a vortex in each component by multiplying the phase factor $e^{i\theta_{v}(\bm{r})}$ with the profile $\theta_{v}(\bm{r})=\arctan[(y\pm y_{0})/x]$ ($-$ for $\psi_{1}$ and $+$ for $\psi_{2}$) and make additional imaginary time evolution. Then, the vortex separation gradually decreases from $2y_{0}$, but after some time, the decreasing rate reaches a quasi-stationary behavior characterized by the slow linear decrease of the energy, which has been also reported in Ref. Eto and Nitta (2018). We stop the imaginary time evolution at a certain time in this quasi-stationary regime and use this configuration with the vortex separation $L(<2y_{0})$ as the initial state of the simulation.

In this situation, we have to consider the finite size effect for the transverse instability. It is necessary for the unstable modes to grow, a half of their wave length should be smaller than the length $L$ of the domain wall. The condition is given by $\lambda/2=\pi/k_{\text{max}}<L$, where $k_{\text{max}}$ can be written as a linear function of $\omega_{\text{R}}$ [see Fig. 2(b)]. Using this condition, we get the boundary for the disintegration of the finite-size domain wall, as depicted in Fig. 5. With increasing the Rabi frequency $\omega_{\text{R}}$, the critical length $L$ of the wall decreases, which is consistent with the observation in Ref. Tylutki et al. (2016); Eto and Nitta (2018). The stability diagram is weakly dependent on the intercomponent coupling $\gamma$.

Starting from the initial state obtained by the above procedure, we monitor the real time dynamics of the domain wall with the length $L$ for given sets of $\gamma$ and $\omega_{\text{R}}$. The typical disintegration process is shown in Fig. 6 for $L=45$, $\gamma=0$, and $\omega_{\text{R}}=0.1$. The wall initially deforms into S-shape and breaks into small pieces at the curved points. After that, the remaining wall again deforms into inverted S-shape, repeating the similar breaking process. Eventually, the wall breaks into five pieces. In this simulation, we take $\gamma=0$ so that there is no repulsive vortex-vortex interaction. Thus, the fragmented walls eventually shrink as a result of the tension of the sine-Gordon domain wall.

When the domain wall disintegrates, the additional vortices are nucleated in each component through the snake instability. Thus, monitoring the number of vortices (phase defects) during the time developments provides a clear criterion for the disintegration; the number of vortices in each component is kept to be unity when the domain wall is stable, otherwise it is unstable. Using this criterion, we calculate the critical length of the wall as a function of $\omega_{\text{R}}$ for $\gamma=0,\pm 0.5$ and plot in Fig. 5. We see that the behavior is qualitatively consistent with the BdG prediction, but quantitatively the disintegration takes place even in the region of the Rabi frequencies smaller than the critical values.

This quantitative difference may attribute to the rotation of the domain wall with the finite size Tylutki et al. (2016); Eto and Nitta (2018). As mentioned before, the domain wall with the finite size exhibits a rotational motion due to the attractive tension by the Rabi coupling. If we consider the rotating frame co-moving with the domain wall, it is seen as a static configuration in this frame.

To see the stability of moving domain wall in a simple situation, we first consider the energetic stability of a moving domain wall with a constant velocity $V$ along the $x$-direction. In the co-moving frame of the velocity $V$, the GP equation reads

[TABLE]

Here, the velocity $V$ is scaled by the sound velocity $\xi/\tau$. Using Eq. (21), we analyze the energetic stability of the domain wall solution satisfying the boundary condition similar to Sec. III.1 through the imaginary time evolution. We show the stability boundary for $V=0.2$ and $V=0.4$ in Fig. 1; the stability region becomes narrower as $V$ increases.

Next, we turn to the rotating frame with the frequency $\Omega_{\text{prec}}$. When we suppose the domain wall along the $y$-axis and neglect the influence due to the vortices at the edge of the wall, we can approximately employ Eq. 21 with the $y$-dependent velocity $V=-\Omega_{\text{prec}}y$ to discuss the stability problem. The precession frequency of the vortex molecule is approximately given by $\Omega_{\text{prec}}\simeq 8\sqrt{2\omega_{\text{R}}}/(\pi L)$ Tylutki et al. (2016). For $\omega_{\text{R}}=0.1$ we find the velocity of the wall at the edge is $V_{\text{prec}}\sim 0.5$; the intermediate region of the wall precesses with the velocity below 0.5. Thus, when we take account of the rotational motion of the domain wall, the stability is significantly reduced for outer region away from the center of the wall. This behavior is clearly seen in Fig. 6, where the disintegration occurs first near the edges.

V Conclusion

We study the transverse instability and disintegration dynamics of the domain wall of the relative phase in Rabi-coupled two-component BECs, motivated by the numerical observation in Refs. Tylutki et al. (2016); Eto and Nitta (2018). Using the exact solutions, we construct the phase diagram representing the stability of the domain wall. In the unstable regime, the domain wall exhibits dynamical instability associated with the transverse displacement, known as the snake instability. The instability causes the disintegration of the domain wall into the small pieces, namely the vortex molecules. The growth time and the size of the resulting vortex molecules are in fairly agreement with the prediction by the BdG analysis.

Acknowledgements.

The work of K.K. is supported by KAKENHI from the Japan Society for the Promotion of Science (JSPS) Grant-in- Aid for Scientific Research (KAKENHI Grant No. 18K03472).