Interfaces between Bose–Einstein and Tonks–Girardeau atomic gases

Solutions of stationary equations (10) and (11) for the DW separating semi-infinite domains occupied by the BEC and TG components (or the spatial domains occupied by the copropagating waves in the above-mentioned optics model) are specified by the following boundary conditions (b.c.), which include the respective asymptotic densities, ${\phi }_{1}^{2}(x=-\infty )\equiv {({n}_{1})}_{\mathrm{asympt}}$ and ${\phi }_{2}^{2}(x\;=\;+\infty )\equiv {({n}_{2})}_{\mathrm{asympt}}$:

$$\begin{eqnarray}{({n}_{1})}_{\mathrm{asympt}} & = & {\mu }_{1}/g,\;{\phi }_{2}(x=-\infty )=0,\\ {\phi }_{1}(x\quad =\quad +\infty ) & = & 0,\;{({n}_{2})}_{\mathrm{asympt}}=\sqrt{{\mu }_{2}}.\end{eqnarray} \tag{ 14 }$$

The condition that formal Hamiltonian (12) must take the same values at $x\to -\infty$ and $x\to +\infty$ across the solution imposes a restriction on b.c. (14), which relates the two chemical potentials:

$$\begin{eqnarray}&&{\mu }_{1}^{2}=(4g/3){\mu }_{2}^{3/2}.\end{eqnarray} \tag{ 15 }$$

In terms of the asymptotic densities, condition (14) takes the form of

$$\begin{eqnarray}&&{({n}_{1})}_{\mathrm{asympt}}=(2/\sqrt{3g}){({n}_{2})}_{\mathrm{asympt}}^{3/2},\end{eqnarray} \tag{ 16 }$$

which actually implies the balance of the pressure applied to the DW from the two sides. Naturally, the mass ratio (m) does not appear in equations (15) and (16). Note that, if condition (16) between the densities does not hold initially, the DW in a finite system, which corresponds to real experimental settings, will move to a position at which the condition holds for the accordingly modified densities.

Further, the immiscibility condition for the BEC and TG, which is necessary for the existence of the DW separating the two quantum gases, is that the demixed configuration must provide a smaller energy density (defined as per equation (7)) than a uniformly mixed state with densities which are equal to half of asymptotic densities (14):

$$\begin{eqnarray}&&{({n}_{1})}_{\mathrm{uni}}=\displaystyle \frac{1}{2}{({n}_{1})}_{\mathrm{asympt}}\equiv \displaystyle \frac{{\mu }_{1}}{2g},\quad {({n}_{2})}_{\mathrm{uni}}=\displaystyle \frac{1}{2}{({n}_{2})}_{\mathrm{asympt}}\equiv \displaystyle \frac{\sqrt{{\mu }_{2}}}{2}.\end{eqnarray} \tag{ 17 }$$

The uniform state corresponds to values of the chemical potentials different from ${\mu }_{1}$ and ${\mu }_{2}$, namely, ${({\mu }_{1})}_{\mathrm{uni}}=g{({n}_{1})}_{\mathrm{uni}}+{({n}_{2})}_{\mathrm{uni}}$ , ${({\mu }_{2})}_{\mathrm{uni}}={({n}_{2})}_{\mathrm{uni}}^{2}+{({n}_{2})}_{\mathrm{uni}}$. Then, the comparison of the average energy densities of the demixed and uniform states gives rise to the immiscibility condition in the following form:

$$\begin{eqnarray}&&\displaystyle \frac{{\mu }_{1}^{2}}{g}+{\mu }_{2}^{3/2}\lt \displaystyle \frac{2{\mu }_{1}\sqrt{{\mu }_{2}}}{g}.\end{eqnarray} \tag{ 18 }$$

Further, the substitution of relation (15) in equation (18) transforms it into an inequality for the self-repulsion strength of the BEC component, g:

$$\begin{eqnarray}&&g\lt \displaystyle \frac{48}{49}\displaystyle \frac{1}{\sqrt{{\mu }_{2}}}=\displaystyle \frac{48}{49}\displaystyle \frac{1}{{({n}_{2})}_{\mathrm{asympt}}}\equiv {g}_{\mathrm{max}}.\end{eqnarray} \tag{ 19 }$$

