Interfaces between Bose-Einstein and Tonks-Girardeau atomic gases

We consider one-dimensional mixtures of an atomic Bose-Einstein condensate (BEC) and Tonks- Giradeau (TG) gas. The mixture is modeled by a coupled system of the Gross-Pitaevskii equation for the BEC and the quintic nonlinear Schroedinger equation for the TG component. An immiscibility condition for the binary system is derived in a general form. Under this condition, three types of BEC-TG interfaces are considered: domain walls (DWs) separating the two components; bubble-drops (BDs), in the form of a drop of one component immersed into the other (BDs may be considered as bound states of two DWs); and bound states of bright and dark solitons (BDSs). The same model applies to the copropagation of two optical waves in a colloidal medium. The results are obtained by means of systematic numerical analysis, in combination with analytical Thomas-Fermi approximations (TFAs). Using both methods, families of DW states are produced in a generic form. BD complexes exist solely in the form of a TG drop embedded into the BEC background. On the contrary, BDSs exist as bound states of TG bright and BEC dark components, and vice versa.


INTRODUCTION
Binary systems, whose behavior crucially depends on the underlying condition of immiscibility or miscibility 1 , play a fundamentally important role in many areas of physics. In the case of immiscibility, a major effect is the formation of domain walls (DWs) between regions occupied by immiscible components. Commonly known are DWs in media featuring a vectorial order parameter, such as ferromagnets 2 , ferroelectrics 3 , and liquid crystals 4 . In self-defocusing optical media, DWs separate regions occupied by electromagnetic waves with orthogonal circular polarizations of light 5,6 . Similar interface patterns were predicted in arrays of nonlinear optical waveguides, modeled by discrete nonlinear Schrödinger equations (NLSEs) 8 .
DWs are known in superfluids too, where they are formed by immiscible binary Bose-Einstein condensates (BECs), as predicted theoretically 9 and demonstrated in experiments 10 . In the mean-field approximation 11 , such settings are modeled by systems of nonlinearly coupled Gross-Pitaevskii equations (GPEs) with the cubic self-repulsive nonlinearity, which are similar to coupled NLSEs describing the above-mentioned optical DWs 5, 6 . In their stationary form, these equations coincide with coupled cubic Ginzburg-Landau equations modeling DWs in dissipative patterns, such as interfaces between rolls with different orientations in large-area Rayleigh-Benard convection 7 .
The analysis of the DWs in BEC was extended for broader settings, including linear interconversion between the immiscible components (this is possible when they represent two different hyperfine states of the same atom coupled by a resonant radiofrequency wave) 12 , dipolar 13 and spinor (three-component) condensates 14 , as well as the BEC discretized by trapping in a deep optical-lattice (OL) potentials 15 . Furthermore, the study of the DWs was recently extended for immiscible binary BECs with three-particle collisions 16 , in the case when the related losses may be neglected, the respective coupled GPEs featuring the cubic-quintic repulsive nonlinearity 17 .
In the effectively one-dimensional (1D) setting, ultracold bosonic gases with strong inter-atomic repulsion may be cast in the Tonks-Girardeau (TG) state, which emulates the gas of non-interacting fermions 18 , provided that the energy of the repulsive interaction between bosons exceeds their kinetic energy, while the opposite situation corresponds to the BEC phase in the bosonic gas (a review of the TG model was given in Ref. 21 ). The TG gas of hard-core bosons has been realized experimentally, using tight transverse confinement 19,20 . In particular, a longitudinal OL potential was used to increase the effective mass in the trapped state, thus making the kinetic energy small enough 19 .
It is commonly known that GPEs furnish very accurate description of the BEC in atomic gases. A similar macroscopic model of the TG gas is offered by the NLSE with the quintic self-repulsion term 22 . In a rigorous form, the relevance of the corresponding sextic term in the free-energy density of the three-dimensional bosonic gas in its ground state, which reduces to the quasi-1D TG phase, was demonstrated in Ref. 23 , under condition GL ≫ N , where G, L, and N are, respectively, the inter-atomic repulsion strength, system's length, and the total number o atoms. The quintic model was used in various contexts, including shock waves 24 , dark 25 and gap-mode 26 solitons, as well as bright solitons supported by dipole-dipole interactions 27 , and, recently, DWs in immiscible binary TG gases 16 . Further, oscillation frequencies derived from fermionic hydrodynamic equations, which apply to the hard-core TG gas, were found to be close to their counterparts predicted by the quintic NLSE 28 . Coupled quintic NLSEs also arise in works aimed at constructing the ground state of a binary TG mixture in the harmonic-oscillator potential by means of the density-functional method 29 . On the other hand, this approach may not apply to TG gases beyond the framework of static configurations and hydrodynamic regimes. In particular, it fails for strongly non-equilibrium problems, such as merger of distinct gas clouds 30 .
As concerns the mixtures, it may be interesting to consider binary systems including the TG gas and another quantum-gas component. In particular, exact solutions were found for the ground state of TG-Fermi mixtures 31 . The binary gas of impenetrable bosons is solvable too 32 . The objective of the present work is to introduce basic nonlinear complexes, such as DWs, bubble-drop (BD) modes (bound pairs of two DWs), and dark-bright solitons (DBSs), in an immiscible system of TG and BEC gases. In the experiment, the system may be realized, in particular, as a bosonic gas composed of two atomic species under tight transverse confinement, with a longitudinal OL potential acing on (being relatively close to a resonance with) one species only. Then, as the experimental setting presented in Ref. 19 suggests, a large effective mass of the near-resonant component will bring it into the TG state, while the other component may stay in the BEC phase.
As a model for this system, in Section II we adopt the cubic GPE for the self-repulsive BEC component coupled by the cubic (collisional) repulsive term to the quintic NLSE for the TG species. The use of the latter equation is appropriate, as we study only static configurations of the system. The same model may find a realization in optics as a model of colloidal waveguides. It has been recently demonstrated that, selecting the size of metallic nanoparticles in the colloid and their concentration, one can engineer desirable coefficients of the corresponding cubic and quintic nonlinearity 33 . In particular, it is possible to design a waveguide which features a nearly pure quintic nonlinearity at a particular wavelength, while the cubic response dominates at a different wavelength. The copropagation of optical signals carried by these wavelengths will then emulate the TG-BEC system.
DWs separating the BEC and TG phases are addressed in Section III. We derive the respective immiscibility condition (see Eq. (19) below), and then generate DW states in a systematic way, using both numerical solutions and the analytical Thomas-Fermi approximation (TFA). The latter method makes it possible to obtain some DWs in an explicit analytical form, as given below by Eqs. (24) - (26). DB and DBS complexes are considered in Sections IV and V, respectively, again using a combination of analytical and numerical methods. The analysis predicts that the DB states exist solely in the form of the TG drop embedded into the BEC background (bubble), but not in the opposite case; on the other hand, the DBS are predicted in either case of the bright TG soliton embedded into the BEC dark soliton, or vice versa. These predictions are fully corroborated by numerical results. The paper is concluded by Section VI.

