Expansion Dynamics of a Two-Component Quasi-One-Dimensional Bose-Einstein Condensate: Phase Diagram, Self-Similar Solutions, and Dispersive Shock Waves

TL;DR

This paper studies the expansion behavior of a two-component Bose-Einstein condensate in a quasi-one-dimensional trap, classifies initial states, derives self-similar solutions, and explores dispersive shock waves, providing both analytical and numerical insights.

Contribution

It introduces a phase diagram for initial states, derives self-similar solutions, and analyzes dispersive shock waves in two-component condensates, advancing understanding of their expansion dynamics.

Findings

Classified initial states with a phase diagram based on interaction constants.

Derived self-similar solutions for specific cases of condensate expansion.

Numerically validated analytical solutions and explored dispersive shock wave formation.

Abstract

We investigate the expansion dynamics of a Bose-Einstein condensate that consists of two components and is initially confined in a quasi-one-dimensional trap. We classify the possible initial states of the two-component condensate by taking into account the non-uniformity of the distributions of its components and construct the corresponding phase diagram in the plane of nonlinear interaction constants. The differential equations that describe the condensate evolution are derived by assuming that the condensate density and velocity depend on the spatial coordinate quadratically and linearly, respectively, what reproduces the initial equilibrium distribution of the condensate in the trap in the Thomas-Fermi approximation. We obtained self-similar solutions of these differential equations for several important special cases and wrote out asymptotic formulas describing the condensate motion on long time scales, when the condensate density becomes so low that the interaction between atoms can be neglected. The problem on the dynamics of immiscible components with the formation of dispersive shock waves was also considered. We compare the numerical solutions of the Gross-Pitaevskii equations with their approximate analytical solutions and study numerically the situations when the analytical method admits no exact solutions.

Figures

Equations

$$\begin{array}{r}i\hslash \frac{\partial {\psi }_i}{\partial t}=-\frac{{\hslash }^2}{2m_i}\Delta {\psi }_i+g_{ii}|{\psi }_i{|}^2{\psi }_i+g_{ij}|{\psi }_j{|}^2{\psi }_i+U_i{\psi }_i,\end{array}$$

$$g_{ii}=\frac{4\pi {\hslash }^2{a}_{ii}}{m_i},g_{ij}=\frac{2\pi {\hslash }^2{a}_{ij}}{m_{ij}},$$

$$\int |{\psi }_i{|}^{2}dV=N_i,$$

$${\text{u}}_i=\frac{\hslash }{m_i}\nabla {\phi }_i.$$

$$g_{12}^2<g_{11}{g}_{22}.$$

$${\psi }_i=\sqrt{{\rho }_i(\text{r},t)}\exp \left(i\frac{m_i}{\hslash }{\int }^{\text{r}}{\text{u}}_i({\text{r}}^{\prime },t)d{\text{r}}^{\prime }-i\frac{{\mu }_i}{\hslash }t\right),$$

$$\frac{\partial {\rho }_i}{\partial t}+\nabla ({\rho }_i{\mathbf{u}}_i)=0,$$

$$\begin{array}{rl}\frac{\partial {\mathbf{u}}_i}{\partial t}& +({\mathbf{u}}_i\nabla ){\mathbf{u}}_i+\frac{g_{ii}}{m_i}\nabla {\rho }_i+\frac{g_{ij}}{m_i}\nabla {\rho }_j+\frac{\nabla U_i}{m_i}+\frac{{\hslash }^2}{4m_i^2}\nabla \left(\frac{(\nabla {\rho }_i{)}^2}{2{\rho }_i^2}-\frac{\Delta {\rho }_i}{{\rho }_i}\right)=0.\end{array}$$

$$\xi \sim \frac{\hslash }{\sqrt{g\rho m}}.$$

$$\frac{\partial {\rho }_i}{\partial t}+\nabla ({\rho }_i{\mathbf{u}}_i)=0,$$

$$\frac{\partial {\mathbf{u}}_i}{\partial t}+({\mathbf{u}}_i\nabla ){\mathbf{u}}_i+\frac{g_{ii}}{m_i}\nabla {\rho }_i+\frac{g_{ij}}{m_i}\nabla {\rho }_j+\frac{\nabla U_i}{m_i}=0.$$

$$U_i=\frac{1}{2}{m}_i[{\omega }_x^2{x}^2+{\omega }_{\perp }^2(y^2+z^2)],$$

$$\begin{array}{r}i\hslash \frac{\partial {\psi }_i}{\partial t}=-\frac{{\hslash }^2}{2m_i}\frac{{\partial }^2{\psi }_i}{\partial x^2}+g_{ii}|{\psi }_i{|}^2{\psi }_i+g_{ij}|{\psi }_j{|}^2{\psi }_i+\frac{1}{2}{m}_i{\omega }_i^2{x}^2{\psi }_i,\end{array}$$

$$\frac{\partial {\rho }_i}{\partial t}+\frac{\partial }{\partial x}({\rho }_i{u}_i)=0,$$

$$\frac{\partial u_i}{\partial t}+u_i\frac{\partial u_i}{\partial x}+\frac{g_{ii}}{m_i}\frac{\partial {\rho }_i}{\partial x}+\frac{g_{ij}}{m_i}\frac{\partial {\rho }_j}{\partial x}+{\omega }_i^{2}x=0.$$

$$\begin{split}H=&\int\Big{[}\frac{\hbar^{2}}{2m_{1}}|\nabla\psi_{1}|^{2}+\frac{\hbar^{2}}{2m_{2}}|\nabla\psi_{2}|^{2}+\frac{1}{2}(g_{11}|\psi_{1}|^{4}+2g_{12}|\psi_{1}|^{2}|\psi_{2}|^{2}+g_{22}|\psi_{2}|^{4})+\\ &+\frac{1}{2}m_{1}\omega_{1}^{2}x^{2}|\psi_{1}|^{2}+\frac{1}{2}m_{2}\omega_{2}^{2}x^{2}|\psi_{2}|^{2}\Big{]}dx.\end{split}$$

$$\begin{array}{r}H=\frac{1}{2}\int (g_{11}{\rho }_1^2+2g_{12}{\rho }_1{\rho }_2+g_{22}{\rho }_2^2+m_1{\omega }_1^2{x}^2{\rho }_1+m_2{\omega }_2^2{x}^2{\rho }_2)dx.\end{array}$$

$${\rho }_i^o(x)=\frac{2(g_{jj}{\mu }_i-g_{ij}{\mu }_j)-(g_{jj}{m}_i{\omega }_i^2-g_{ij}{m}_j{\omega }_j^2)x^2}{2(g_{ii}{g}_{jj}-g_{ij}^2)},$$

$${\rho }_i^s(x)=\frac{2{\mu }_i-m_i{\omega }_i^2{x}^2}{2g_{ii}},$$

$$\int {\rho }_{1}dx=N_1,\int {\rho }_{2}dx=N_2,$$

