Dispersive hydrodynamics of nonlinear polarization waves in two-component Bose-Einstein condensates

We study one dimensional mixtures of two-component Bose-Einstein condensates in the limit where the intra-species and inter-species interaction constants are very close. Near the mixing-demixing transition the polarization and the density dynamics decouple. We study the nonlinear polarization waves, show that they obey a universal (i.e., parameter free) dynamical description, identify a new type of algebraic soliton, explicitly write simple wave solutions, and study the Gurevich-Pitaevskii problem in this context.


Introduction
As demonstrated in various physical contexts, the interplay between dispersive and nonlinear effects can lead to a number of spectacular phenomena as, for instance, the formation of solitons and vortices. Bose-Einstein condensates (BECs) display both effects: (i) dispersion which is due to the so-called quantum pressure and (ii) nonlinear properties due to the interaction between the condensed atoms. Already in a pioneering paper, Bogoliubov [1] showed that the combination of these two features yields reconstruction of the ground state of the many-particle system, with formation of new types of elementary excitations-Bogoliubov quasiparticles. The generalization of the Bogoliubov method to nonuniform time-dependent systems by Gross [2] and Pitaevskii [3] permitted to develop the theory of quantum vortices and later Tsuzuki [4] demonstrated the existence of dark solitons in a one dimensional model of weakly interacting bosons. After the experimental realization of BEC in ultracold gases, dark solitons were observed first in a one-dimensional geometry under the form of dips propagating along a stationary background [5,6] and then in two dimensions under the form of stationary oblique solitons [7][8][9] generated by the flow of an exciton-polariton condensate past an obstacle [10,11]. More complicated nonlinear wave structures were experimentally observed [12] and interpreted as dispersive shock waves, the description of which can be developed in the framework of Whitham's theory of modulations of nonlinear waves [12,13] (for a recent review on modulation theory of nonlinear waves see, e.g., Ref. [14]).
The experimental realization of condensates consisting of two (or more) species has opened the possibility of studying even richer dynamics triggered by the additional degree(s) of freedom consisting in the relative motion of the components. These are new modes that can interact with each other leading, in particular, to different types of solitons. For two-component systems these new modes can be visualized as pertaining to two types of waves: "density waves" with in-phase motion of the two components and "polarization waves" with counterphase motion of the components. In the simplest situations these two types of excitations decouple: the first type does not involve relative motion of the components and the second type of waves does not affect the total density of the condensate. In the small amplitude limit these two types of waves and the distinction between density and polarization excitations were studied, e.g., in Ref. [15].
It has been recently noticed [16] that the polarization dynamics can be separated from the density dynamics even for the case of large amplitude waves if the difference between intraand inter-species interaction constants is small, and this observation was applied to the theory of polarization solitons-which were denoted as "magnetic solitons". In the present paper we extend the method of Ref. [16] to the general case of polarization dynamics in two-component BECs with small difference between the nonlinear interaction constants. In section 2 we derive the general equations of the polarization dynamics. In section 3 we study their traveling wave solutions that include, as a limiting case, the soliton solution found in [16] and in Sec. 4 we study the dispersionless limit of the nonlinear polarization waves. This forms the basis for discussing in section 5 the evolution of initial discontinuities in the polarization distribution. We show that such discontinuities evolve into a wave structure consisting in a rarefaction wave separated from a dispersive shock wave by a plateau with constant polarization and relative flow velocity. The main characteristics of this structure are calculated with the use of Whitham theory and are shown to compare very well with the results of numerical simulations. The relevance of our results for experimental studies is discussed in section 6. Our conclusions are summarized in Sec. 7 and some technical aspects are detailed in Appendixes A and B.
We consider a one-dimensional BEC system described by a two-component spinor order parameter Ψ(x, t) = (ψ ↑ , ψ ↓ ) t (where the superscript t denotes the transposition) obeying the following coupled Gross-Pitaevskii equations In Eqs. (1), it has been assumed that the two intra-species non linear coefficients g ↑↑ and g ↓↓ have the same value, denoted as g. For instance, this is exactly realized in the mixture of the two hyperfine states |F = 1, m F = ±1〉 of 23 Na [19], and, to a good approximation, in the mixture of hyperfine states of 87 Rb considered in Ref. [20] (|F, m F 〉 = |1, 1〉 and |2, 2〉). The inter-species coefficient g ↑↓ is written as g − δg, and we assume that Both conditions are realized in the above presented cases of 23 Na (δg/g 0.07) and 87 Rb (δg/g 0.01). The left condition is the mean-field miscibility condition of the two species (see, e.g., Refs. [21,22]). The right condition will be shown later to lead to important simplifications in the dynamics of the system.
We represent the spinor wave function as In this expression, ρ(x, t) is the total density and θ (x, t) governs the linear densities of the two components: ρ ↑ (x, t) = |ψ ↑ | 2 and ρ ↓ (x, t) = |ψ ↓ | 2 (cf. Eqs. (18) in the case of a constant total density ρ 0 ). Φ(x, t) and φ(x, t) act as potentials for the velocity fields v ↑ and v ↓ of the two components. Namely By means of the substitution (3) the Gross-Pitaevskii system (1) is cast into the form where it is assumed that at equilibrium both components are at rest and both have the same uniform density. The total density is denoted as ρ 0 . In this case the chemical potentials take the same value: µ ↑ = µ ↓ = (g − δg)ρ 0 /2. As is known, in such a system there are two types of waves that can be called "density" and "polarization" waves. In the small amplitude and long wavelength limit the velocity of polarization waves that correspond to the (mainly) relative motion of the components is equal to In the limit (2) c p is very small compared to the long wavelength velocity c d of density waves [mc 2 d = ρ 0 (g−δg/2)]. Following Ref. [16], we introduce also the "polarization healing length" Then the characteristic time scale for the polarization dynamics can be measured in units of T p and ξ p are much larger than the corresponding characteristic time and length associated with density waves, and for the study of the polarization nonlinear waves it is thus appropriate to pass to the non-dimensional variables Then a very important consequence can be inferred from the second equation (5) that, in new non-dimensional variables, can be written as We see that for θ ∼ 1 at space and time scales of order (7) and (8), correspondingly, the righthand side becomes small if δg/g 1. In this case we can assume at the leading order that ρ ≈ ρ 0 , so that the polarization hydrodynamics is decoupled from the density dynamics. This important feature of the two-component BEC dynamics with a small difference of the inter and intra-nonlinear constants was first indicated in Ref. [16] for the case of polarization solitons. The appearance of the cot θ -function in the first term in the braces shows that the condition δg/g 1 should be complemented by another condition: θ should not be too close to zero or π so that the right-hand side of (10) remains small. Thus, in addition to (2), we assume also that This condition implies that the densities ρ ↑ or ρ ↓ are not too close to ρ 0 or 0, cf. Eqs. (18). If the conditions (2) and (11) are fulfilled, then the density and polarization dynamics are decoupled and we can study the polarization dynamics separately assuming that ρ(x, t) = ρ 0 = const and disregarding the second equation in the system (5). This approximation greatly simplifies the other equations. The first one reduces to If we choose to work in a reference frame in which there is no flux of the total density, Eq. (12) simplifies to and Φ ζ can then be excluded from the remaining two equations. This yields the system This closed system of nonlinear equations shows that, for time scales of order T p and length scales of order ξ p , the polarization degree of freedom decouples from the density degree of freedom, even in the nonlinear regime. All dimensional parameters have been scaled out from Eqs. (14) which thus correspond to a universal behavior of polarization waves. Note here that the relevant characteristic time (8) and length (7) have been identified for equal densities of both components (and will keep the same value throughout the paper), but the validity of the system (14) does not rely on this assumption: it describes the polarization dynamics in the limit (11), for a system verifying (2). In this case, we see from Eq. (10), that the ratio of the amplitude of density waves with respect to the one of polarization waves is roughly of order δg/g. The system (14) can be derived from the Hamilton principle of extremal action [16] for a Lagrangian Λ = dζ with a Lagrangian density From this expression we can write the (correctly dimensioned) energy of the system under the form where is the energy density corresponding to the Lagrangian (15). This expression coincides with the energy of ferromagnetic bodies in dissipationless Landau-Lifshitz theory [17] with account of dispersion and uniaxial easy-plane anisotropy.
The system (14) can be cast into other forms that may be more convenient in some instances. In particular, the angle θ is related to the density of each component by the formulas hence is the variable describing the variations of the relative density. On the other hand, represents the non-dimensional relative velocity of the components. In terms of the two variables (w, v) which have clear physical meanings, the system (14) takes the form This system coincides with the one-dimensional version of the system derived in the recent preprint [18] for hydrodynamic description of magnetization dynamics in ferromagnetic thin films. For subsonic flows with velocities |v| < 1 we can introduce a variable σ such that v = cos σ, (22) and then in terms of (θ , σ)-variables the system of equations of the polarization dynamics reads θ τ + 2 cos θ · cos σ · θ ζ − sin θ · sin σ · σ ζ = 0, The importance of distinguishing subsonic from supersonic flows-an essential assumption for being able to write the relation (22)-can be seen from the following observation: consider a stationary uniform background characterized by a relative density w 0 and a relative velocity v 0 . Linear perturbations of the form obey the following dispersion relation: By definition we always have |w 0 | ≤ 1, however v 0 can have any value, and for |v 0 | > 1 the frequency ω is complex for small enough wavevectors k. This implies a long wavelength instability of supersonic relative motions of two-component superfluids, more precisely for a background relative velocity v ↓ − v ↑ larger than 2c p . This mechanism of instability has been first theoretically studied in Ref. [23], and the regime (2) we consider here corresponds to what is denoted as "strong coupling" in this reference. We note here for future use that, for subsonic flows with w 0 = cos θ 0 and v 0 = cos σ 0 , the dispersion relation (24) can be written as ω = 2 cos σ 0 cos θ 0 ± sin 2 σ 0 sin 2 θ 0 + k 2 k .
The long wave length behavior of the dispersion relations (24) and (25) is linear and corresponds to a velocity of sound in the laboratory frame For a uniform system in which both components have equal densities (w 0 = 0) and no relative velocity (v 0 = 0) one gets c = ±1, i.e., going back to dimensional quantities, the speed of the polarization sound is c p as expected.