II. THE MODEL
The system of coupled NLSEs with the cubic and quintic nonlinear terms for the BEC and TG components, ψ 1 and ψ 2 , with respective scaled masses m 1 and (1) (i.e., m is the relative effective mass of TG components which, as said above, may be made larger than the actual atomic mass 19 ) is where real parameters γ ≥ 0 and Γ > 0 are strengths of the self-repulsion of the BEC component, and repulsion between the BEC and TG ones, respectively, while and the coefficient of the effective quintic self-repulsion of the TG component are scaled to be 1 (in the notation of Ref. 22 , the natural value of the latter coefficient is π 2 ). Substitution makes it possible to further fix Γ ≡ 1 and m 1 ≡ 1, thus simplifying Eq. (2) to a system with two free coefficients, m and g: For the above-mentioned spatial-domain optical model (with t replaced by the propagation distance, Z, and transverse coordinate X), the scaled propagation equations are derived, using the standard procedure 37 , as Here Ψ 1 and Ψ 2 represent amplitudes of the copropagating electromagnetic waves with relative wavenumber and frequency k = k 2 /k 1 and ω = ω 2 /ω 1 , cf. Eq. (1). Further, G 3,5 in Eq. (6) are the cubic and quintic SPM coefficients for the two waves, while the cubic XPM coefficient is normalized to be 1. Additional rescaling, transforms Eq. (6) into the system of equations (4) and (5), with m ≡ k.
The Hamiltonian corresponding to Eqs. (4) and (5) is In addition to H, the system preserves the norms (scaled numbers of atoms in the ultracold gas, or total powers of the two waves, in terms of the optical model), of the two components, and, for dynamical solutions, also the total momentum, P = i +∞ −∞ 2 n=1 ψ n (∂ψ * n /∂x) dx, although, as mentioned above, the use of the quintic NLS equation for the description of dynamics of the TG gas may be impugnable.
Normalizing the stationary wave functions as one may fix µ 2 = 1, which implies that the density of the uniform TG component is also fixed to be 1, see Eq. (14) below. Numerical results for DWs are presented in the following section chiefly for this case, while in Sections IV and V other normalizations are used for BD and DBS states.