$$\frac{g_{12}}{g_{22}}=\frac{\frac{N_1}{N_2}\frac{g_{11}}{g_{12}}-\frac{1}{2}+\sqrt{(1+\frac{N_1}{N_2}\frac{m_1{\omega }_1^2}{m_2{\omega }_2^2})\frac{N_1}{N_2}\frac{g_{11}}{g_{12}}+\frac{1}{4}}}{(\frac{N_1}{N_2}{)}^2\frac{m_1{\omega }_1^2}{m_2{\omega }_2^2}+2\frac{N_1}{N_2}-(\frac{N_1}{N_2}{)}^2\frac{g_{11}}{g_{12}}}.$$

$$\frac{g_{12}}{g_{22}}=\frac{(\frac{N_1}{N_2}\frac{g_{11}}{g_{12}}+1{)}^2}{(\frac{N_1}{N_2}{)}^2\frac{g_{11}}{g_{12}}+\frac{m_2{\omega }_2^2}{m_1{\omega }_1^2}+2\frac{N_1}{N_2}}.$$

$$\frac{g_{11}}{g_{12}}=\frac{m_1{\omega }_1^2}{m_2{\omega }_2^2}\text{or}\frac{g_{12}}{g_{22}}=\frac{m_1{\omega }_1^2}{m_2{\omega }_2^2}.$$

$$\frac{{\mu }_1}{{\mu }_2}=\frac{m_1{\omega }_1^2}{m_2{\omega }_2^2}.$$

$$\frac{g_{12}}{g_{22}}=\frac{1}{\frac{N_1}{N_2}\frac{m_2{\omega }_2^2}{m_1{\omega }_1^2}\frac{g_{11}}{g_{12}}+\frac{m_2{\omega }_2^2}{m_1{\omega }_1^2}-\frac{N_1}{N_2}}.$$

$${\rho }_i^n(x,t)={\alpha }_{i,0}^n(t)-{\alpha }_i^n(t)x^2,u_i^n(x,t)={\beta }_i^n(t)x,$$

$$-{\dot{\alpha }}_{i,0}^n={\alpha }_{i,0}^n{\beta }_i^n,-{\dot{\alpha }}_i^n=3{\alpha }_i^n{\beta }_i^n,$$

$$\begin{array}{rl}& -{\dot{\beta }}_i^o=({\beta }_i^o{)}^2-2\frac{g_{ii}}{m_i}{\alpha }_i^o-2\frac{g_{ij}}{m_i}{\alpha }_j^o+{\omega }_i(t{)}^2,\\ & -{\dot{\beta }}_i^s=({\beta }_i^s{)}^2-2\frac{g_{ii}}{m_i}{\alpha }_i^s+{\omega }_i(t{)}^2\end{array}$$

$${\alpha }_i^n=\frac{m_i}{2g_{ii}({\zeta }_i^n{)}^3}.$$

$$\begin{array}{r}{\ddot{\zeta }}_i^o=\frac{1}{({\zeta }_i^o{)}^2}+\frac{m_j}{m_i}\frac{g_{ij}}{g_{ii}}\frac{{\zeta }_i^o}{({\zeta }_j^o{)}^3}-{\omega }_i(t{)}^2{\zeta }_i^o,{\ddot{\zeta }}_i^s=\frac{1}{({\zeta }_i^s{)}^2}-{\omega }_i(t{)}^2{\zeta }_i^s.\end{array}$$

$$\begin{array}{rl}& {\int }_0^{R_2^o}{\rho }_1^{o}dx+{\int }_{R_2^o}^{R_1}{\rho }_1^{o}dx=\frac{1}{2}{N}_1,\\ & {\int }_{R_2^o}^{R_1}{\rho }_2^{o}dx+{\int }_{R_1}^{R_2^s}{\rho }_2^{s}dx=\frac{1}{2}{N}_2;\\ & {\rho }_2^o(R_1,t)={\rho }_2^s(R_1,t).\end{array}$$

Full text

1 Introduction

The dynamics of a Bose-Einstein condensate is the subject of active current research. A multitude of experimental and theoretical works aimed at studying the solitons, vortices, dispersive shock waves, and other structures that determine the characteristic features of the behavior of a condensate in various experimental conditions have been performed by now (see, e.g., [1]). One of the main problems referring to this direction of research is to study the expansion dynamics of a condensate after the trap confining this condensate has been switched off, because in many experiments the results are recorded after the condensate expansion to a state when the cloud sizes are large enough for the measurements to be made (see, e.g., the experiments in [2, 3]). This problem was first investigated theoretically in [4] in the hydrodynamic approximation, where the equations admit a simple self-similar solution. This approach was then developed in [5, 6, 7, 8] for a condensate consisting of one component, and good agreement between theory and experiment was found. However, the situation changes significantly for the case of a condensate consisting of several components, where, for example, atoms of two different species (see [9, 10]), two different isotopes of one species of atoms (see [11]), or one species of atoms in two different quantum states, such that the difference of the energy levels of these states is much smaller than the condensate temperature (see [12, 13, 14]), are condensed. In particular, in two-component condensates the cases of relatively strong mutual repulsion between the components, where they are immiscible, and relatively weak mutual repulsion, where they are miscible, i.e., occupy the same volume, should be distinguished. This difference between the condensate phase states affects significantly the dynamics of the condensate, including the dynamics of its expansion. So far this dynamics has not been studied comprehensively enough. Some partial results illustrating the difference between the expansions of one-component and two-component condensates were obtained in [15]. However, in this paper the author used predominantly numerical methods. In our paper we will show that there are interesting situations where a comprehensive analytical study can be carried out within the hydrodynamic approximation used previously in the one-component case in [4, 5, 6, 7, 8]. Comparison with numerical calculations shows that although the dispersion effects can play some role at intermediate expansion stages, nevertheless, these effects become extremely small and they may be neglected at the asymptotic stage, which is most interesting for an experiment. The assumption that the evolution of the component density and velocity profiles is self-similar, which generalizes the approach from [4], plays a significant role in such cases favorable for the analytical theory. However, we will show that this assumption does not always adequately describe the dynamics if the initial state of the condensates before their release from the trap is near the boundary of the phase transition between component miscibility and immiscibility, and the problem requires a numerical solution in this case. Nevertheless, even in the case of immiscible components one can distinguish a characteristic case where one of the components may be considered as a piston moving the other component. For such an idealized situation the condensate expansion is accompanied by the formation of a dispersive shock wave in one component and a rarefaction wave in the other component. The theory developed for this case agrees well with the numerical results. The results of this paper allow the characteristic features of the phenomenon depending on the condensate parameters to be predicted.

2 The hydrodynamical form of the Gross-Pitaevskii equations

The dynamics of a Bose-Einstein condensate under the action of a potential $U$ is described with a high accuracy by the Gross-Pitaevskii equations. In the two-component case that we will be concerned with here, these can be written as

[TABLE]