In particular, for ${\mu }_{2}=1$ (as said above, this value will be fixed by rescaling), equation (19) yields ${g}_{\mathrm{max}}\quad =\quad 48/49\approx 0.9796$.

The DW may be characterized by an effective width of its core, which may be naturally defined by the following integral expression:

$$\begin{eqnarray}&&W=\displaystyle \frac{1}{{({n}_{1})}_{\mathrm{asympt}}{({n}_{2})}_{\mathrm{asympt}}}{\displaystyle \int }_{-\infty }^{+\infty }{\phi }_{1}^{2}(x){\phi }_{2}^{2}(x){\rm{d}}x.\end{eqnarray} \tag{ 20 }$$

The analysis should produce W as a function of parameters m and g, once the TG asymptotic density is fixed by setting ${\mu }_{2}=1$, see figure 3 below.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Close to the existence boundary of the DW, i.e., at $0\lt {g}_{\mathrm{max}}-g\ll {g}_{\mathrm{max}}$, the DW becomes very wide. In this limit, the dependence between the width and proximity to the threshold may be estimated as follows: the density of the gradient energy in equation (7) scales as $1/{W}^{2}$, hence the full gradient energy scales as $1/W$, and the respective effective force may be estimated as $-\partial H/\partial W$ $\sim 1/{W}^{2}$. It must be balanced by the effective bulk force vanishing at $g={g}_{\mathrm{max}}$, which scales as ${g}_{\mathrm{max}}-g$. Thus, the equilibrium condition predicts that the DW's width diverges at ${g}_{\mathrm{max}}-g\to 0$ as

$$\begin{eqnarray}&&W\sim 1/\sqrt{{g}_{\mathrm{max}}\quad -\quad g}.\end{eqnarray} \tag{ 21 }$$

Accurate analytical results for the DW can be obtained in the limit case corresponding to $g\gg 1$, or to $m\ll 1$, when the derivative term in equation (10) may be neglected (which actually implies the use of the TFA [11]), yielding

$$\begin{eqnarray}&&{\phi }_{1}^{2}=\displaystyle \frac{{\mu }_{1}-{\phi }_{2}^{2}}{g}.\end{eqnarray} \tag{ 22 }$$

Strictly speaking, the assumption of $m\ll 1$ may contradict the experimentally relevant way of the realization of the TG gas, based on making the respective effective mass large [19]. Nevertheless, the other option justifying the applicability of the TFA, $g\gg 1$, is quite relevant.

The substitution of expressions (22) and (15) in equation (11) leads to the single stationary equation for ${\phi }_{2}$, with the cubic–quintic nonlinearity:

$$\begin{eqnarray}&&{\mu }_{2}^{3/4}\left({\mu }_{2}^{1/4}-\displaystyle \frac{2}{\sqrt{3g}}\right){\phi }_{2}+\displaystyle \frac{1}{2m}\displaystyle \frac{{{\rm{d}}}^{2}{\phi }_{2}}{{\rm{d}}{x}^{2}}-{\phi }_{2}^{5}+\displaystyle \frac{1}{g}{\phi }_{2}^{3}=0.\end{eqnarray} \tag{ 23 }$$

(Equation (15) was used here to eliminate ${\mu }_{1}$.) Then, using a known particular exact solution for the DW solution of equation (23) [38], we obtain a solution at the following values of the chemical potentials:

$$\begin{eqnarray}&&{({\mu }_{2})}_{\mathrm{TFA}}={\left(\displaystyle \frac{3}{4g}\right)}^{2},\;{({\mu }_{1})}_{\mathrm{TFA}}=\displaystyle \frac{3}{4g},\end{eqnarray} \tag{ 24 }$$

i.e., ${({n}_{1})}_{\mathrm{asympt}}=3/(4{g}^{2})$, ${({n}_{2})}_{\mathrm{asympt}}=3/(4g)$, in the form of