Traveling waves and solitons of polarization
In this section we consider traveling wave for which the physical variables θ and v depend on ξ = ζ − V τ only, V being the phase velocity of the wave. In the framework of the system (14) this corresponds in making the ansatz that the velocity potential φ(ζ, τ) and θ (ζ, τ) can be represented as φ(ζ, τ) = qζ +φ(ξ), θ = θ (ξ).
Substitution into the first equation of the system (14), multiplication by sin θ and integration give at onceφ where B is an integration constant. Substituting this expression into the second equation of the system (14) gives after simple transformations the equation Multiplication by θ ξ and integration yield the final equation for the variable w = cos θ : where C is an integration constant. The four parameters V, q, B, C can be expressed in terms of the four zeroes w 1 ≤ w 2 ≤ w 3 ≤ w 4 of the polynomial where the s i 's are standard symmetric functions of the zeroes w i 1 . In particular, we obtain and The solution of Eq. (30) can be expressed in terms of Jacobi elliptic functions and, without going into well-known details (see, e.g., [24]), we shall present here the final results. The variable w can oscillate between two zeroes of the polynomial Q(w) where Q(w) ≤ 0 provided these two zeroes are located in the interval [−1, 1]. There are two possibilities, labeled as (A) and (B) below.
(A) In the first case the periodic solution corresponds to oscillations of w in the interval The solution of Eq. (30) can be written as To simplify the notations, we put in (35) (and in all subsequent similar equations) the integration constant ξ 0 equal to zero. A standard calculation yields where cn and sn being Jacobi elliptic functions [25]. The wavelength is given by where K(m) is the complete elliptic integral of the first kind [25]. In the limit w 3 → w 2 (m → 1) the wavelength tends to infinity and the solution (36) transforms to a soliton This is a "dark soliton" for the variable w = cos θ . The limit m → 0 can be reached in two ways.
This is a small-amplitude limit describing propagation of a harmonic wave.
(ii) If w 4 = w 3 but w 1 = w 2 , then we get a nonlinear wave represented in terms of trigonometric functions: If we take the limit w 2 − w 1 w 3 − w 1 in this solution, then we return to the small-amplitude limit (40) with w 4 = w 3 . On the other hand, if we take here the limit w 2 → w 3 = w 4 , then the trigonometric functions in (41) have a small argument and can be approximated by the first terms of their series expansions. This yields a solution which we denote as an "algebraic soliton": (B) In the second case the variable w oscillates in the interval so that instead of (35) we get Again, a standard calculation yields with the same definitions for W , m, and L as in Eqs. (37) and (38). In the soliton limit This is a "bright soliton" for the variable w = cos θ . Again, the limit m → 0 can be reached in two ways.
(i) If w 4 → w 3 , then we obtain a small-amplitude harmonic wave This is a small-amplitude limit describing a harmonic wave.
(ii) If w 2 = w 1 , then we obtain another nonlinear trigonometric solution, If we assume in this solution w 4 − w 3 w 4 − w 1 , then we reproduce the small-amplitude limit (47) with w 2 = w 1 . On the other hand, in the limit w 3 → w 2 = w 1 we obtain the algebraic soliton solution: This ends the general presentation of the different solutions of Eq. (30).
It is now interesting to discuss in more detail the soliton solutions which play a special role in the description of dispersive shock waves (Sec. 5.2). The bright soliton solution (46) corresponds to an increased number of particles in the "up" component: where ρ is the background density of the up component and ρ sol ↑ (ζ, τ) = ρ 0 (1 + w)/2, w(ξ) being given by (46). One gets The soliton is characterized by the three zeros w 1 , w 2 (= w 3 ) and w 4 which relate to the physical variables w 0 (relative background density of the components), V (velocity of the soliton) and v 0 (relative background velocity of the components) through The energy of the soliton is the difference between the energy (16) of the system in the presence and in the absence of the soliton. It reads u(ξ) being here the energy density (17) computed for the distribution (46). It is shown in Appendix A that The soliton solution found in [16] is reproduced from Eqs. (39) and (46) if we consider the situation where the two components have equal background densities (w 0 = 0), and no relative velocity (v 0 = 0). In this case, one gets from Eq. (52) which agrees with formula (32). As a result we obtain and Eqs. (18) give the corresponding densities of each component. From (51) and (54), one sees that this soliton corresponds to an increase of the number particles of the up component ∆N ↑ = π 2 ρ 0 ξ p and to an energy E sol = ħ hρ 0 c p 1 − V 2 , in agreement with the findings of Ref. [16]. Note however that the existence of polarization solitons of the form (39) and (46) is not restricted to the condition of equal background densities ρ ↑0 = ρ ↓0 considered in Ref. [16].
Our approach made it possible to identify new algebraic solitons (42) and (49) with unique properties which we now briefly discuss. The algebraic soliton (49) can be obtained as the limit w 2 (= w 3 ) → w 1 of (46). It corresponds to an increased number of "up" particles ∆N ↑ = πρ 0 ξ p . At variance with the case of dark/bright solitons, once the background parameters w 0 and v 0 are fixed, the velocity V of an algebraic soliton is not free. One finds that it is fixed to be exactly the sound velocity (26). For an algebraic soliton, one has w 2 → w 1 and thus the energy (54) of such a soliton is zero, as can be checked directly from (49) and (53).
Also note that the dark/bright solitons (39) and (46) are of a quite different nature than the one identified by Busch and Anglin in Ref. [26] and observed in Ref. [27]. It can be shown that if one considers the limit of a stationnary soliton of type (46) with no pedestal (w 0 → −1), then one does not reach the limit of the dark-bright solitons of Ref. [26], but instead one obtains an algebraic soliton of the form ρ ↑ (ζ, τ) = ρ 0 (1 + ζ 2 ) −1 .