where $i,j=1,2\leavevmode\nobreak\ (i\neq j)$ is the number of the corresponding condensate component, $(\psi_{1},\psi_{2})$ are the wave functions of the components, $g_{ii}$ are the interaction constants between atoms of component $i$, and $g_{ij}$ are the interaction constants between atoms of different species. Usually, $g_{12}=g_{21}$, which we will assume in the subsequent discussion. The interaction constants can be expressed via the scattering lengths $a_{ij}$ of atoms by one another as

[TABLE]

where $m_{ij}^{-1}=m_{i}^{-1}+m_{j}^{-1}$ is the reduced mass of the atoms being scattered by one another. Each of the wave functions is normalized to the number of particles of a given species in the condensate:

[TABLE]

so that $|\psi_{i}|^{2}=\rho_{i}$ is the number density of atoms in the $i$th component. The gradient of the phase $\varphi_{i}$ of the wave function $\psi_{i}=\sqrt{\rho_{i}}\exp(i\varphi_{i})$ is related to the flow velocity $\textbf{u}_{i}$ of $i$th component by the relation (see [1])

[TABLE]

The condensate components are miscible, i.e. their uniform distribution over space (in the absence of an external potential) is stable, if the interaction constants satisfy the condition (see [16])

[TABLE]

If, alternatively, the sign of this inequality is opposite, then the condensate is unstable with respect to the separation into regions containing the components of only one of the condensate species. However, this condition is valid only for a uniform distribution of the condensate. In contrast, for the case of a condensate confined in a trap, the miscibility condition requires a modification, which we will dwell on in more detail in the next section of our paper.

If the phase $\varphi_{i}$ is a single-valued function of coordinates, which physically means the absence of vortices in the condensate, then the wave functions of the two-component condensate can be represented as

[TABLE]

where $\mu_{i}$ is the chemical potential of the $i$th component (see [1]). Substituting (6) into (1), separating the real and imaginary parts, and differentiating one of the equations with respect to $r$ bring the Gross-Pitaevskii equations to the so-called hydrodynamic form:

[TABLE]

[TABLE]

Equations (7) are responsible for the conservation of the number of particles in the corresponding condensate component. If there were no last term proportional to $\hbar^{2}$ in Eqs. (8), then these equations would correspond to ordinary Eulerian hydrodynamics with a pressure gradient $\nabla p_{i}=(g_{ii}\nabla\rho_{i}+g_{ij}\nabla\rho_{j})/m_{i}$. However, the last term of Eqs. (8) attributable to the dispersion of quantum particles introduces new properties if the condensate characteristics change rapidly enough. Let us make an estimate for the distance $\xi$ at which the pressure and dispersion contribute identically. We assume that the masses of atoms and the number densities of particles in the components are of the same order of magnitude ($m_{i}\sim m_{j}$, $\rho_{i}\sim\rho_{j}$) for both components, so that by $m$, $\rho$ and $g$ we can understand the corresponding parameter of any component. We will then estimate the pressure as $p\sim g\rho^{2}/2m$ and obtain $g\rho\sim\hbar^{2}/(m\xi^{2})$, for $\xi$, whence $\xi\sim\hbar/\sqrt{g\rho m}$. Thus, the condensate has an intrinsic characteristic size that is called the correlation length and can be defined as

[TABLE]

If the characteristics change weakly at distances $\sim\xi$, then the last term in Eqs. (8) can be neglected, and the system will then take the form

[TABLE]

[TABLE]

Equations (11) correspond to Eulerian hydrodynamics, and this form of hydrodynamic equations describes fairly smooth solutions, in particular, the component density distributions in the trap before it is switched off. If, however, solitions with sizes $\sim\xi$, are generated or dispersive shock waves are formed in the course of evolution, then the dispersion effects should also be taken into account. Such a case will be considered in Section 5.

To describe the characteristic features of the phenomenon, we will dwell on the example of traps in which the motion of particles in two directions is frozen and is reduced to zero-point oscillations. In an experiment such a quasi-one-dimensional condensate acquires a highly elongated cigar shape. The potential of such a trap for the $i$th component can be written as

[TABLE]

where $\omega_{y}=\omega_{z}\equiv\omega_{\bot}\gg\omega_{x}$. Owing to the latter inequality, the motion of the condensate in the transverse direction is frozen, i.e., the transverse wave function is reduced to the ground state in a transverse potential with frequency $\omega_{\bot}$. The Gross-Pitaevskii equation can then be averaged over the transverse direction, and the dynamics of the condensate is reduced to its motion in the longitudinal $x$ direction (for details, see, e.g., [8]). Introducing the effective nonlinear constants of longitudinal condensate dynamics (but retaining, for simplicity, the previous notation for them), we arrive at the equation

[TABLE]

while Eqs. (10) and (11) will transform to

[TABLE]

[TABLE]

Having established the basic equations for the dynamics of a binary condensate, let us first turn to the problem of classifying the possible initial density distributions of the components before their release from the trap.

3 The phase diagram for a binary condensate confined in a quasi-one-dimensional trap

Numerical calculations (see, e.g., [17, 18]) and experimental works (see, e.g., [19]) show that various particle number density profiles can be realized, depending on the relation between the interaction constants, particle masses, numbers of particles in each component, and trap frequencies. Obviously, the condensate loaded into a trap will be distributed over space so as to minimize the total energy

[TABLE]

These distributions can have different forms, depending on the nonlinear constants and trap frequencies, and in this section we will classify the possible forms in the Thomas-Fermi approximation, where the dispersion properties of the condensate may be neglected:

[TABLE]

This approximation will allow us to establish the main types of possible distributions on a qualitative level. For the Thomas-Fermi approximation to be applicable, the size of each condensate cloud must be much greater than the correlation length $\xi$ what we will assume below. First of all, note that in the distribution there can be regions of space where both components are present (“overlap” regions) and regions where only one of the components is present (“singlet” region). Therefore, let us write out the stationary solution of Eqs. (14) and (15) that corresponds to the Thomas-Fermi approximation for these two possible cases:

[TABLE]

where $\rho_{i}(x)=|\psi_{i}(x,0)|^{2}$. Here, the index “$o$” (overlap) denotes the particle number density in the overlap region and the index “$s$” (singular) denotes the densities in the singlet regions where only one of the components is located. The chemical potentials $\mu_{i}$ are functions of the number of particles, particle masses, interaction constants, and trap frequencies. These functions are defined by the equations

[TABLE]