$$\begin{eqnarray}&&{({\phi }_{2}(x))}_{\mathrm{TFA}}=\displaystyle \frac{1}{2}\sqrt{\displaystyle \frac{3}{g}\displaystyle \frac{1}{1+\mathrm{exp}(\pm \sqrt{3\;m/2}x/g)}}.\end{eqnarray} \tag{ 25 }$$

Note that ${({\mu }_{2})}_{\mathrm{DW}}$ satisfies condition (19), and rescaling (13) transforms values (24) into ${({\mu }_{2})}_{\mathrm{DW}}={({\mu }_{1})}_{\mathrm{DW}}=1$. The width of this solution, defined according to equation (20), is

$$\begin{eqnarray}&&{W}_{\mathrm{TFA}}=\sqrt{2/(3m)}g.\end{eqnarray} \tag{ 26 }$$

Comparison of analytical solution (25) with its numerical counterpart is displayed in figure 2.

An analytical approach can also be developed in the opposite limit of $m\to \infty$ (very heavy TG atoms), when the TFA may be applied to equation (11) (i.e., the derivative term may be dropped in this equation), reducing it to

$$\begin{eqnarray}&&{\phi }_{2}^{2}(x)=\sqrt{{\mu }_{2}-{\phi }_{1}^{2}(x)}.\end{eqnarray} \tag{ 27 }$$

In this case, equation (10) takes the form of

$$\begin{eqnarray}&&{\mu }_{1}{\phi }_{1}=-\displaystyle \frac{1}{2}\displaystyle \frac{{{\rm{d}}}^{2}{\phi }_{1}}{{\rm{d}}{x}^{2}}+g{\phi }_{1}^{3}+\sqrt{{\mu }_{2}-{\phi }_{1}^{2}}{\phi }_{1},\end{eqnarray} \tag{ 28 }$$

see equations (22) and (23) derived above in the opposite limit of $m\to 0$. In the particular case of g = 0 (no intrinsic interaction in the BEC component), equation (28) coincides with the equation considered in [39], in the context of a completely different physical model, for a resonantly absorbing Bragg reflector in optics. The formal Hamiltonian for equation (28) is

$$\begin{eqnarray}&&h=\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{\rm{d}}{\phi }_{1}}{{\rm{d}}x}\right)}^{2}+{\mu }_{1}{\phi }_{1}^{2}-\displaystyle \frac{g}{2}{\phi }_{1}^{4}-\displaystyle \frac{2}{3}{({\mu }_{2}-{\phi }_{1}^{2})}^{3/2},\end{eqnarray} \tag{ 29 }$$

see equation (12). The DW solution corresponds to a solution of equation (28) with the b.c. produced by equation (14): ${\phi }_{1}^{2}(x=-\infty )={\mu }_{1}/g,\quad {\phi }_{1}(x\quad =\quad +\infty )=0$, and, simultaneously, ${\phi }_{1}^{2}(x=-\infty )={\mu }_{2}$, the latter relation following from equation (27). Obviously, this b.c. set is self-consistent solely for ${\mu }_{1}=g{\mu }_{2}$. In the combination with this relation, equating values of h corresponding to $x=-\infty$ and $x\;=\;+\infty$, as per equation (29), yields $g\sqrt{{\mu }_{2}}=-4/3$, which is impossible, as g cannot be negative. Thus, the DW solution does not exist in this approximation. Nevertheless, it readily produces analytical solutions for complexes of other types, BD and DBS, as shown below.

Generic DW structures obtained from numerical solutions of equations (4) and (5) are displayed in figure 1. Stationary solutions have been obtained by means of the imaginary-time-propagation method. The validity and stability of the solutions was then checked by simulations in real time. The simulations in both imaginary and real time were run by means of the split-step Crank–Nicholson method [40].