A. Analytical considerations
Solutions of stationary equations (10) and (11) for the domain wall (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, φ 2 1 (x = −∞) ≡ (n 1 ) asympt and φ 2 2 (x = +∞) ≡ (n 2 ) asympt : The condition that formal Hamiltonian (12) must take the same values at x → −∞ and x → +∞ across the solution imposes a restriction on b.c. (14), which relates the two chemical potentials: In terms of the asymptotic densities, condition (14) takes the form of 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 Eqs. (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 Eq. (7)) than a uniformly mixed state with densities which are equal to half of asymptotic densities (14): The uniform state corresponds to values of the chemical potentials different from µ 1 and µ 2 , namely, (µ 1 ) uni = g (n 1 ) uni + (n 2 ) uni , (µ 2 ) uni = (n 2 ) 2 uni + (n 2 ) 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: Further, the substitution of relation (15) in Eq. (18) transforms it into an inequality for the self-repulsion strength of the BEC component, g: In particular, for µ 2 = 1 (as said above, this value will be fixed by rescaling), Eq. (19) yields g max = 48/49 ≈ 0.9796. The DW may be characterized by an effective width of its core, which may be naturally defined by the following integral expression: The analysis should produce W as a function of parameters m and g, once the TG asymptotic density is fixed by setting µ 2 = 1, see Fig. 3 below. Close to the existence boundary of the DW, i.e., at 0 < g max − g ≪ g 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 Eq. (7) scales as 1/W 2 , hence the full gradient energy scales as 1/W , and the respective effective force may be estimated as −∂H/∂W ∼ 1/W 2 . It must be balanced by the effective bulk force vanishing at g = g max , which scales as g max − g. Thus, the equilibrium condition predicts that the DW's width diverges at g max − g → 0 as Accurate analytical results for the DW can be obtained in the limit case corresponding to g ≫ 1, or to m ≪ 1, when the derivative term in Eq. (10) may be neglected (which actually implies the use of the TFA 11 ), yielding Strictly speaking, the assumption of m ≪ 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 ≫ 1, is quite relevant. The substitution of expressions (22) and (15) in Eq. (11) leads to the single stationary equation for φ 2 , with the cubic-quintic nonlinearity: (Eq. (15) was used here to eliminate µ 1 ). Then, using a known particular exact solution for the DW solution of Eq. (23) 38 , we obtain a solution at the following values of the chemical potentials: i.e., (n 1 ) asympt = 3/ 4g 2 , (n 2 ) asympt = 3/ (4g), in the form of (φ 2 (x)) TFA = 1 2 Note that (µ 2 ) DW satisfies condition (19), and rescaling (13) Comparison of analytical solution (25) with its numerical counterpart is displayed below in Fig. 2.
An analytical approach can also be developed in the opposite limit of m → ∞ (very heavy TG atoms), when the TFA may be applied to Eq. (11) (i.e., the derivative term may be dropped in this equation), reducing it to In this case, Eq. (10) takes the form of cf. Eqs. (22) and (23) derived above in the opposite limit of m → 0. In the particular case of g = 0 (no intrinsic interaction in the BEC component), Eq. (28) coincides with the equation considered in Ref. 39 , in the context of a completely different physical model, for a resonantly absorbing Bragg reflector in optics. The formal Hamiltonian for Eq. (28) is cf. Eq. (12). The DW solution corresponds to a solution of Eq. (28) with the b.c. produced by Eq. (14): φ 2 1 (x = −∞) = µ 1 /g, φ 1 (x = +∞) = 0, and, simultaneously, φ 2 1 (x = −∞) = µ 2 , the latter relation following from Eq. (27. Obviously, this b.c. set is self-consistent solely for µ 1 = gµ 2 . In the combination with this relation, equating values of h corresponding to x = −∞ and x = +∞, as per Eq. (29), yields g √ µ 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.

B. Numerical findings
Generic DW structures obtained from numerical solutions of Eqs. (4) and (5) are displayed in Fig. 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 < 1, and no solutions could be obtained at g > 1, which is readily explained by Eq. (19) (recall we set µ 2 = 1). Actual numerical results for the DWs have been obtained for g ≤ 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 < g < 48/49 (see Eq. (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 Fig. 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 terms of the optics model, the BEC and TG densities correspond to power densities of the waves subject to the action of the cubic and quintic self-defocusing nonlinearity, respectively (in captions to the following figures, they are referred to as "cubic" and "quintic" components, respectively).
In addition, Fig. 2 displays the comparison of the analytical approximation (TFA) based on Eqs. (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 Eq. (20) has been numerically evaluated, as shown in Fig. 3(a), vs. the mass ratio, m, in the interval of 0.05 ≤ m ≤ 50. It is seen that the dependence of the width on m is rather weak, in agreement with examples displayed in Fig. 1. Further, the dependence of the DW's width on the BEC self-repulsion constant, g, is shown in Fig. 3(b). The range of values of g displayed in the figure is bounded by existence limits of the DW solutions, as given by Eq. (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 Eq. (21) and solution (25).]
(a BEC drop embedded into the TG background), or (a TG drop embedded into the BEC background), cf. Eq. (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 µ 2 = 1 for the TG background, see Eq. (30), and selecting several fixed values of N 1 for the BEC drop, see Eq. (8); • setting µ 1 = 1 for the BEC background, see Eq. (31), and selecting several fixed values of N 2 for the TG drop, see Eq. (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 Eq. (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 → ∞ (heavy TG atoms), which corresponds to Eqs. (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 Eq. (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 → ±∞, unlike dark solitons 34 ). Further, the analysis of the corresponding Hamiltonian (29) demonstrates that Eq. (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 → ∞ or m → 0, when the derivative term may be neglected in Eq. (10), leading to Eq. (22), as shown above. The substitution of this approximation into Eq. (11) leads to the single equation with the cubic-quintic nonlinearity. For µ 1 not linked to µ 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 Eq. (23) derived above: A straightforward consideration demonstrates that Eq. (32) gives rise to bubbles, i.e., effectively, the BDs of type (30), in the case of and to bright solitons, i.e., the BDs of type (31), in the range of In fact, the present limit of g → ∞ or m → 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.

B. Numerical findings
For the state defined as per Eq. (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 ≤ g ≤ 5 and 0.1 ≤ m ≤ 10, and 0.1 ≤ N 1 ≤ 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 1 replaced by N 2 ) stable solutions have been readily found for the BD state defined by Eq. (31) (a TG/quintic drop embedded into the BEC/cubic background), see a typical example in Fig. 4. The BD may be characterized by the missing mass of the bubble void, i.e., Computations demonstrate that, for the family of the BD states, M 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 Fig. 5 may be approximated by relation M void ≃ N 2 /g, which is easy to understand: according to Eq. (10), the TG component with density φ 2 2 replaces the missing BEC field with effective density φ 2 1 void = φ 2 2 /g. This argument also explains an approximately linear increase of M void with the growth of the norm of the TG/quintic drop, N 2 , as clearly seen in Fig. 5. Missing mass 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 (a bright BEC component embedded into the TG background), or (a bright TG component embedded into the BEC background), cf. Eqs. (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 Eq. (15) for the DW.
The limit case of m → ∞ (heavy TG atoms), which is represented by Eqs. (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 Eq. (28). The analysis of the Hamiltonian (29) demonstrates that Eq. (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 → ∞ or m → 0, which amounts to the consideration of Eq. (32). Dark solitons of that equation correspond to the DBS defined by b.c. (36), and it is easy to see that Eq. (32) gives rise to the dark solitons in the range of −∞ < m (µ 1 /g − µ 2 ) < 3/8. Obviously, this condition holds in the present limit of g → ∞ or m → 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.

B. Numerical results
Both types of the DBS have been produced by the numerical computations, as shown in Fig. 6. 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 Fig. 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.

VI. CONCLUSION
We have introduced a binary system in the form of an immiscible pair of BEC and TG quantum gases, in the effectively 1D setting. Using the system of the GPE for the BEC, nonlinearly coupled to the quintic NLSE for the TG component, we have analyzed the possibility of the existence of three types of interfaces in this binary system: DWs (domain walls), BDs (bubble-drops), and BDSs (bright-dark solitons). The same model applies to a bimodal light propagation in a specially designed colloidal medium in optics. The immiscibility condition for the binary system of the present type was obtained in the general form. Analytical results were produced by means of the TFAs (Thomas-Fermi approximations), in the combination with the systematic numerical analysis. The DWs have been found an explicit analytical form, by means of the TFA, and in the generic form numerically. The analysis and numerical results demonstrate that the BDs exist solely in the form of the TG drop immersed into the BEC background, while the BDS complexes may be built equally well of the TG bright and BEC dark components or vice versa. Thus, the predicted results suggest new experiments with tightly confined two-species mixtures of ultracold bosonic gases. The analysis can be extended for a chain of BDs and/or DBSs in a long system, including a circular one, which corresponds to the binary gas loaded into a tight toroidal trap 35 .