where $N_{1}$ and $N_{2}$ are the numbers of particles in the first and second components, respectively, and the integration is over the region where the corresponding condensate component is located. It is easy to find the sizes of each of the components of the condensate confined in the trap from Eqs. (18) and (19). Following [17], we will classify the possible configurations by associating them with points on the plane with coordinate axes $(g_{11}/g_{12},g_{12}/g_{22})$. These variables characterize the relative value of the interaction constants. The phase diagram arising in this way is shown in Fig. 1, while the typical distributions corresponding to the points on this plane are shown in Fig. 2. Let us introduce the following terminology for the various phases that can be identified in these figures. We will call the configuration where both components overlap at the trap center the miscibility phase (Fig. 22-2,2-2) and the configuration where the components are separated and one of the components is surrounded by the other the symmetric immiscibility phase (Fig. 2, 2-2). The lines in Fig. 1 separate the regions by the following attributes. It follows from (5) that the diagonal $g_{11}/g_{12}=g_{12}/g_{22}$ (dotted line) separates the miscibility/immiscibility regions of a homogeneous condensate: the components in the homogeneous condensate are immiscible above the diagonal and miscible below it. Because of the influence of the trap, this line now plays a slightly different role—it separates the condensates with nonzero (Fig. 2 and 2) and zero (Fig. 2 and 2) widths of the overlap region. These two sets of figures differ by the numbers of the components located at the trap center and outside the central component (in the “shell”).

On the solid lines the density of the external component becomes zero at the trap center, i.e., according to our definition, these lines separate the configurations with miscibility and immiscibility. The equation for the curves on which the first component forms an external shell and its density becomes zero at the trap center is analytically expressed by the formula $g_{12}/g_{22}=\mu_{1}/\mu_{2}$ which can be derived from the condition $\rho_{1}^{o}(0)=0$, and the other branch for which the first component is internal and the second one is external is defined in a similar way. The ratio of the chemical potentials can be found from system (20), and we obtain the following equation for the first component to become zero:

[TABLE]

The second component becomes zero at the trap center when passing through the curve

[TABLE]

The solid lines in Fig. 1 indicate examples of the curves, where the external component at the center of symmetry becomes zero in the Thomas-Fermi approximation, for the number of particles in the second component that is twice that in the first one and identical masses and trap frequencies. The passage through these lines is illustrated by a qualitative difference between the distributions in Figs. 2 and 2 (the first component at the center) and Fig.2 and 2 (the second component at the center).

The distributions in 2 and 2 for the external components are concave at the trap center. However, as one recedes from the solid curves, the distributions of the external components at the center become flatter and, at some moment, they become horizontal. The equation $\partial\rho_{1}^{o}(x)/\partial x=0$ gives the condition for the first component being constant, while the equation $\partial\rho_{2}^{o}(x)/\partial x=0$ gives the condition for the second component being constant. Accordingly, we find the relations between the constants of our problem:

[TABLE]

These straight lines are indicated in Fig. 1 by the dashed lines parallel to the coordinate axes. In particular, below the line $g_{12}/g_{22}=m_{1}\omega_{1}^{2}/m_{2}\omega_{2}^{2}$ the first component has an upward-convex distribution (the dashed line in Fig. 2) in the overlap region, the distribution becomes flat on this line (Fig. 2), and slightly above this line it becomes downward-concave (Fig. 2).

The dash-dotted curve in Fig. 1 separates the diagram into two regions: the region where the first component is external, while the second one is internal (the region above the dash-dotted curve), as shown in Fig. 22-2, and the region where the second component is external, while the first one is internal (the region below the dash-dotted curve), as shown in (Fig. 22-2). Equating the coordinates where the densities of the external and internal components become zero, we will obtain the following relation for the miscible components:

[TABLE]

From this condition and Eqs. (20) we will derive the following equation for the dash-dotted curve for the region where $g_{12}^{2}<g_{11}g_{22}$:

[TABLE]

As we see, it is a hyperbola in the $(g_{11}/g_{12},g_{12}/g_{22})$ plane. In the case of immiscible components, comparison of the energies for symmetric configurations shows (see Fig. 3) that the internal and external components also change places (see Figs. 2 and 2) when passing through the hyperbola on which the equality $g_{22}=g_{11}$ holds. Our numerical calculations show that this hyperbola does not depend on the number of particles in the components. This completes the construction of a phase diagram in the Thomas-Fermi approximation. In general terms,

In general terms, the constructed diagram gives a correct idea of the pattern of the component distributions in traps, except for the region near the part of the hyperbola that separates the distributions of types Fig. 2 and 2. The point is that on this curve not only the energies of the symmetric distributions in Fig. 2 and 2 but also the energy of the asymmetric distribution, where the components are on different sides of the trap center (see Fig. 2), are equal to the same value. As a result of such a degeneracy of the energies, which is illustrated in Fig. 3, even a small perturbation makes one of the distributions energetically more favorable. As our numerical calculations show, allowance for the dispersion gives an advantage to the asymmetric phase in Fig. 2. This difference is not captured by the Thomas-Fermi approximation and requires a more accurate calculation. The above classification of the possible initial states that the condensate has before the trap is switched off is sufficient for our purposes.

It should be noted that the distributions found have breaks at the transition points from the overlap regions to the singlet ones. Clearly, the dispersion effects at these points also become significant and lead to a smoothing of the curves. In particular, on the solid curves the sharpening in the distribution at the center, where the density of the second component in the Thomas-Fermi approximation is zero, is smoothed out, and a numerical solution of the Gross-Pitaevskii equation gives a relatively small, but nonzero density at the center (see, e.g., [15]).

Fig. 1 shows a general phase diagram for identical masses and trap frequencies. The number of particles in the second component is twice that in the first one. If we change the ratio of the trap frequencies and masses, then the point of intersection between the perpendicular straight lines will move along the diagonal $g_{11}/g_{12}=g_{12}/g_{22}$. For example, the point of intersection will move upward as the parameters of the first component increase and downward as the constants of the second component increase. When changing the number of particles, the point of intersection between the straight lines will be stationary, but the dash-dotted curve the passage through which interchanges the internal and external components in the miscibility region will change. In contrast, for immiscible components the dash-dotted curve will remain unchanged. The regions between the curves where the external component becomes zero at the center of symmetry and the diagonal will also change. In particular, as the number of particles in the first component increases, the region where the second component is expelled from the trap center will grow, while the region where the first component becomes zero will be reduced.

As a result, we have arrived at a complete classification of the possible initial distributions and can now turn to our main problem on the condensate expansion after the trap has been switched off.

4 A self-similar solution for the condensate expansion dynamics

As was found in [4, 5, 6, 7, 8], during the expansion of a one-component condensate it can be assumed with a good accuracy that the dependence of the density distribution on the spatial coordinate does not change in pattern, and the entire time dependence consists only in the evolution of the parameters of this distribution and the emergence of a distribution of the flow velocity proportional to the coordinate. As a result, the problem can be reduced to the solution of ordinary differential equations for the distribution parameters, and the solution can be found in a closed form in the most interesting characteristic cases. In the two-component case, this approach can have only a limited applicability. For example, if the repulsion between atoms in the internal component is much greater than the interaction forces between atoms in the external one, so that the initial distribution of the internal component is shaped mainly by the trap potential, then after the trap has been switched off, the pressure in the internal component will be a dominant force and the internal component will act on the external one like a piston . Nevertheless, if the difference between the parameters of the two components is not too large, then the time-evolving distributions will retain their initial shape with a good accuracy during the expansion, and, as in the one-component case, the problem can be reduced to solving the equations for the distribution parameters. The condition for this approximation to be applicable is that each component evolves predominantly under the action of its own pressure. In addition, if the components are separated, then the condition of mechanical equilibrium at the boundary between them must be fulfilled. This means that mechanical equilibrium is established in a time much shorter than the characteristic expansion time until the stage of motion by inertia, i.e., $R/c_{s}\ll 1/\omega$, were $R$ is the characteristic size of the condensate, and $c_{s}=\sqrt{g\rho/m}$ is the sound velocity in the condensate component. We will begin our discussion of the expansion dynamics precisely with this case.