Stable solutions have been found in a wide range of values of m at $g\lt 1$, and no solutions could be obtained at $g\gt 1$, which is readily explained by equation (19) (recall we set ${\mu }_{2}=1$). Actual numerical results for the DWs have been obtained for $g\leqslant 0.85$. In other words, a natural conclusion is that the DW exists if the BEC–TG repulsion is stronger than the intrinsic repulsion of the BEC component. In the remaining interval of $0.85\lt g\lt 48/49$ (see equation (19)), it is difficult to generate DWs numerically, as the initial guess for finding the solution by means of the imaginary-time propagation should be very close to the true one, in case the solution is sought for close to its existence boundary. Unlike g, the dependence of the DW solutions on m is very weak, as seen in figure 1 which displays examples for m = 0.1 and m = 10 (the comparison of the results obtained for large and small values of m, presented here and below, is appropriate even if very small values of m may be experimentally irrelevant, as mentioned above).

In addition, figure 2 displays the comparison of the analytical approximation (TFA) based on equations (22), (24), and (25) with the numerical solution obtained for the same values of parameters. It is observed that the analytically predicted density profiles are virtually identical to their numerical counterparts.

To characterize the entire family of the DW solutions, their width defined as per equation (20) has been numerically evaluated, as shown in figure 3(a), versus the mass ratio, m, in the interval of $0.05\leqslant m\leqslant 50$ . It is seen that the dependence of the width on m is rather weak, in agreement with examples displayed in figure 1. Further, the dependence of the DW's width on the BEC self-repulsion constant, g, is shown in figure 3(b). The range of values of g displayed in the figure is bounded by existence limits of the DW solutions, as given by equation (19). The increase of the width with g is a natural property of the system, as stronger self-repulsion stretches the transient layer (in the analytical form, this is clearly shown by equation (21) and solution (25).)

BD solutions are even ones, with ${\phi }_{\mathrm{1,2}}(-x)={\phi }_{\mathrm{1,2}}(x)$, singled out by b.c.

$$\begin{eqnarray}{\left.\displaystyle \frac{{\rm{d}}{\phi }_{1}}{{\rm{d}}x}\right|}_{x=0} & = & {\left.\displaystyle \frac{{\rm{d}}{\phi }_{2}}{{\rm{d}}x}\right|}_{x=0}=0\\ {\phi }_{1}(x\quad =\quad +\infty ) & = & 0,\;{\phi }_{2}^{2}(x\quad =\quad +\infty )\equiv {({n}_{2})}_{\mathrm{asympt}}=\sqrt{{\mu }_{2}},\end{eqnarray} \tag{ 30 }$$

(a BEC drop embedded into the TG background), or

$$\begin{eqnarray}{\left.\displaystyle \frac{{\rm{d}}{\phi }_{1}}{{\rm{d}}x}\right|}_{x=0} & = & {\left.\displaystyle \frac{{\rm{d}}{\phi }_{2}}{{\rm{d}}x}\right|}_{x=0}=0\\ {\phi }_{1}^{2}(x\quad =\quad +\infty ) & \equiv & {({n}_{1})}_{\mathrm{asympt}}={\mu }_{1}/g,\;{\phi }_{2}(x\quad =\quad +\infty )=0,\end{eqnarray} \tag{ 31 }$$

(a TG drop embedded into the BEC background), see equation (14). In fact, the BDs may be considered as bound pairs of the DWs, with the layer BEC or TG trapped between the two semi-infinite TG or BEC domains, respectively.

The existence and stability of these solutions has been numerically checked by fixing the chemical potential of the background to unity, and varying the norm of the drop, or more precisely:

• setting ${\mu }_{2}=1$ for the TG background, see equation (30), and selecting several fixed values of N₁ for the BEC drop, see equation (8);
• setting ${\mu }_{1}=1$ for the BEC background, see equation (31), and selecting several fixed values of N₂ for the TG drop, see equation (8).

It is relevant to stress that, unlike the DWs, the chemical potentials of the bubble and drop components are not related by any condition similar to equation (15). In this setting, the remaining free parameters are relative mass m and BEC self-repulsion coefficient g.