Dispersionless hydrodynamics and Riemann equations
If the velocity and density distributions v and w are smooth enough, that is, if they experience little change over one polarization healing length (7), then we can neglect the dispersion effects described by the last terms in the second equations of the systems (21) and (23) 2 This corresponds to the so-called dispersionless approximation. We shall present the corresponding equations in two forms-for the variables (w, v), and for the variables (θ , σ), These are equations of hydrodynamic type which can be studied by means of well documented methods.
First of all, we find at once from the system (58) that the variables where or in terms of the variables (v, w) 2 In this regime, the dispersion relation (24) can be approximated by a straight line of slope c [c being the speed of sound, as given by (26)], which is legitimate when k 1, i.e., for wave lengths large compared to ξ p .
The characteristic velocities V 1 and V 2 are the velocities of propagation of small disturbances along a non-uniform background (θ , σ) or (w, v), correspondingly. In the case of a uniform background w = w 0 = cos θ 0 , v = v 0 = cos σ 0 they coincide with the sound velocities (26). The variables r 1,2 are called Riemann invariants, and Eqs. (60) are the hydrodynamic equations written in the Riemann invariant form (see, e.g., Ref. [28]). They have the familiar form of equations of compressible gas dynamics written in terms of the Riemann invariants, however the relationships between the Riemann invariants and the physical variables are more complicated here than for a gaseous system. Once r 1 and r 2 have been found, the physical variables w, v are given by At this point, we have reduced the polarization hydrodynamic equations to the symmetric Riemann form (60). We shall now study a special class of solutions of these equations.

