{'fraction_non_alphanumeric': 0.03402215242188791, 'fraction_numerical': 0.018548520416597786, 'mean_word_length': 4.532467532467533, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spin modulation instabilities and phase separation dynamics in trapped two-component Bose condensates

In the study of trapped two-component Bose gases, a widely used dynamical protocol is to start from the ground state of a one-component condensate and then switch half the atoms into another hyperfine state. The slightly different intra-component and inter-component interactions can then lead to highly nontrivial dynamics. We study and classify the possible subsequent dynamics, over a wide variety of parameters spanned by the trap strength and by the inter- to intra-component interaction ratio. A stability analysis suited to the trapped situation provides us with a framework to explain the various types of dynamics in different regimes.

In a number of two-component BEC experiments reported over more than a decade, a standard technique has been to start from the equilibrium state of a single-component BEC, e.g. populating a single hyperfine state of 87 Rb and then using a π/2 pulse to switch half the atoms to a different hyperfine state [1][2][3][4][5][6][7][8][9]. This results in a binary condensate where the two intra-species interactions (g 11 and g 22 ) and one inter-species interaction (g 12 ) are all slightly different from each other, but the starting state is the ground state determined by g 11 alone. Since it has been realized several times in several different laboratory setups, this is a paradigm non-equilibrium initial state for binary condensate dynamics. A thorough and general analysis of the dynamics subsequent to such a π/2 pulse is thus clearly important. In this paper, we present such an analysis, clarifying the combined role of the inter-species interaction (g 12 ) and the strength λ of the trapping potential. We provide a stability analysis mapping out regions of the λ-g 12 parameter space hosting different types of dynamics. Since it is now routine to monitor realtime dynamics in such experiments (e.g. [6]), we also directly analyze the real-time evolution after a π/2 pulse.
It is widely known that the ground state of a uniform two-species BEC is phase separated or miscible depending on whether or not the inter-species repulsion dominates over the selfrepulsions of the two species, i.e. if then the ground state is phase separated [15]. This criterion is also a key ingredient in understanding dynamical features such as pattern dynamics in the density difference between the two species-such 'spin patterns' emerge when the phase separation condition is satisfied. (Since it is common to refer to the components of a two-component Bose gas as 'spin' states, e.g. [29][30][31][32], we refer to the density difference as 'spin density' and patterns in the density difference as 'spin patterns'.) The emergence of spin patterns whenever equation (1) is satisfied can be understood as the onset of a modulation instability [16][17][18], identified by the appearance of an unstable mode in the excitation spectrum around a reference stationary state. For a homogeneous situation, linear stability analysis shows that modulation instability sets in when the condition of equation (1) is satisfied [16][17][18]. The situation is different in the presence of a trapping potential. Phase separation in the ground state, as well as the appearance of modulation instability when starting from a mixed state, now requires larger inter-species repulsion [14,19]. This suggests that the region of parameter space where pattern dynamics occurs also depends on the trap. A trap is almost always present in cold-atom experiments, and it is easy to imagine experiments where the trapping potential is not extremely shallow but varies between tight and shallow limits. ('Tight' and 'shallow' will be specified more precisely in section 2.) It is thus necessary to examine the relevance of equation (1) for trapped binary BECs. To this end, we explore different trap strengths spanning several orders of magnitude, and identify the appropriate extensions of equation (1) for the type of spin dynamics resulting from the π/2 protocol described above.
We focus on the effects of two parameters. Firstly, we study the effects of changing cross-species interaction g 12 , thus generalizing equation (1) for trapped situations. Secondly, we explore the role of the relative strength of the trap with respect to the interactions. Our analysis, performed for a one-dimensional (1D) geometry, sheds light on the situation where g 11 and g 22 are close but unequal: (a) the stability analysis is performed for g 11 = g 22 and their difference serves only to select appropriate instability modes; (b) the simulations are performed with g 22 /g 11 = 1.01.
We restrict ourselves to repulsive interactions (g i j > 0). Attractive interactions (g i j < 0) are known to cause collapse of the condensate for a large enough number of particles, even for a single-component BEC [33].
In section 2, we introduce the formalism and geometry. In section 3, we show the results from a linear stability analysis for a sequence of trap strengths, and identify and analyze relevant modulation instabilities. Through an analysis of unstable modes, we present a classification of the parameter space into dynamically distinct regions, in relation to the prototypical initial state explained above. This may be regarded as a dynamical 'phase diagram'. A remarkable aspect is that the 'phase transition' line most relevant to spin pattern dynamics does not arise from the first modulation instability (studied in [14]). This first instability mode is antisymmetric in space, and as a result is not naturally excited in a symmetric trap with symmetric initial conditions. Complex dynamics (not due to collective modes but rather due to modulation instability) is generated only when the first spatially symmetric mode becomes unstable, which occurs at a higher value of g 12 .
In section 4, we provide a relatively detailed account of the time evolution. For each trap strength λ, for values of g 12 not much larger than g 11 , we observe simple collective modes. Above a threshold value of g 12 , the oscillation amplitude becomes sharply stronger, and at the same time the motion becomes notably aperiodic, signaling that the dynamics is more complex than a combination of a few modes. Dynamical spin patterns start appearing at this stage and become more pronounced as g 12 is increased further. The threshold value at which the dynamics changes sharply corresponds to the second modulation instability line rather than the first, as we demonstrate through a careful choice of parameters in each region of the phase diagram derived from stability analysis.
Some further connections between the stability analysis and dynamical features, relating to the length scale of the generated patterns, appear in section 5. In the concluding section 6, we place our results in context and point out open questions.

Geometry and formalism
The relevant time-resolved experiments have been performed in both quasi-1D geometries (highly elongated traps with strong radial trapping) [6] and in a three-dimensional BEC of cylindrical symmetry with the radial variable playing an analogous role as the 1D coordinate [2,5]. Since the basic phenomena are very similar, we expect the same theoretical framework to describe the essential features of each case. For definiteness, in this work we show results for 1D geometry. We expect the general picture emerging from this work to be qualitatively true also for other geometries exhibiting the same type of spin dynamics.
We will restrict ourselves to the mean field regime. Bosonic systems in elongated traps can of course also be in regimes beyond the applicability of mean field descriptions, e.g. when the particle number is small. In such a case a Lieb-Liniger or Tonks-Girardeau description might be more appropriate. Dynamics in such regimes is beyond the scope of this paper. The mean field regime is generally valid when the particle number is large enough [34].
In the mean field framework at zero temperature, the dynamics is described by two coupled Gross-Pitaevskii equations [35][36][37]: Condensate wave functions ψ 1 (x, t) and ψ 2 (x, t) are normalized to unity, and λ is the strength of the harmonic trap. Factors of particle number and radial trapping frequency are absorbed as appropriate into the effective 1D interaction parameters g i j [6,37,38]. We consider purely non-dissipative dynamics, i.e. we do not attempt to model experimental loss rates with a phenomenological dissipative term as was done in, e.g., [5][6][7]. The equations above are in a dimensionless form, because we measure lengths in units of trap oscillator length and time in units of inverse trapping frequency, for a hypothetical trap of unit strength (λ = 1). The scale for trap strengths is itself fixed by imposing g 11 = 1. With this convention, small values of λ correspond to a BEC in the Thomas-Fermi limit. For comparison, we note that the parameters of the experiment of [6] correspond to λ of the order of 10 −5 in these units. For the purposes of this paper, 'tight' and 'shallow' mean, respectively, 10 −2 and 10 −4 in our units. Note that the trap oscillator strength is λ −1/2 . Of course, one can switch between different units via the transformation: where l is an appropriately chosen scale. The initial state after a π/2 pulse involves both components occupying the ground state of a single-component system of interaction 2g 11 , because the atoms were all in the first hyperfine state before the pulse. We model this initial situation as a two-component BEC with g 11 = g 22 = g 12 . The π/2 pulse may then be regarded as a sudden change (a quantum quench [39]) of the interaction parameters g 22 and g 12 .
We use g 11 = g 22 for the stability analysis of section 3. For the explicit time evolution reported in section 4, we use g 11 and g 22 values close but unequal: g 11 = 1, g 22 = 1.01. This choice of close values is convenient for illustrating the structure of the phase diagram, especially for shallow traps. In rubidium experiments the difference between g 11 and g 22 is somewhat larger (in the common case using 87 Rb hyperfine states |1 = |F = 1, m F = −1 and |2 = |F = 2, m F = 1 ); however, our insights should be relevant to a broad regime of possible experiments. A full exploration of the regime of arbitrary differences (g 11 − g 22 ) remains an open task beyond the scope of the present paper.
Numerical simulations presented in section 4 were performed using a semi-implicit Crank-Nicolson method [40,41].

Stability analysis and dynamical 'phase diagram'
We provide in this section a stability analysis for g 11 = g 22 that maps out the regions of λ-g 12 parameter space which support pattern formation instabilities.
Ideally, one might wish to perform a stability analysis around the initial state. However, in contrast to the homogeneous case [16], we are faced with the situation that the initial state is not a stationary state of the final Hamiltonian. The choice of reference state is therefore a somewhat subtle aspect of the present analysis.
We use as the reference state ψ 0 (x) the lowest-energy spatially symmetric stationary state of the case g 11 = g 22 , with parameter g 12 set to its final value. (For large g 12 , this is not the ground state for these parameters, which is phase-separated.) This reference state has the advantage of looking relatively similar to our actual initial state (two components totally overlapping in space), and of being a stationary state of the Hamiltonian for which we analyze linear stability. Our reference state can be regarded as placing both components in the single-component ground state for interaction g 11 + g 12 . We are not aware of a suitable stationary state even more similar to the actual initial state. We will see that our stability analysis around this reference state will predict remarkably well the main observed time-evolution features described in section 4.
Note that it is not natural to use g 11 = g 22 , because stationary states for such a case typically do not overlap completely in space. Instead, in our approach the difference between g 11 and g 22 will play the important role of selecting certain instability modes over others. For this reason, inferences from the present analysis apply only to small relative differences between g 11 and g 22 .
We linearize equations (2) and (3) around the reference stationary state ψ 0 (x): where µ is the chemical potential corresponding to the reference state. By keeping only terms of the first order in δψ 1 (x, t) and δψ 2 (x, t), we obtain a system of linear equations which can be cast in the form Here M is a matrix differential operator which, upon discretization or upon expansion in a set of orthogonal functions, becomes the so-called stability matrix (e.g. [22,24]). We analyze below the eigenmodes of the stability matrix, which we have obtained by numerically calculating the reference stationary state ψ 0 (x) and expanding in the basis of harmonic trap (non-interacting) eigenstates.
Since we use g 11 = g 22 for the stability analysis, eigenmodes will have well-defined 'species parity', i.e. will all be either even [δψ 1 (x, t) = δψ 2 (x, t)] or odd [δψ 1 (x, t) = −δψ 2 (x, t)] with respect to the interchange of species. Even modes describe in-phase motion of the two components and simply correspond to the excitation spectrum of a single-component BEC with interaction constant g 11 + g 12 . Odd modes are more interesting-they describe out-ofphase motion of two components and are therefore reflected in the spin dynamics. Additionally, due to the spatial inversion symmetry x → −x, the solutions will also have well-defined spatial reflection symmetry, and we can distinguish spatially symmetric and antisymmetric modes. Typical eigenspectra are presented in figure 1. In the case of a tight trap λ = 0.2, we notice two modes whose frequencies are nearly constant. These are even modes encoding singlecomponent or in-phase physics. The lower one is the dipole (Kohn) mode with frequency equal to the trap frequency λ. The second nearly constant mode is the breathing mode, which for elongated traps takes a value close to ω 2 = 3λ 2 . The breathing mode (oscillations of cloud size) is visible in the plots of figure 3 (section 4) as a fast oscillation of the total condensate widths.
The two lowest-lying eigenmodes are odd modes encoding out-of-phase physics. For g 12 1, their frequencies are significantly below the breathing mode, and therefore lead to relatively slow oscillations in the spin density. This will also be visible in the real-time dynamics presented in section 4 (the first two columns of figures 3 and 4). The forms of the corresponding eigenvectors are shown in the lower right of figure 1. The nature of the eigenvectors shows that the motion related to the lowest mode corresponds to the left-right oscillations of the two species, while the next odd mode corresponds to spatially symmetric spin motion. The frequencies of these two modes become imaginary at certain values of g 12 , thus leading to the onset of modulation instabilities. The antisymmetric mode becomes unstable at smaller value of g 12 (g a 12 ≈ 1.6 for λ = 0.2) in comparison with the symmetric mode (g s 12 ≈ 2.4 for λ = 0.2). In a spatially symmetric trap, there is no natural mechanism for exciting the spatially antisymmetric mode. On the other hand, any difference between g 11 and g 22 naturally excites the second (spatially symmetric) mode. Thus, the second mode, occurring at larger g 12 , is the relevant instability for understanding the dynamics observed in experiments and explored numerically in section 4.
We find similar excitation spectra for trap strengths λ spanning several orders of magnitude. The spatially antisymmetric mode becomes unstable before the spatially symmetric mode, and both instabilities get closer to 1 as the trap gets shallower. For example, for λ = 10 −3 (also shown in figure 1) the lowest instability sets in for g a 12 ≈ 1.02, while the next one appears at g s 12 ≈ 1.05. The distinction between two instabilities becomes ever smaller as we go toward a uniform system λ → 0, where the phase-separation condition, equation (1), becomes exact. Nevertheless, even for shallow traps, the issue is not purely academic as the precision in experimental measurement and control of scattering lengths continues to improve [6,42].
In figure 2 (main panel), the results of the stability analyses are combined to present a dynamical 'phase diagram'. The two lines show the two instabilities (g a 12 and g s 12 ) as a function of trap strength λ. For very shallow traps, the two transition lines merge as g s 12 ≈ g a 12 ≈ 1. The lower transition line (g a 12 ) was previously introduced in [14]. However, for a trap and initial state with left-right spatial symmetry, this is not the relevant dynamical transition line, because the first even mode only becomes unstable at some higher g 12 value, given by the g s 12 line. In section 4, we will see that spin pattern dynamics is indeed only generated when the inter-component repulsion g 12 exceeds the second instability line (g 12 > g s 12 ), and that crossing the first instability (g a 12 < g 12 < g s 12 ) is not enough for pattern formation in a spatially symmetric trap.

Dynamical features across the parameter space
In this section, we present and analyze the dynamics obtained from direct numerical simulation of the Gross-Pitaevskii equations (2) and (3), after the system is initially prepared in the ground state of the situation g 11 = g 22 = g 12 = 1. The subsequent dynamics is performed with g 11 = 1, g 22 = 1.01, and several different values of g 12 for each trap strength λ.
It is difficult to show the full richness of pattern dynamics through plots of a few quantities. We choose to show the dynamics through two types of plots (figures 3 and 4). Figure 3 shows the time dependence of the root mean square widths of the two components while figure 4 shows density plots of the density difference (spin density), |ψ 1 (x, t)| 2 − |ψ 2 (x, t)| 2 . In both figures, each row corresponds to a different trap strength (λ), and we approach the shallow trap (Thomas-Fermi) limit going from top to bottom.
For each λ the four values of g 12 (table 1).
one, and the third on each row is at g 12 values just above the second, relevant, instability. The fourth panel on each row is at higher g 12 values. The choice of g 12 values with respect to instability lines is clear in the tighter traps of the top three rows, as also shown by squares in figure 2. For shallow traps (lower rows), the instability lines are too close together and too close to g 12 = 1, so making such choices is not meaningful. In the following, as we compare the features of the different columns, we implicitly exclude the lowest row (smallest λ). This is also indicated by the fact that the schematic instability lines in figures 3 and 4 are not extended to the lowest row.
Broadly speaking, we note that there is only regular (collective-mode) dynamics in the second-column figures (g a 12 < g 12 < g s 12 ) even though an instability is present. There is generally a sharp difference between the second and third figures in each row, indicating that the second instability (g s 12 ) is the relevant one. The fourth panel on each row is at higher g 12 values, showing more rich dynamics. Tighter to shallower traps  While the first two columns look qualitatively the same and show regular oscillatory dynamics, in the third column we observe aperiodic motion of stronger amplitude that we relate to the onset of spin pattern dynamics. The spin dynamics is even more pronounced in the fourth column.
In figure 3, we show the time dependence of the individual widths (w 1 , w 2 ) and also of the total root-mean-square width, w(t) = (w 2 1 (t) + w 2 2 (t))/2. Consistent with our observation that spatially symmetric modes (and not the antisymmetric ones) are naturally excited in the current setup, the dynamics shows signatures of the two most prominent spatially symmetric modes noted in figure 1. The breathing mode is the easiest to note and most ubiquitous-it shows up in almost every parameter choice as oscillations in the total density (in-phase in the two components), with a typical period given by 2π/ √ 3λ ≈ 3.63/λ. This follows from the frequency of this mode being almost constant near √ 3λ.   We also see out-of-phase motion of the two components, associated with the lower spatially symmetric mode in figure 1, which has odd species parity. In the first two columns of figure 3, corresponding to smaller values of g 12 such that this mode has small real frequencies, this is excited as a regular 'spin' mode. For example, at λ = 10 −3 and g 12 = 1.04, we observe an out-ofphase oscillation with the period of approximately ≈ 30, much slower than the breathing mode.
In addition, the oscillation period of the out-of-phase motion is slower in the second than in the first column of each row, corresponding to the decreasing frequency of the mode, as seen in stability analysis (figure 1). Once g 12 becomes large enough that the instability threshold for this mode is crossed, the oscillation amplitudes increase sharply and the width dynamics becomes strongly aperiodic and irregular (the third column of figure 3). This signifies the onset of pattern dynamics, as opposed to the excitation of a regular collective mode around a stable state. Irregularity of the width dynamics at stronger g 12 is even more apparent in the fourth column of figure 3.
It is noteworthy that the spatially antisymmetric modes play no role and do not show up in these dynamical simulations. We see no signature of the Kohn mode. Nor do we see any sharp change associated with the instability of the antisymmetric mode, i.e. there is no sharp difference between the first two columns of figure 3.
In figure 4, we show the dynamics of the 'spin density' |ψ 1 (x, t)| 2 − |ψ 2 (x, t)| 2 . The case of very shallow traps (last row) resembles the data in [17,18]. As in figure 3, the first two columns show regular oscillations, corresponding to collective modes without instability. A sharp change occurs, not across the first instability line (between the first and second columns), but instead across the second instability line (the second and third columns), especially for tighter traps (top three rows) where the comparison with instability lines is meaningful. The sharp change can be noted through the color scales, which is dramatically different between the second and third columns in the upper rows.

Length scales of patterns
In homogeneous stability analysis, the length scale of patterns is inferred from the wavevector (momentum) at which an instability first occurs. Since we perform our stability analysis specifically for trapped systems, we do not have a momentum quantum number. Nevertheless, the eigenvectors of the unstable modes contain information about the form of patterns generated in the dynamics of the trapped system. This is illustrated in figure 5, where the eigenvectors of the lowest unstable spatially symmetric modes are shown for several values of g 12 , together with the spin patterns generated in the non-equilibrium dynamics. There is a close match between the distance between nodes of the eigenvectors (rough analogue of 'wavevector') and the length scales involved in the patterns. As in previous figures, the lengths are shown scaled with the trap oscillator length, λ −1/2 , which plays the role of setting the overall length scale for the cloud.
In figure 4, we see that the patterns contain more spatial structure in shallow traps. The top two rows (tight traps) only show in-out-in type of patterns. This can be understood from the idea that the interactions induce length scales ('healing lengths') in the problem, which are smaller for larger interactions and which set the length scale of spatial structures. For tight traps, the healing length set by the interactions is large or comparable to the cloud size, so that only global dynamical patterns are generated. In such traps, generation of complex patterns with many spatial oscillations would require much higher values of (g 12 /g 11 − 1). For shallow traps, the healing length becomes much smaller than the cloud size; as a result one can have a multitude of dynamical spin structures in the system, of the type seen in experiments and prior simulations [5,6,18]. This heuristic explanation can be made quantitative by counting the number of nodes appearing in the eigenmodes (as in figure 5). Generally, there are as many nodes in the eigenmodes as there are crossings between positive and negative spin densities seen in the dynamical patterns. Accordingly, as we have more nodes at larger g 12 , we get more intricate dynamical patterns.

Conclusions and open problems
In this paper, we have analyzed a widely used dynamical protocol for two-component BECs, which involves starting from the ground state of one component and switching half the atoms to a different component through a π/2 pulse. We have presented a stability analysis suitable to the trapped situation and also presented results from extensive dynamical simulations. Through an analysis of unstable modes, we have presented a classification of the parameter space into a number of dynamically distinct regions, in relation to the prototypical initial state. This may be regarded as a dynamical 'phase diagram'.
In the 'stable' regime of parameter space (no modulation instabilities), our stability analysis explains the observed slow spin oscillations compared to the fast breathing mode oscillations of the total density. We demonstrate that the important 'phase transition' line for spatially symmetric situations relevant to most experiments is not the first instability (studied in [14]), but a second transition line. The first instability is antisymmetric in space, and as a result is not naturally excited in a symmetric trap.
Our stability analysis is performed relative to a stationary state of the situation g 11 = g 22 . The π/2 pulse of the experiments (in the cases where g 11 = g 22 ) can be considered as turning on a nonzero (g 11 − g 22 ), i.e. turning on 'buoyancy' such that one component gains more energy by being in the interior of the trap compared to the other. This helps to select instability modes which are symmetric in space.
Since we have used a stability analysis with g 11 = g 22 to analyze dynamics with g 11 = g 22 , an obvious question is how the ratio g 22 /g 11 affects the regime of applicability of this scheme. We expect that features of this (g 11 = g 22 ) stability analysis are useful for dynamical predictions as long as g 12 /g 11 − 1 is roughly more than g 22 /g 11 − 1. For example, for shallow traps (small λ), the instabilities occur at g 12 /g 11 − 1 values comparable to 0.01, which is why the placement of parameters in the three dynamical regions of the 'phase diagram' (figure 2) is not meaningful for the smallest λ values (the lowest rows of figures 3 and 4).
For the stability analysis, we used a reference stationary state which is of course not the initial state: the initial state is the ground state for g 11 = g 22 = g 12 , while the reference state is the lowest-energy spatially symmetric stationary state corresponding to the final value of g 12 . The instability lines found in this stability analysis would describe even better a situation where the dynamics is triggered by a small quench of g 12 , as opposed to the changes of g 12 that we consider here, which can be relatively large. We have looked at some examples of this type of dynamics and indeed find instabilities matching the stability analysis extremely well. However, although the initial state in the π/2 dynamics is somewhat different from the reference state of our stability analysis, our results show that this stability analysis does provide an excellent overall picture of the dynamics generated by the π/2 protocol.
Our work opens up a number of questions deserving further study. Firstly, we have thoroughly explored the λ-g 12 parameter space, while assuming that the intra-component interactions g 11 and g 22 are unequal but close in value. The regime of large difference (g 11 − g 22 ) clearly might have other interesting dynamical features which are yet to be explored.
Secondly, in this work we have restricted ourselves to the mean field regime. While the mean field description captures well the richness of pattern formation phenomena (see [5,14,17,18,20] in addition to this work), it may be worth asking whether quantum effects beyond mean field might have interesting consequences for the patter dynamics generated by a π/2 pulse. For bosons in elongated traps, regimes other than mean field (such as the Lieb-Liniger or Tonks regimes) may occur naturally in experiments [29,34,[43][44][45]. Dynamics subsequent to a π/2 pulse in strongly interacting 1D gases outside the mean field regime is an open area for investigation.
Thirdly, we have assumed a spatially symmetric trap and an initial condition with spatial symmetry, and this plays a crucial role in the selection of instability channels. In a reallife experiment, the trap will have some left-right asymmetry. Also, thermal and quantum fluctuations can initiate spatially antisymmetric excitations. The extent to which a small spatial asymmetry affects spin dynamics remains unexplored; in such a case we would have some type of competition between the two types of instabilities. Navarro et al [14] have studied dynamical effects of fluctuations (noise), but the effects of thermal and quantum fluctuations are yet to be studied in the context of a π/2 protocol. Finally, one could consider time evolution and spatiotemporal patterns generated by a π/2 pulse in the presence of an optical lattice, described by the dynamics of a two-component Bose-Hubbard model. This is a situation easy for us to imagine realizing experimentally. One could speculate that there is a complex interplay between spin dynamics and the spatial arrangement of Mott and superfluid regions.