Analytical results for the BD structures are available in the limit case of $m\to \infty$ (heavy TG atoms), which corresponds to equations (27) and (28). Indeed, in this approximation the BD solutions, defined by b.c. (30) or (31) correspond, respectively, to bright solitons or bubbles produced by equation (28) ('bubbles' are solutions of NLSEs in the form of a local drop of the density supported by the flat background, without zero crossing, and without a phase shift at $x\to \pm \infty$, unlike dark solitons [34]). Further, the analysis of the corresponding Hamiltonian (29) demonstrates that equation (28) cannot generate bright solitons, but it readily gives rise to bubbles, i.e., the BD patterns in the form of the TG drop embedded into the BEC background. The fact that the BDs exist in the form of the TG drop embedded into the BEC bubble, but do not exist in the reverse form, is confirmed by numerical results presented below.

Analytical results for the BD structures are also available in the limit case of $g\to \infty$ or $m\to 0$, when the derivative term may be neglected in equation (10), leading to equation (22), as shown above. The substitution of this approximation into equation (11) leads to the single equation with the cubic–quintic nonlinearity. For ${\mu }_{1}$ not linked to ${\mu }_{2}$ by relation (15), which was specific for the DW, but is not relevant in the present case, the cubic–quintic equation takes a form slightly more general than equation (23) derived above:

$$\begin{eqnarray}&&\left({\mu }_{2}-\displaystyle \frac{{\mu }_{1}}{g}\right){\phi }_{2}=-\displaystyle \frac{1}{2m}\displaystyle \frac{{{\rm{d}}}^{2}{\phi }_{2}}{{\rm{d}}{x}^{2}}+\left({\phi }_{2}^{4}-\displaystyle \frac{1}{g}{\phi }_{2}^{2}\right){\phi }_{2}.\ \end{eqnarray} \tag{ 32 }$$

A straightforward consideration demonstrates that equation (32) gives rise to bubbles, i.e., effectively, the BDs of type (30), in the case of

$$\begin{eqnarray}&&\displaystyle \frac{3}{8}\lt m\left(\displaystyle \frac{{\mu }_{1}}{g}-{\mu }_{2}\right)\lt \displaystyle \frac{1}{2},\end{eqnarray} \tag{ 33 }$$

and to bright solitons, i.e., the BDs of type (31), in the range of

$$\begin{eqnarray}&&0\lt m\left(\displaystyle \frac{{\mu }_{1}}{g}-{\mu }_{2}\right)\lt \displaystyle \frac{3}{8}.\end{eqnarray} \tag{ 34 }$$

In fact, the present limit of $g\to \infty$ or $m\to 0$ implies that interval (33) does not exist, while condition (34) may hold. Thus, this consideration again predicts that the BD exist solely in the form of a TG drop (bright soliton) embedded into the BEC bubble, which is exactly confirmed by numerical results following below.

In terms of the above-mentioned optics model, this conclusion means that a bright soliton driven by the quintic self-defocusing may be embedded into the background filled by the optical wave subject to the cubic self-defocusing, but not vice versa.

For the state defined as per equation (30) (a BEC/cubic drop embedded into the TG/quintic background), it was not possible to find solutions in the entire parametric region investigated, which amounts to intervals $0.2\leqslant g\leqslant 5$ and $0.1\leqslant m\leqslant 10$, and $0.1\leqslant {N}_{1}\leqslant 10$. In imaginary-time simulations, this configuration, if taken as an input, transforms into a flat one. The nonexistence of this BD species is readily explained by the analytical results presented above.

On the other hand, also in agreement with the above analysis, in the same parameter region (with N₁ replaced by N₂) stable solutions have been readily found for the BD state defined by equation (31) (a TG/quintic drop embedded into the BEC/cubic background), see a typical example in figure 4. The BD may be characterized by the missing mass of the bubble void, i.e.

$$\begin{eqnarray}&&{M}_{\mathrm{void}}={[{({n}_{1})}_{\mathrm{asympt}}]}^{-1}{\displaystyle \int }_{-\infty }^{+\infty }[{({n}_{1})}_{\mathrm{asympt}}-{\phi }_{1}^{2}(x)]{\rm{d}}x.\end{eqnarray} \tag{ 35 }$$