Thus, we will seek a time-dependent solution in a form analogous to the initial distributions (18) and (19). More specifically, suppose that the density and flow velocity depend on the coordinate quadratically and linearly, respectively, with time-dependent coefficients:

[TABLE]

where $n=o,s$. Here, as before, the indices o and s denote the quantities corresponding to the overlap and singlet regions, respectively.

Substituting (26) into the continuity equation (14) and the Euler equation (15) gives

[TABLE]

[TABLE]

(the dot denotes a time derivative). Introducing a new variable $\zeta_{i}^{n}$ defined by the relation

[TABLE]

we simplify considerably these equations: from (27) we find that $\beta_{i}^{n}=\dot{\zeta}_{i}^{n}/\zeta_{i}^{n}$, and Eqs. (28) will then take the form

[TABLE]

This system of six second-order differential equations defines the motion of the Bose-Einstein condensate components. The last terms of the equations reflect the influence of the confining potential on the motion, while the other terms arise from the interaction between atoms. In what follows, we will be interested in the condensate dynamics after the trap has been switched off, i.e., we should set $\omega_{i}(t)=0$ at $t>0$. The initial conditions for these equations are determined by the original configurations that were described in the previous section. At fixed nonlinear constants they depend on the trap parameters and the number of particles in each component.

Having solved the system of equations (30), we can find the velocities $u_{i}^{n}(x,t)$ and coefficients $\alpha_{i}^{n}(t)$. To find $\alpha_{i,0}^{n}(t)$, we will use the normalization of the wave functions (20) and the fact that at the boundary between the singlet and overlap regions the pressures in them and, consequently, the densities are equal. To be specific, we number the condensate components in such a way that the first component is surrounded by the second one, as shown in Fig. 4 (the first and second components are indicated by the solid and dashed lines, respectively), so that we have

[TABLE]

Here $N_{1}$ and $N_{2}$ are are the numbers of particles in the first and second components, respectively, $R_{1}$ is the coordinate where the first component becomes zero ($\rho_{1}^{o}(R_{1},t)=\alpha_{1,0}^{o}(t)-\alpha_{1}^{o}(t)R_{1}^{2}=0$), $R_{2}^{o}$ and $R_{2}^{s}$ are the coordinates where the second condensate component becomes zero ($\rho_{2}^{n}(R_{2}^{n},t)=\alpha_{2,0}^{n}(t)-\alpha_{2}^{n}(t){R_{2}^{n}}^{2}=0$). The point $R_{2}^{o}$ corresponds to zero density of the second component in the overlap region, while $R_{2}^{s}$ denotes the coordinate at which the density of the second component becomes zero in the singlet region (see Fig. 4). Consequently, these parameters are defined by the relations

[TABLE]

The coordinates $R_{1}$, $R_{2}^{o}$ and $R_{2}^{s}$ are functions of time. Thus, we have reduced the problem to integrating the ordinary differential equations (30) with their initial conditions determined by the original component density distributions in the trap. In general, this system must be solved numerically, which is considerably easier than the solution of the Gross-Pitaevskii equations. However, an important case where system (30) is simplified considerably and admits the derivation of some relations in a closed form is noteworthy.

4.1 The Case of Miscibility

Let we have an initial overlap configuration where the first component has no singlet region (see Fig. 2-2 and Fig. 2-2), i.e., the components are miscible. System (31) will then be written as

[TABLE]

Here, $R_{2}$ is the coordinate where the particle number density of the second component becomes zero (see Fig. 4). In this case, the relations for $\alpha_{i,0}^{n}(t)$ expressed via $\alpha_{i}^{n}(t)$ can be found analytically, and the solution of this system will be

[TABLE]

The case of the expansion of a condensate with an asymmetric initial profile (see Fig. 2) will be considered separately below.

The differential equations (30) are Newton-type equations that have the total energy of the system (16) the conservation of the number of particles in the singlet region and the overlap region of each of the components as the integrals of motion. Generally, these integrals are not enough to find the analytical solution of the system. Therefore, we turn to its numerical solution. Figure 5 shows examples of comparing the numerical solution of the Gross-Pitaevskii equations (13) (dashed curves) with the numerical solution of the equations of motion (30) (solid curves) for different interaction constants at fixed times. As can be seen from the figure, in the region where the components are well miscible the particle density distributions retain their shape during the expansion, and the self-similar solution quantitatively describes the dynamics of the system excellently. A change in the nonlinear interaction constant between the components by a factor of $10$ does not affect significantly the accuracy of the approximation as long as the miscibility criterion (5) holds with a margin (in the case of Fig. 5, $g_{12}^{2}=4$ is smaller than $g_{11}g_{22}=12$ by a factor of $3$).

In practice the interaction constants have almost the same value in many cases. For example, for an 87Rb atom in different states of the hyperfine structure ($|1,-1\rangle$ and $|2,-2\rangle$) the scattering lengths are $a_{11}=98.98a_{0}$, $a_{12}=98.98a_{0}$ and $a_{22}=100.4a_{0}$, where $a_{0}$ the Bohr radius (see, e.g., [21, 22]) i.e., the interaction constants are also equal, $g_{11}=g_{12}$. The component masses and trap frequencies can often be also equal ($m_{1}=m_{2}$, $\omega_{1}=\omega_{2}$). The second (external) component of the Bose-Einstein condensate will then be constant in the overlap region and, consequently, $\alpha_{2}^{o}(t)=0$. Consider a slightly more general case where the condition (23) is fulfilled. The self-similar solution will then be written as

[TABLE]

After the substitution into Eqs. (10) and (11) the self-similar solution will give equations analogous to (27) and (28):

[TABLE]

By substituting $\alpha_{1}^{o}={m_{1}}/{2g_{11}(\zeta_{1}^{o})^{3}}$ and $\alpha_{2}^{s}={m_{2}}/{2g_{22}(\zeta_{2}^{s})^{3}}$, the system of equations at $\omega_{i}(t)=0$ is reduced to

[TABLE]

From the first equation (38) we will find

[TABLE]

where we introduce an effective frequency of the potential in which the first condensate component is located:

[TABLE]

While from the second equation (38) we obtain

[TABLE]