Simple wave solutions
In the framework of the hydrodynamic approximation a special role is played by the so-called simple wave solutions that are characterized by the fact that one of the Riemann invariants (59) is constant along the solution, so that the system (60) reduces to a single equation of the Hopf type. For example, let r 2 = r 0 2 = const; then we get the equation for the variable r 1 . This equation admits the well-known solution where f (r 1 ) is an arbitrary function. Equation (65) determines the dependence of r 1 on ζ and τ in an implicit form. The function f (r 1 ) can be thought of as the inverse function of the initial distribution of r 1 at the moment τ = 0, i.e., f −1 (ζ) = r 1 (ζ, τ = 0). The simple wave solution with constant Riemann invariant r 1 = r 0 1 = const can be easily written in a similar form. The importance of the simple wave solutions is related to the fact that, generally speaking, a hydrodynamic solution of a typical problem consists of different functions defined on several regions in the (ζ, τ)-plane separated by lines of discontinuity of the fields (here θ and σ). Along the so-called weak discontinuities one has discontinuities of the derivatives while the functions remain continuous. In particular, if the fluid flow has a boundary with adjacent quiescent fluid, then this boundary is a weak discontinuity and the neighboring flow is described by a simple wave solution (see, e.g., [28]).
A special role is played by self-similar solutions, for which r 1,2 depend on the self-similar variable z = ζ/τ only. In particular, such solutions appear in problems where the initial distributions do not contain parameters with dimension of a length, e.g., in the case of initial discontinuities with abrupt jumps of the variables w and/or v (θ and/or σ). The jump occurs at some coordinate that can be taken as the origin of the ζ-coordinate frame. In this case, r 1,2 = r 1,2 (z) and the hydrodynamic equations (60) take the form Their solutions are evidently r 2 = r 0 2 = const, and V 1 (r 1 , r 0 2 ) = z, or r 1 = r 0 1 = const, and V 2 (r 0 1 , r 2 ) = z.
where the values of the constants (r 0 2 or r 0 1 and n ∈ ) and the signs are to be determined from the boundary conditions. Let us consider here such solutions in the case where a dispersionless polarization flow is neighboring a condensate at rest. We shall first consider a self-similar simple wave matching at its right side a quiescent condensate (i.e., with σ = π/2) where θ = θ R = const. It is easy to see from simple considerations [28] that its right edge, being a weak discontinuity, must propagate to the right with the sound velocity c = sin θ R [cf., (26)]; that is, this self-similar flow has to satisfy the boundary condition θ = θ R at z = z R = sin θ R . Simple inspection shows that this is achieved by the first of solutions (68) (where r 2 = σ + θ = π/2 + θ R = const) with a lower sign and n = 0. Hence, owing to the relation arccos x = π/2 − arcsin x, we obtain and, consequently, by virtue of constancy of r 2 = r 0 2 = π 2 + θ R , It is usually supposed that θ takes values in the interval 0 ≤ θ ≤ π, however any interval of same length is suitable for the description of the physical variable w = cos θ . We shall use here the equivalent interval which is more suitable for the solution (69). The solution (69) does not cover all the interval (71) over which one has The resulting plot is displayed in Fig. 1 for a value of θ R chosen in the interval 0 < θ R < π.
The left edge of this wave must have a boundary either with one of the general solutions of equations (60), or with another simple wave with constant values of σ and θ (that is, with a plateau in the density distribution). For future applications we shall confine ourselves to the second possibility and demand that the left edge of the solution corresponds to θ = θ 0 and, consequently, to σ = σ 0 = π/2 + θ R − θ 0 , since r 2 is constant across our simple wave. Here we have to distinguish two main typical situations denoted as (a) and (b) below.
Case (a): If then the constant left flow characterized by σ 0 and θ 0 is connected with the quiescent condensate at the right by a rarefaction wave shown in Fig. 2 (region z 0R < z < z R of this figure) whose left edge propagates with velocity The corresponding distributions of the density ρ ↑ and the flow velocity v = cos σ = sin(θ − θ R ) are shown in Figs. 3 and 4, respectively.
Case (b): We will see that in this case there is an interval on the z-axis where the formal solution of the hydrodynamic equations becomes three-valued. Although such a solution does not have a direct physical meaning, it provides important relations remaining correct after replacement of the nonphysical multi-valued parts of the flow by a dispersive shock wave. To be definite, we illustrate such a situation in Fig. 5 which is drawn in the subcase we denote as (b1) in which   Figure 6: Distribution of ρ ↑ (z) in the simple wave solution; case (b1). Here z R = sin θ R , z 0R = 3 2 sin(2θ 0 − θ R ) − 1 2 sin θ R .
In this case, in the region of the simple wave, θ (z) is given by the single-valued solution (69), but the matching with the left and right boundaries can only be performed at the price of overlapping the region of validity of the single wave solution with the ones of the plateau at the boundary. This corresponds to an overall multi-valued solution, as shown in Fig. 5. The corresponding plot of the density is shown in Fig. 6 and a similar graph can be plotted for the flow velocity v(z).
In the subcase we denote as (b2) for which the simple wave solution obtained from (72) already corresponds to a multi-valued θ (z) and the graphs of the formal hydrodynamic solutions can be easily plotted.
Let us now turn to a self-similar simple wave propagating to the left into a quiescent condensate with σ = π/2, θ = θ L = const. This problem is obviously symmetric to the one just studied: the left edge of the wave propagates here to the left with the sound velocity c = − sin θ L that is, we have to satisfy the boundary condition θ = θ L at z = z L = − sin θ L . This time we have to consider the second of solutions (68) (where r 1 = σ − θ = π/2 − θ L = const) with an upper sign and n = 0. Hence, we obtain and, consequently, It is clear that the plots for this case can be obtained from the previous ones by the change z → −z replacing the notation θ R → θ L , etc. Therefore we shall illustrate such a situation only by the plot of θ (z) which is displayed in Fig. 7. Thus, we have obtained simple wave solutions which match on one boundary with a quiescent uniform condensate, and on the other with a flow with constant density and velocity-the "plateau solution".
Two important typical situations have been identified in this section. First, in some cases, the plateau solution can be connected to a simple wave solution joining a quiescent condensate on its other boundary. This is the situation illustrated in Figs. 3 and 4. For such flows the dispersionless hydrodynamic approach is indeed legitimate, and it is just expected that a more precise treatment of the weak discontinuities should exhibit a small amount of linear radiation (on both sides of the simple wave). Such flows are called rarefaction waves. Second, in some instances, the solution of the dispersionless hydrodynamic approach is multi-valued in some regions of space, cf. Fig. 6. In these regions, the physical flow is expected to be a dispersive shock wave, as commonly encountered in similar situations. In the next section we shall consider a configuration where these two possibilities are realized.