Computations demonstrate that, for the family of the BD states, ${M}_{\mathrm{void}}$ very weakly depends on relative mass m. Figure 5 displays the missing mass of the BEC component as a function of g. The dependence observed in figure 5 may be approximated by relation ${M}_{\mathrm{void}}\simeq {N}_{2}/g$, which is easy to understand: according to equation (10), the TG component with density ${\phi }_{2}^{2}$ replaces the missing BEC field with effective density ${({\phi }_{1}^{2})}_{\mathrm{void}}={\phi }_{2}^{2}/g$. This argument also explains an approximately linear increase of ${M}_{\mathrm{void}}$ with the growth of the norm of the TG/quintic drop, N₂, as clearly seen in figure 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Bright–dark soliton (DBS) solutions represent complexes with one even localized (bright) component, and the other delocalized spatially odd zero-crossing one (the dark soliton), the difference from the BD structure being that in the latter case both components were patterned as even ones, without zero crossing in the bubble structure. The corresponding b.c. for the DBS complex are given by

$$\begin{eqnarray}{\left.\displaystyle \frac{{\rm{d}}{\phi }_{1}}{{\rm{d}}x}\right|}_{x=0} & = & {\phi }_{2}(x=0)=0\\ {\phi }_{1}(x\quad =\quad +\infty )=0,\;{\phi }_{2}^{2}(x\quad =\quad +\infty ) & \equiv & {({n}_{2})}_{\mathrm{asympt}}=\sqrt{{\mu }_{2}},\end{eqnarray} \tag{ 36 }$$

(a bright BEC component embedded into the TG background), or

$$\begin{eqnarray}{\phi }_{1}(x=0) & = & {\left.\displaystyle \frac{{\rm{d}}{\phi }_{2}}{{\rm{d}}x}\right|}_{x=0}=0\\ {\phi }_{1}^{2}(x\quad =\quad +\infty ) & \equiv & {({n}_{1})}_{\mathrm{asympt}}={\mu }_{1}/g,\;{n}_{2}(x\quad =\quad +\infty )=0,\end{eqnarray} \tag{ 37 }$$

(a bright TG component embedded into the BEC background), see equations (30) and (31). Note that, similar to the case of the BD, the DBS existence condition does not impose any relation on the chemical potentials of its components, unlike equation (15) for the DW.

The limit case of $m\to \infty$ (heavy TG atoms), which is represented by equations (27) and (28), makes it possible to produce analytical results for the DBS defined by b.c. (37), which corresponds to a dark soliton of equation (28). The analysis of the Hamiltonian (29) demonstrates that equation (28) indeed gives rise to dark solitons, i.e., the DBS patterns in the form of the bright TG soliton embedded into the dark BEC soliton.

The consideration of the DBS is also possible in the opposite limit of $g\to \infty$ or $m\to 0$, which amounts to the consideration of equation (32). Dark solitons of that equation correspond to the DBS defined by b.c. (36), and it is easy to see that equation (32) gives rise to the dark solitons in the range of $-\infty \lt m({\mu }_{1}/g-{\mu }_{2})\lt 3/8$. Obviously, this condition holds in the present limit of $g\to \infty$ or $m\to 0$.

Thus, the analysis of the limit cases readily predicts the existence of both species of the DBSs, which correspond to b.c. (36) and (37). This prediction is confirmed by the following numerical results.

Both types of the DBS have been produced by the numerical computations, as shown in figure 6.

Zoom In Zoom Out Reset image size

The DBS complexes have their existence boundaries. They cease to exist when the self-repulsion BEC coefficient, g, becomes too large. This is shown in figure 7, that displays the boundaries in the plane of $(g,m)$. It is seen that a lower chemical potential of the dark component favors the existence of the DBS, as the boundaries move towards higher values of g with the decrease of the chemical potential. Relative mass m plays a different role: the existence area for the DBS complex with the bright BEC (cubic) component shrinks with the increase of m, while for the case of the bright TG (quintic) component the increase of m favors the DBS existence.

Zoom In Zoom Out Reset image size