Equations (39) and (40) implicitly specify $\zeta_{i}^{n}$ as a function of $t$. The densities derived from these equations and by numerically solving the Gross-Pitaevskii equations are compared in Fig. 6.

The velocities expressed via $\zeta_{i}^{n}$ can also be easily found from system (38):

[TABLE]

At $t\gg\omega_{i}^{-1}$ from (39) and (40) we find the asymptotic formulas

[TABLE]

These solutions correspond to the motion of condensate atoms by inertia, when the density in the condensate cloud becomes so low that the pressure no longer accelerates the condensate.

To find the velocity of the second component in the overlap region, we will use (42) and the second equation (36) and find a differential equation for $\beta_{2}^{o}$

[TABLE]

From this equation we will obtain

[TABLE]

Knowing the asymptotic solutions for $\zeta_{i}$ and the velocities expressed via these $\zeta_{i}$ we can easily find the asymptotic solution for the velocities of each of the components $u_{i}^{n}$:

[TABLE]

The extreme points of the distribution of each condensate component will move with the greatest velocities:

[TABLE]

The first formula (46) corresponds to the velocity of the distribution boundary for the first condensate component, when the motion occurs by inertia, while the second formula represents the same velocity for the second component. From the asymptotics found for $\zeta_{i}$ we can also derive simple formulas for the density distributions at $t\gg\omega_{i}^{-1}$:

[TABLE]

Thus, in the case of strong miscibility of the condensate components with initial distributions like those in Fig. 2,2,2,2, the self-similar solution gives quite a satisfactory description of the condensate expansion after the trap has been switched off. If, however, we approach the miscibility-immiscibility boundary with initial distributions like those in Fig. 2,2,2,2, then the Thomas-Fermi approximation loses its accuracy even when calculating the stationary distributions due to the appearance of large jumps in derivatives in the density distributions. The expansion dynamics in such cases also differs significantly from the predictions of the self-similar theory. In particular, characteristic regions of nonlinear oscillations can be formed in the density and velocity distributions during the evolution. This means that wave breaking occurs under their deformation due to nonlinear effects, so that the dispersion effects can no longer be neglected. In contrast, allowance for the simultaneously nonlinear and dispersion effects gives rise to a region of oscillations connecting the regions of flows with different parameters. These regions of oscillations are analogous to shock waves in low-dissipation systems, and they were called dispersive shock waves . However, before turning to their study, let us consider the case where the components are strongly immiscible, which is also described by a self-similar solution analogous to those obtained above for miscible components.

4.2 The Case of Immiscible Components

Here, we assume that the Bose-Einstein condensate components are absolutely immiscible, i.e., the overlap region is of the order of the correlation length. As has been pointed out in Section 3, symmetric (Fig. 2,2) distributions can be realized in this case. The self-similar solution will now take the form

[TABLE]

and out of Eqs. (30) only those responsible for the singlet regions will remain. The initial distributions $\alpha_{i,0}^{s}$ and $\alpha_{i}^{s}$ are defined by the stationary solutions

[TABLE]

At the boundary between the immiscible components the pressures in them are equal, $p_{i}=p_{j}$, which serves as a condition defining the boundary coordinate. In the hydrodynamic approximation the pressures are $p_{i}=g_{ii}{\rho_{i}^{s}}^{2}/2m_{i}$. So, from this condition we find that the coordinate $R_{b}$ of the boundary between the components is given by the formula

[TABLE]

Thus, the problem has been reduced to determining six functions of time that are the coefficients in (50).

The expressions for $\alpha_{i}^{s}(t)$ can be found in a closed form. At $t>0$, where $\omega_{i}(t)=0$, we will obtain simple equations of motion for the parameters of the singlet distributions:

[TABLE]

($\alpha_{i}^{s}=m_{i}/2g_{ii}(\zeta_{i}^{s})^{3}$). It is easy to find the solution of this system:

[TABLE]

These equations implicitly specify $\zeta_{i}^{s}$ as functions of $t$. The functions $\beta_{i}^{s}(t)={\dot{\zeta}_{i}^{s}}/{\zeta_{i}^{s}}$ are expressed via the functions $\zeta_{i}^{s}$ found. As a result, we obtain the expressions for the velocity of each component $u_{i}^{s}$ expressed via $\zeta_{i}^{s}$:

[TABLE]

Finally, the relation between the parameters $\alpha_{i,0}^{s}$ and $\alpha_{i}^{s}$ follows from the normalization conditions (here as before $R_{1}=\sqrt{\alpha_{1,0}^{s}/\alpha_{1}^{s}}$ and $R_{2}=\sqrt{\alpha_{2,0}^{s}/\alpha_{2}^{s}}$):

[TABLE]

for the symmetric profile and from the equations

[TABLE]

for the asymmetric one. In principle, these equations allow $\alpha_{i,0}^{s}$ to be expressed via the already known functions $\alpha_{i}^{s}$ through numerically easily solvable algebraic equations (which we do not write out here, because they are cumbersome).

As in the case of $g_{11}=g_{12}$, from the system of equations (54) we can find the asymptotic solution for $\zeta_{i}^{s}$ corresponding to times $t\gg\omega_{i}^{-1}$:

[TABLE]

As above, they describe the flow of the condensate cloud by inertia once the potential energy of the condensate compressed in the trap has been converted into the kinetic energy of its flow. The extreme points of the distributions of each condensate component move with maximum velocities:

[TABLE]

for the symmetric and asymmetric particle number density profiles, respectively. The expressions for the particle number densities and velocities at asymptotically large $t$ are

[TABLE]

The functions $\alpha_{i,0}^{s}(t)$ The functions (56) in the symmetric case and (57) in the asymmetric case for each time $t$. Knowing $\alpha_{1,0}^{s}$ and $\alpha_{2,0}^{s}$, we can find the forms of the particle number densities and the maximum velocities of each of the condensate components at this time.

Thus, the derived formulas give the solution of our problem. It is illustrated in Fig. 7, where the particle number densities are shown for the symmetric and asymmetric phases at a fixed time.

5 Dispersive shock wave

A characteristic feature of the expansion of a twocomponent Bose-Einstein condensate is the possibility of the appearance of dispersive shock waves in it (see [23, 24]). Fig. 8 shows an example of such waves, where the internal component expels the external one so strongly that this leads to a breaking of the particle number density distribution of the external component, which implies the beginning of the formation of a dispersive shock wave. In the case where the external component is much larger than shock wave size, we can roughly assume that in the place of shock wave formation the external component is homogeneous, while the boundary of the internal component is a piston pushing the external component. Similarly, the formation of dispersive shock wave is also possible for immiscible components. For example, in a condensate confined in a trap a boundary is formed between the condensates at which the pressures of both components are equal. If the frequencies of the traps confining the various components differ greatly, then a pressure jump arises at the boundary between the condensates immediately after the traps have been switched off, so that one component will again act on the other component as a piston. At the initial stage the region near the boundary between the condensates is of greatest interest, while the condensate parameters far from the boundary may be deemed to be specified by constants. Such problems on a piston both with a constant velocity of its motion and with a uniformly accelerated or arbitrary law of its motion were considered in [25, 26] and can be applied to the problem on the condensate expansion after the traps have been switched off.

