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

###### 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.

## 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 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

(1) |

where is the number of the corresponding condensate component, are the wave functions of the components, are the interaction constants between atoms of component , and are the interaction constants between atoms of different species. Usually, , which we will assume in the subsequent discussion. The interaction constants can be expressed via the scattering lengths of atoms by one another as

(2) |

where 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:

(3) |

so that is the number density of atoms in the th component. The gradient of the phase of the wave function is related to the flow velocity of th component by the relation (see [1])

(4) |

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])

(5) |

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 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

(6) |

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

(7) |

(8) |

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 in Eqs. (8), then these equations would correspond to ordinary Eulerian hydrodynamics with a pressure gradient . 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 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 (, ) for both components, so that by , and we can understand the corresponding parameter of any component. We will then estimate the pressure as and obtain , for , whence . Thus, the condensate has an intrinsic characteristic size that is called the correlation length and can be defined as

(9) |

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

(10) |

(11) |

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 , 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 th component can be written as

(12) |

where . 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 . 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 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

(13) |

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

(14) |

(15) |

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

(16) |

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:

(17) |

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 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:

(18) | |||

(19) |

where . Here, the index “” (overlap) denotes the particle number density in the overlap region and the index “” (singular) denotes the densities in the singlet regions where only one of the components is located. The chemical potentials are functions of the number of particles, particle masses, interaction constants, and trap frequencies. These functions are defined by the equations

(20) |

where and 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 . 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. 2\subrefTomasFermiA-\subrefTomasFermiC,\subrefTomasFermiH-\subrefTomasFermiJ) and the configuration where the components are separated and one of the components is surrounded by the other the symmetric immiscibility phase (Fig. 2, \subrefTomasFermiD-\subrefTomasFermiG). The lines in Fig. 1 separate the regions by the following attributes. It follows from (5) that the diagonal (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 which can be derived from the condition , 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:

(21) |

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

(22) |

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 gives the condition for the first component being constant, while the equation gives the condition for the second component being constant. Accordingly, we find the relations between the constants of our problem:

(23) |

These straight lines are indicated in Fig. 1 by the dashed lines parallel to the coordinate axes. In particular, below the line 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. 2\subrefTomasFermiA-\subrefTomasFermiE, 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. 2\subrefTomasFermiF-\subrefTomasFermiJ). Equating the coordinates where the densities of the external and internal components become zero, we will obtain the following relation for the miscible components:

(24) |

From this condition and Eqs. (20) we will derive the following equation for the dash-dotted curve for the region where :

(25) |

As we see, it is a hyperbola in the 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 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 . 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., , were is the characteristic size of the condensate, and 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:

(26) |

where . 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

(27) |

(28) |

(the dot denotes a time derivative). Introducing a new variable defined by the relation

(29) |

we simplify considerably these equations: from (27) we find that , and Eqs. (28) will then take the form

(30) |

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 at . 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 and coefficients . To find , 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

(31) |

Here and are are the numbers of particles in the first and second components, respectively, is the coordinate where the first component becomes zero (), and are the coordinates where the second condensate component becomes zero (). The point corresponds to zero density of the second component in the overlap region, while 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

(32) |

The coordinates , and 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-\subrefTomasFermiC and Fig. 2-\subrefTomasFermiJ), i.e., the components are miscible. System (31) will then be written as

(33) |

Here, is the coordinate where the particle number density of the second component becomes zero (see Fig. 4). In this case, the relations for expressed via can be found analytically, and the solution of this system will be

(34) |

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 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, is smaller than by a factor of ).

In practice the interaction constants have almost the same value in many cases. For example, for an Rb atom in different states of the hyperfine structure ( and ) the scattering lengths are , and , where the Bohr radius (see, e.g., [21, 22]) i.e., the interaction constants are also equal, . The component masses and trap frequencies can often be also equal (, ). The second (external) component of the Bose-Einstein condensate will then be constant in the overlap region and, consequently, . Consider a slightly more general case where the condition (23) is fulfilled. The self-similar solution will then be written as

(35) |

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

(36) |

By substituting and , the system of equations at is reduced to

(37) |

From the first equation (37) we will find

(38) |

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

While from the second equation (37) we obtain

(39) |

Equations (38) and (39) implicitly specify as a function of . The densities derived from these equations and by numerically solving the Gross-Pitaevskii equations are compared in Fig. 6.

The velocities expressed via can also be easily found from system (37):

(40) |

At from (38) and (39) we find the asymptotic formulas

(41) |

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 (41) and the second equation (36) and find a differential equation for

(42) |

From this equation we will obtain

(43) |

Knowing the asymptotic solutions for and the velocities expressed via these we can easily find the asymptotic solution for the velocities of each of the components :

(44) |

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

(45) |

The first formula (45) 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 we can also derive simple formulas for the density distributions at :

(46) | |||

(47) | |||

(48) |

Thus, in the case of strong miscibility of the condensate components with initial distributions like those in Fig. 2,\subrefTomasFermiB,\subrefTomasFermiI,\subrefTomasFermiJ, 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,\subrefTomasFermiD,\subrefTomasFermiG,\subrefTomasFermiH, 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,\subrefTomasFermiF) distributions can be realized in this case. The self-similar solution will now take the form

(49) |

and out of Eqs. (30) only those responsible for the singlet regions will remain. The initial distributions and are defined by the stationary solutions

(50) |

At the boundary between the immiscible components the pressures in them are equal, , which serves as a condition defining the boundary coordinate. In the hydrodynamic approximation the pressures are . So, from this condition we find that the coordinate of the boundary between the components is given by the formula

(51) |

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

The expressions for can be found in a closed form. At , where , we will obtain simple equations of motion for the parameters of the singlet distributions:

(52) |

(). It is easy to find the solution of this system:

(53) |

These equations implicitly specify as functions of . The functions are expressed via the functions found. As a result, we obtain the expressions for the velocity of each component expressed via :

(54) |

Finally, the relation between the parameters and follows from the normalization conditions (here as before and ):

(55) |

for the symmetric profile and from the equations

(56) |

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