Evolution of a step-like discontinuity
As a typical application of the theory, let us consider an initial step-like distribution of polarization θ (ζ, τ = 0) = θ L , when ζ < 0 , and we assume here that the left and right asymptotic regions are both initially at rest, σ(ζ, τ = 0) = σ L = π 2 , when ζ < 0 , σ R = π 2 , when ζ > 0 . (81) We shall consider this problem in the framework of the polarization dynamics governed by Eqs. (14), (21) or (23). We shall begin with the dispersionless hydrodynamic approximation corresponding to Eqs. (57) or (58) that can be written in the Riemann invariant form (60).

Hydrodynamic approximation
The step-like discontinuity evolves into a wave whose edges propagate into quiescent regions located at ζ → ±∞. If such an edge is represented by a weak discontinuity, then the adjacent flow is described by a simple wave solution. The step-like initial distribution (80) does not include any parameter having the dimension of a length and, consequently, the solution has to depend only on the self-similar variable z = ζ/τ (and of course also, parametrically, on θ L and θ R ). One cannot find a single simple wave joining its right and left boundaries with asymptotic regions corresponding to the initial conditions (80) and (81). Instead, the initial discontinuity evolves, for τ > 0, into a more complex structure: an expanding self-similar wave consisting of two simple waves separated by a plateau characterized by the constant parameters θ 0 and σ 0 . One edge of each simple wave has a boundary with a condensates whose parameters are given by one (the left or the right) of the boundary conditions (80) and (81), the other edge matching the plateau distribution. As was discussed in the preceding section, along the simple wave solution [matching with the left asymptotic region σ = π/2, θ = θ L ] we have r 1 = σ − θ = π/2 − θ L = σ 0 − θ 0 , and along the other simple wave solution [matching with the right asymptotic region σ = π/2, θ = θ R ] we have r 2 = σ + θ = π/2 + θ R = σ 0 + θ 0 . These two conditions determine the parameters of the flow on the plateau: Combining with the simple wave solutions (whose characteristics are discussed in the previous section), we find the full solution of the problem -determined within the dispersionless approach -under the form where The edge at ζ = −z L · τ propagates to the left at velocity − sin θ L which is the sound velocity in the left condensate. The edge at ζ = z R · τ propagates to the right with velocity sin θ R which is the sound velocity in the right condensate [cf. (26)], and the plateau is located between the edges z 0L · τ ≤ ζ ≤ z 0R · τ.
Thus, for given values of the densities at both sides of the initial discontinuity (i.e. for given values of θ L and θ R ) one can calculate the parameters θ 0 , σ 0 defining the plateau distribution from (82) and determining the "left" and "right" simple wave solutions joining the quiescent condensates with the plateau. One of these simple waves represents a rarefaction wave and the other one describes a formal non-physical multi-valued solution. This means that the hydrodynamic approximation fails in the region where the flow is multi-valued and we have there to take into account the dispersion effects neglected in the long wavelength hydrodynamic theory. As a result of dispersion effects, the multi-valued region is replaced by a dispersive shock wave which is an oscillatory nonlinear wave structure. Such a situation is illustrated in Fig. 8. There the orange line describes the hydrodynamic approximation (83), for  (14) for an initial profile given by v(ζ, τ = 0) = 0 , and θ (ζ, with θ L = 0.2 π , θ R = 0.4 π , and ζ 0 = 1 . The value of ζ 0 is not negligibly small, and the argument previously invoked for justifying the self-similar nature of the flow does not hold for all times. Instead, the structure of the flowwith a well defined plateau joined to both asymptotic regions by specific structures -does not appear instantaneously, but takes a finite amount of time to get formed. As a result, the flow can be considered as self-similar only for times large compared with this set-up time, which we numerically evaluate to be of order of τ setup 8 in the case of the initial conditions specified by (85) and (86).
It is clearly seen from Fig. 8 that both the right rarefaction wave and the plateau region are very well described by the hydrodynamic theory, the dispersion effects leading only to small oscillations in vicinity of the weak discontinuities located at the interface between these two regions. On the contrary, the region of large amplitude oscillations on the left side of the wave pattern is completely beyond reach of the dispersionless approach and in the next subsection we shall use a theory able to describe such dispersive shock wave (DSW) structures with account of dispersion effects.

Whitham modulation theory and Gurevich-Pitaevskii problem
As seen in Fig. 8, the numerical solution suggests that the dispersive shock wave can be seen as a nonlinear periodic solution of the polarization equations -such as those studied in Section 3 -which is however modulated, as shown by the fact that the amplitude of the oscillations is not constant. This modulation is gentle, in the sense that the parameters (amplitude, velocity, period, etc.) of the wave change little over one wavelength and one period of oscillation. This means that we can apply the Whitham averaging method for the description of this structure. In his original paper [29], Whitham assumed that the evolution of slowly modulated nonlinear waves can be described by equations obtained by averaging the densities and fluxes of the conservation laws over the rapid oscillations of the wave. He derived these averaged equations for several nonlinear wave equations, in particular, for the case of cnoidal wave solutions of the celebrated Korteweg-de Vries (KdV) equation, and-what was most remarkable from a mathematical point of view-he succeeded in transforming these equations into a diagonal Riemann form analogous to equations (60) obtained in the dispersionless approximation of hydrodynamic flows. As it became clear later, this success was related to the specific mathematical properties-complete integrability-of the KdV equation.
For the case we are interested in, a most important application of the Whitham theory was suggested by Gurevich and Pitaevskii [30]. In their approach it was assumed that the expanding DSW which develops after wave breaking can be described by the nonlinear periodic solution of the wave equation provided the parameters of this solution change slowly with time and space coordinate. They illustrated the method by applying it to the evolution of an initial step-like discontinuity and to the formation of a DSW after the wave breaking moment for the KdV wave dynamics.
Since the publications of the work of Whitham and Gurevich and Pitaevskii, the Whitham theory has been considerably developed in different directions and has found many applications in nonlinear physics. In particular, it was shown that many problems can be reduced to the consideration of the evolution of an initial step-like discontinuity. It was therefore of great importance to discover [31] that, for this specific step-like problem, the main characteristics of DSWs can be obtained by a simple method applicable to both completely integrable and non-integrable nonlinear wave equations. In our case the polarization wave dynamics is governed by the 1D version of the dissipationless Landau-Lifshitz equation which is completely integrable (see, e.g., [32]). However, the Whitham theory is not developed well enough for this equation and therefore El's method [31] seems the most appropriate for the description of the DSW observed in Fig. 8.
We thus assume that, instead of the multi-valued solutions found in the dispersionless approximation in the preceding subsection, a DSW is generated that joins the neighboring quiescent condensate at the left side of the wave structure with the plateau region. For definiteness, and in accordance with the example shown in Fig. 8, we consider the case where the Riemann invariant r 1 = σ − θ is constant across the multi-valued region. As was assumed by Gurevich and Meshcherkin [33] -and confirmed in many particular cases -one of the Riemann invariants preserves its value even after replacement of the multi-valued solution by the oscillatory DSW: in a sense, an equality of the type r 1 | − = r 1 | + replaces in the case of DSWs the well-known Rankine-Hugoniot relation of the theory of viscous shocks. It is then natural to assume that this relation is preserved by the Whitham averaging method, which yields an appropriate interpolation between the two edges of the DSW.
As we know, at the small-amplitude edge the DSW can be approximated by a modulated linear wave (47), however now propagating along a non-uniform background corresponding to the simple wave solution with r 1 = σ − θ = π/2 − θ L = const, where we have used the values of the parameters at the left edge that matches with the left boundary conditions. With help of this relation we can write σ = π/2 − (θ L − θ ) in the dispersion relation (25), leading to In (87) we have chosen the minus sign in front of the square root of (25) because we consider wave propagating to the left with respect to the background condensate. Equation (87) is the dispersion of linear waves propagating along a non-uniform θ -distribution. During the smooth evolution of the oscillatory structure the local "number of waves" is preserved [34] which is expressed by the equation Following El [31], we make a simple-wave type of assumption: in the DSW the wave number k is a function of θ only, k = k(θ ). Then, with account of (87), the law (88) of conservation of number of waves can be written under the form On the other hand, substitution of σ = π/2 + θ − θ L into the first of equations (23) yields Imposing consistency of (89) and (90) considered as equations for θ , we get This is El's equation that can be extrapolated into the large amplitude nonlinear region by imposing the condition that the wavelength tends to infinity at the soliton edge, that is Introducing the function makes it possible to cast equation (91) into the form whose solution-with account of the boundary condition (92)-reads This yields Consequently, the left edge of the DSW propagates with the group velocity evaluated at k(θ L ): At the soliton edge of the DSW, we use the "soliton dispersion law" [31] Ω relating the velocity V = Ω/κ of the soliton with the inverse width κ that describes the exponential profile w ∼ = w 3 + 1 2 (w 4 − w 3 ) exp{−κ|ζ + V τ|} of the soliton far away from its center (in the regime |ζ| → ∞). The relation (98) follows from the remark that the soliton's tail propagates with the same velocity as the soliton itself and therefore the soliton's velocity can be found from the asymptotic behavior of its profile, see, e.g., [35,36]. Again following El, we assume that along the shock κ = κ(θ ). Then the following equation can be derived (see [31]) for this function: If we extrapolate the solution of (99) to the small amplitude region where κ tends to zero, we obtain the boundary condition κ(θ L ) = 0.
Similarly to what has been done for the leading edge of the DSW [Eq. (91)], it is convenient for solving Eq. (99) to introduce the auxiliary functioñ Inserting (101) into (99) and taking into account the boundary condition (100) one obtains Then, at the soliton edge,α is equal tõ and, consequently, this edge propagates with velocity The comparison of the analytic predictions (97) and (103) for the velocities of the edges of the dispersive shock wave with our numerical simulations is easily done for the well defined soliton edge, because, indeed, a leading soliton is easily identified at this edge of the numerically determined DSW. The velocity of this soliton tends for large time to the theoretical value, as illustrated in Fig. 9. In this figure, the numerical result for the velocity V (τ) of the soliton at the interface between the DSW and the plateau region is fitted with the empirical formula As one can see in Fig. 8, it is difficult from the numerical solution to unambiguously locate the dispersive edge of the shock. Hence, at variance with the situation for the soliton edge, the velocity of the dispersive edge cannot be precisely extracted from the numerical simulation. However, one can reasonably argue that the value v gr = −1.54 obtained from the theoretical formula (97) for the initial datas (86) matches quite well with the numerical results (cf. Fig. 8).

Discussion
In this section we discuss the accuracy of the polarization description of the dynamics of a twocomponent BEC [Eqs. (14)] and also the relevance of our approach to experimental studies.
A first question can be asked: in which extend does the assumption of decoupled dynamics apply? In other words, how small should δg/g be in order for the approach followed in the present work to apply? A simple way for answering this question is to compare the results obtained from (14) with the ones obtained from the numerical solution of the full Gross-Pitaevskii system (1). This is done in Fig. 10 which displays the evolution of an initial profile of type (85). As one can see from this plot, the agreement is reasonable already for δg/g = 0.2 and becomes quite good for δg/g = 0.05. The lower part of the Figure shows that the assumption of constant total density is verified with an accuracy of order of 0.5% for δg/g = 0.05. We note that the largest departure of the total density from a constant occurs when ρ ↑ /ρ 0 is close to unity, i.e., when θ is close to 0, as anticipated in Eq. (11). Note also that the spatial and time scales (ξ p and T p ) are quite relevant: the Gross-Pitaevskii system is solved for quite different values of these characteristic scales (the value of ξ p is multiplied by a factor 2 and the one of T p by a factor 4 when one goes from δg/g = 0.2 to δg/g = 0.05), but after the same time expressed in units of T p (24 T p in the case of Fig. 10  Another question naturally arises: since Bose-Einstein condensation of ultra-cold atomic vapors is always realized in trapped systems, it is important to evaluate the experimental relevance of the infinitely extended configuration studied in the present work. One can first state that the theory has a physical meaning as long as its characteristic length ξ p (7) is much less that the size X of spatial overlap of the two components which can be estimated in the framework of the Thomas-Fermi approximation presented in Appendix B : where ω is the longitudinal trapping angular frequency and ρ 0 N /X , N being the total number of atoms. The condition (104) combined with (2) reads where ξ = ħ h/ 2mρ 0 g is the healing length (ξ p = ξ g/δg ). The first inequality of (105) can be also rewritten as ω ξ c p or ξ that is the polarization sound velocity must be much greater than the healing length divided by the period of oscillations of atoms in the trap, or, in other words, the polarization wave passes the healing length in a time much less that the period of oscillations in the trap.
It is also worthwhile to address another point: it is known [37][38][39] that, in the presence of a trapping potential, the condition of uniform miscibility (which, in our notations, reads δg > 0) is not sufficient to ensure a good spatial overlap of the two components. This point is discussed in Appendix B where it is shown that, close to the mixing-demixing transition, the trapping potential induces a kind of phase separation if the lower of the intra-species nonlinear constants (say g ↓↓ ) is smaller than the inter-species constant g ↑↓ , although the criterion of uniform miscibility g ↑↓ < g ↑↑ g ↓↓ is (weakly) fulfilled.
This phenomenon could explain why, in Ref. [40], a kind of phase separation is observed in the mixture of the two hyperfine states | ↓〉 = |F = 1, m F = −1〉 and | ↑〉 = |F = 1, m F = 0〉 of 87 Rb in spite of fulfilment of the uniform mixing condition. For this system (a ↑↑ , a ↓↓ , a ↑↓ ) = (100.86 a 0 , 100.4 a 0 , 100.41 a 0 ), where a 0 is the Bohr radius. Thus a ↓↓ < a ↑↓ which implies mixing of the components in a uniform case; but non-uniformity caused by the trap potential induces phase separation. Instead, for the mixture of the two hyperfine states | ↑ 〉 = |F = 1, m F = −1〉 and | ↓〉 = |F = 2, m F = −2〉 of 87 Rb one has (a ↑↑ , a ↓↓ , a ↑↓ ) = (100.4 a 0 , 98.98 a 0 , 98.98 a 0 ), that is the criterion on miscibility is also fulfilled, but here a ↓↓ = a ↑↓ and the authors observe a large region of overlap of the two components.
Finally, concerning the comparison of our results with the ones presented in Ref. [20], it is worth noticing that if θ L → 0, that is ρ ↑L → 1, then the left edge group velocity (97) tends to the value v gr = −2 sin θ R which coincides with the limiting value of velocity (74) of the left edge of the rarefaction wave z 0R corresponding to θ 0 = 0. This means that the DSW pattern is represented by small amplitude oscillations around the extrapolation of the rarefaction wave to the region with θ L → 0, ρ ↑L → 1. As a result, the pattern looks like the rarefaction wave connecting two regions of quiescent condensates with different values of θ : θ L = 0 and θ R = 0. This apparently agrees with the numerical simulations of the so-called subcritical regime discussed in [20] where only the rarefaction wave was observed for small enough values of the relative velocity and ρ ↑L = 1.

Conclusion
In vicinity of the mixing/demixing transition, in the limit (2) first identified in Ref. [16], the polarization dynamics decouples from density waves and is described by the universal equations (14). In this paper we have identified new specific polarization structures associated with these equations in the case of a one dimensional system: algebraic solitons, simple waves, dispersive shock waves, etc. But more remains to be done. For instance, the non-monotonous behavior of the Riemann velocities (cf. section 4.2) is typically associated to a rich variety of different types of shocks [24] which remain to be investigated in the case at hand; in particular for situations with large jumps of the parameter θ , when DSWs consisting of combined cnoidal and trigonometric parts are expected. The precise behavior of algebraic solitons in several instances, and a reliable procedure for their physical implementation would also be of great interest. The configuration described by the initial distributions (80) and (81) is too schematic for being able to describe the experiments presented in [20] where regions with different density ratios are colliding with finite initial relative velocities. One should thus consider the case where σ L and σ R are not both equal to π/2, and where the plateau formed after the collision is modulationnally unstable. Finally, the approach developed in this paper can be generalized to include Rabi coupling between the components (see, e.g., [41]) and also to two-or three-dimensional situations [42]. In particular, formation of oblique polarization solitons by the flow of the binary condensate past a polarized obstacle (see, e.g., [43]) can be considered in the framework of the present method. Works in these directions are in progress.

A Computation of the energy of a soliton
We briefly present here the computation leading to the result (54) for the energy of the soliton. From (27) and (28) The integrand being symmetric -since w(ζ) is-one can thus restrict the range of integration to the domain (−∞, 0] over which one can write dζ = +d w/ −Q(w). Using the fact that one can express B, v 0 and V as functions of w 1 , w 2 and w 4 [cf. Eq. (52)], it is then possible to re-write (107) under the form which yields the result (54).

B Effective demixing in a 1D trap
In this appendix we present 1D computations in the framework of the Thomas-Fermi description of the system (1) in the presence of a trapping potential [44]. It is known [45] that the Thomas-Fermi approximation cannot quantitatively describe all the possible configurations encountered the mixture of two BECs, but it will permit to identify specific situations which will then have to be confirmed by a full numerical solution.
We consider here N ↑ and N ↓ atoms of each component placed in a harmonic potential of longitudinal angular frequency ω much smaller than the radial trapping angular frequency ω ⊥ . In the so called "1D mean field regime" [46], the system can be described by the effective 1D Gross-Pitaevskii equation (1) with g ↑↑ = 2ħ hω ⊥ a ↑↑ [47] where a ↑↑ is the 3D intra-species s-wave scattering length of the "up" component (an similar expressions for g ↓↓ and g ↑↓ ). In the situation we are interested in where N ↑ ∼ N ↓ and a ↑↑ ∼ a ↑↓ ∼ a ↓↓ , the 1D mean field regime holds when N ↑ (ω /ω ⊥ )(a ↑↑ /a ⊥ ) 1, where a ⊥ = ħ h/mω ⊥ is the radial harmonic oscillator length.
We chose the parameters so that the mean field condition of miscibility a ↑↑ a ↓↓ > a ↑↓ > 0 is always fulfilled, and in the following we denote as A the parameter having the dimension of length defined by A 2 = a ↑↑ a ↓↓ − a 2 ↑↓ > 0 .
We define the non-dimensional position X = x/a , where a = ħ h/mω is the longitudinal harmonic oscillator length, and the non-dimensional densities n ↑,↓ such that n ↑,↓ (X )d X = situation would be expected in the situation a ↓↓ a ↓↑ a ↑↑ . The point is here that the same effect is observed for a system verifying the miscibility condition (109) provided one remains close to immiscibility and that a ↓↓ a ↓↑ . The parameter governing the expulsion of the up component from the center of the trap is the non-dimensional curvature of its density at X = 0. From (113) this parameter is equal to In the cases presented in Fig. 11 the value of this parameter changes from −1.3 (in the left plot of the figure) to +4.2 (right plot) just by changing a ↓↓ by 2%.