As an example, consider an immiscible condensate at an evolutionary stage when the size of the region occupied by the dispersive shock wave is much smaller than the size of the entire condensate cloud. We will also assume that the internal component plays the role of a piston. The evolution of the external component is then described by a solution analogous to those found in [25, 26], but with allowance made for the fact that now the piston velocity is not specified but must be found in a self-consistent way from the condition that the pressures of the two components at the boundary between them are equal. This means that for the internal component we are also dealing with the piston problem, but, in this case, the piston stretches the condensate rather than compresses it. As a result, a rarefaction wave described by the well-known solution of the hydrodynamic equations (14) and (15). for the internal component will propagate into the internal condensate. A characteristic density profile of the emerging wave structure is shown in Fig. 9 it was obtained by numerically solving the Gross-Pitaevskii equations (13). Our task in this section is to calculate the main characteristics of this structure by expressing them via the condensate parameters and initial conditions.

Thus, we will assume that at the initial time, immediately after the traps have been switched off, the first component is located to the left of the interface at $x=0$ and has a density $\rho=\rho_{1,0}$, while the second component is located to the right of the interface and has a density $\rho=\rho_{2,0}$. The flow velocities of both components are zero, $u_{1}=u_{2}=0$, at the initial time $t=0$. To be specific, we will also assume that the pressure to the left of the interface at the initial time is greater than the pressure to the right of the interface:

[TABLE]

Precisely such a condition corresponds to the formation of the wave structure shown in Fig. 9. Clearly, this condition does not limit the generality of our consideration: if the inequality is opposite, then the waves will propagate in the opposite direction, and it will be sufficient to replace $x$ by $-x$ in the formulas derived below in order that they describe this case.

In accordance with Fig. 9 it follows from the theory of compressible fluid motion [27] that the left rarefaction wave connects two regions of homogeneous flow: the condensate at rest with a density $\rho_{1,0}$ and the “shel” to the right of the rarefaction wave with some density $\overline{\rho}_{1}$ and flow velocity $\overline{u}_{1}$. This shelf extends to the boundary with the second condensate that moves with a velocity $\overline{u}$ equal to the flow velocity $\overline{u}_{1}$ of the left condensate. As has been pointed out above, the boundary between the components serves as a “piston” for the second condensate. Therefore, in accordance with the well-known solution [25], a shelf of the flow of the second condensate with a density $\overline{\rho}_{2}$ and flow velocity $\overline{u}_{2}$ again equal to the boundary velocity $\overline{u}$ also emerges rightward of the boundary. This second shelf is connected by the dispersive shock wave with the region of the second component at rest with a density $\rho_{2,0}$. The boundary between the two condensates is an analog of the contact discontinuity known from the classification of the decays of initial discontinuities in gas dynamics (see. [27]). Owing to the equality of the nonlinear constants, $g_{11}=g_{22}$ the equality of the densities follows from the equality of the pressures at such a discontinuity. For simplicity, precisely such a case is assumed in Fig. 9 The parameters of the shelfs $\overline{\rho}_{1}$, $\overline{\rho}_{2}$, the velocity of the boundary (contact discontinuity) $\overline{u}=\overline{u}_{1}=\overline{u}_{2}$, and the velocities of the rarefaction wave and dispersive shock wave boundaries are required to be found in the problem formulated in this way.

Let us first discuss this problem in terms of the dispersionless hydrodynamic equations

[TABLE]

These are a special case of Eqs. (14) and (15), where we take into account the fact that, in this case, the components are spatially separated, so that the index $i=1$ refers to the component to the left of the interface, while the index $i=2$ refers to the component to the right of this interface. There is no interaction between the components everywhere except the narrow region near the boundary separating them, but its presence is taken into account by an appropriate boundary condition for the pressures to be equal on both sides of the boundary. In addition, we assume the masses of the atoms in both components to be identical and equal to unity, while the trap potential at the expansion stage to play no role in an obvious way. The solution of such problems is simplified considerably if we pass from the ordinary physical variables $\rho_{i}$, $u_{i}$ to the so-called “Riemann invariants”. For Eqs. (63), which coincide in form with the gas-dynamic equations, the Riemann invariants are well known and can be written as (see, e.g., [28])

[TABLE]

In these variables the hydrodynamic equations (63) take a simple symmetric form:

[TABLE]

where the ”Riemann velocities $v_{\pm}^{(i)}=u_{i}\pm c_{i}$ are expressed via the Riemann invariants by the relations

[TABLE]

If the solution of Eqs. (65) has been found and the Riemann invariants are known, then the physical variables are expressed via them by the formulas

[TABLE]

A rarefaction wave is known [27] to belong to the class of “simple waves” characterized by the fact that one of the Riemann invariants has a constant value along the flow. For the case shown in Fig. 9, the rarefaction wave propagates to the left in the first condensate. Hence, the following invariant is constant for it:

[TABLE]

where we set its value equal to the value at the boundary with the condensate at rest. It must have the same value at the rarefaction wave boundary with the shelf to the right of it:

[TABLE]

The pressures of the components at the boundary between them are equal, which gives the relation

[TABLE]

Finally, after the passage through the dispersive shock wave in the second component, the Riemann invariant $r_{-}^{(2)}$ retains its value, which gives the relation

[TABLE]

The three equations (69)-(71) allow the densities of the components on the shelves and the boundary velocity to be found:

[TABLE]

[TABLE]

When inequality (62) holds, we have $\overline{u}>0$, i.e., the boundary between the components moves to the right, as it must be. The derived formulas are simplified considerably if the nonlinear constants are equal, $g_{11}=g_{22}\equiv g$:

[TABLE]

In this case, the wave amplitude $\sqrt{\rho}$ on the plateau between the rarefaction wave and the dispersive shock wave is equal to the arithmetic mean of the component amplitudes in the initial state, while the velocity of the boundary between the components is equal to the difference of the sound velocities in the initial states of the components.

Having established the boundary values of the parameters on both sides of the rarefaction wave and the dispersive shock wave, we can turn to finding the solutions for the waves themselves. We will begin with the rarefaction wave. As we know, the Riemann invariant $r_{+}^{(1)}$ is constant along it (see (68)), so that Eqs. (65) is satisfied identically for this invariant. The second Riemann invariant $r_{-}^{(1)}$ depends only on the self-similar variable $\zeta=x/t$, because our initial condition in the form of a pressure jump at the boundary between the components contains no parameter with the dimensions of length. Therefore, Eq. (65) is reduced to the equation

[TABLE]

with the obvious solution

[TABLE]

Thus, along the rarefaction wave the Riemann invariants are given by the expressions

[TABLE]

Substituting them into (67), we will find the density of the condensate and its flow velocity in the rarefaction wave:

[TABLE]

The flow velocity becomes zero at the boundary with the condensate at rest at $x_{-}^{(1)}=-\sqrt{g_{11}\rho_{1,0}}\cdot t$, i.e., the left edge of the rarefaction wave propagates into the condensate at rest with a sound velocity $s_{-}^{(1)}=-\sqrt{g_{11}\rho_{1,0}}$ in it. The right edge is joined with the shelf when $u_{1}=\overline{u}$, whence it follows that the velocity of the right edge is

[TABLE]

At $g_{11}=g_{22}=g$ velocities of the rarefaction wave edges are

[TABLE]

If we attempt to find the wave between the plateau and the second component at rest $\overline{\rho}_{2}$, $\overline{u}_{2}$ in the dispersionless approximation, then we will arrive at the multivalued solution indicated in Fig. 9 by the dashed line. This means that the dispersionless approximation is inapplicable in this case, and allowance for the dispersion leads to the replacement of the multivalued solution by a region of fast oscillations called a dispersionless shock wave. This region can be roughly described as a modulated periodic solution of the Gross-Pitaevskii equations with slowly changing parameters. The equations describing the evolution of the parameters can be derived by averaging the proper number of conservation laws. This method of deriving the modulation equations for nonlinear waves was proposed by Whitham [29, 30]; for the Korteweg-de Vries equation these equations were transformed by him to a diagonal Riemann form analogous to the gasdynamic equations (65) The method developed by Whitham was applied to the problem on the formation of dispersive shock waves by Gurevich and Pitaevskii in [31]. Since then the theory of dispersive shock waves based on Whitham s method has been elaborated greatly (see, e.g., the recent review [32]). In particular, Whitham s modulation equations were derived for periodic solutions of the Gross-Pitaevskii equation (a nonlinear Schrodinger equation) in [33, 34], the problem on the evolution of a step was analyzed in [35, 36], and the developed theory was applied to the problems on the evolution of a condensate under the action of a moving piston in [25, 26]. In our case, the evolution of the second component to the right of the interface between the components is reduced to the problem on the flow of a condensate under the action of a piston moving with velocity (73). Therefore, the results of the above papers are directly applicable to our problem, and below we will provide the basic formulas that describe the flow of the second condensate.

The periodic solution of the Gross-Pitaevskii equations in a “hydrodynamic form” can be written as

[TABLE]

where

[TABLE]

and the real parameters ?i $\lambda_{i}$ are ordered according to the inequalities

[TABLE]

As we see, the wave phase velocity $V$, amplitude $a=(\lambda_{4}-\lambda_{3})(\lambda_{2}-\lambda_{1})$, and background density $\rho_{0}=(\lambda_{4}-\lambda_{3}-\lambda_{2}+\lambda_{1})^{2}/4$ through which the wave propagates are expressed via these parameters. These parameters are slow functions of $x$ and $t$ in a dispersive shock wave. The periodic solution written in form (81) has the advantage that the parameters $\lambda_{i}$ are Riemann invariants, and their evolution is defined by Whitham s equations in a diagonal Riemann form (see [33, 34])

[TABLE]

Here, $v_{i}$ are the characteristic Whitham velocities:

[TABLE]

where $K=K(m)$, $E=E(m)$ are elliptic integrals of the first and second kinds, respectively. In the limit $m\rightarrow 1$ ($\lambda_{3}\rightarrow\lambda_{2}$) a traveling wave transforms into a soliton solution against the background of a constant density:

[TABLE]

In the other (small-amplitude) limit $m\rightarrow 0$ ($\lambda_{3}\rightarrow\lambda_{4}$ or $\lambda_{2}\rightarrow\lambda_{1}$) the wave amplitude approaches zero, while the density takes its background value. Significantly, the pair of Whitham velocities in these limits transforms (to within a constant coefficient) to the Riemann velocities of the dispersionless limit. In particular, in the soliton limit at $\lambda_{3}=\lambda_{2}$

[TABLE]

and in the small-amplitude limit needed for us at $\lambda_{3}=\lambda_{4}$

[TABLE]

This means that the edges of the dispersive shock wave are joined with the smooth solutions of the hydrodynamic approximation, i.e., in our case with the shelf $\overline{\rho}_{2}$, $\overline{u}_{2}$ at the soliton edge and with the component at rest $\rho_{2,0}$, $u_{2,0}=0$ at the small-amplitude edge.

To elucidate the behavior of the solution of Eqs. (83), we will again use the argument that the initial condition contains no parameters with the dimensions of length, so that the modulation parameters depend only on the self-similar variable $\zeta=x/t$. Therefore, Eqs. (83) are reduced to

[TABLE]

Hence it follows that only one Riemann invariant is variable, while the three remaining ones must have constant values. From the joining condition at the edges of the dispersive shock wave we find that at the soliton edge

[TABLE]

where $\overline{r}_{\pm}$ are the Riemann invariants of the dispersionless theory that are defined by Eqs. (64) and take the values on the plateau bordering the soliton edge of the dispersive shock wave. Similarly, at the small-amplitude edge we find

[TABLE]

Thus, we arrive at the dependence of the Riemann invariants on the spatial coordinate shown in Fig. 10 with the following constant Riemann invariants:

[TABLE]

In contrast, the dependence of the variable Riemann invariant $\lambda_{3}$ on $\zeta=x/t$ is implicitly defined by the equation

[TABLE]

Substituting the Riemann invariants found into the periodic solution (81) gives the coordinate dependence of the condensate density at a fixed time. It is illustrated in Fig. 11 for the same values of the parameters at which our numerical calculation shown in Fig. 9. was made. As we see, there is good agreement of the analytical results with the numerical ones.

The derived formulas also give analytical expressions for the velocities of the edges of the dispersive shock wave. The soliton edge moves with the velocity

[TABLE]

while the velocity of the small-amplitude edge is

[TABLE]

When the nonlinear constants are equal, $g_{11}=g_{22}=g$, these velocities are

[TABLE]

These values also agree well with the results of our numerical calculation.

6 Conclusions

The results of this paper show that the expansion dynamics of a two-component condensate is distinguished by a great diversity compared to the one-component case. First, the initial states of the two-component condensate can form various configurations in a nonuniform trap field, and we constructed a phase diagram of these states by refining the previous studies of other authors. Second, the simple self-similar ansatz that was successfully used in the theory of the expansion of a one-component condensate now has a limited applicability, describing in general terms the expansion of the two-component condensate only far from the transition line from the miscible components to the immiscible ones. Finally, third, if the components are immiscible, then the expansion regimes with the generation of dispersive shock waves are possible. This can be of interest both for the analysis of condensate flows in specific experimental conditions and for nonlinear physics as a whole. The results obtained in this paper allow one to predict the main characteristics of the expansion dynamics and to estimate the main parameters of the emerging wave structures.

Acknowledgments

This work was supported by the Russian Foundation for Basic Research (project no. 16-01-00398).