Spin modulation instabilities and phase separation dynamics in trapped two-component Bose condensates

Two-component Bose–Einstein condensates (BECs) are increasingly appreciated as a rich and versatile source of intricate non-equilibrium pattern dynamics phenomena. In addition to experimental observations [1–13], pattern dynamics in two-component BECs has also attracted significant theoretical interest (see, e.g., [14–28] and references cited in [14]).

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 ⁸⁷Rb and then using a π/2 pulse to switch half the atoms to a different hyperfine state [1–9]. This results in a binary condensate where the two intra-species interactions (g₁₁ and g₂₂) and one inter-species interaction (g₁₂) are all slightly different from each other, but the starting state is the ground state determined by g₁₁ 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₁₂) and the strength λ of the trapping potential. We provide a stability analysis mapping out regions of the λ–g₁₂ parameter space hosting different types of dynamics. Since it is now routine to monitor real-time 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 self-repulsions of the two species, i.e. if

$$\begin{equation} g_{11} g_{22} ~<~ \left(g_{12}\right)^2, \end{equation} \tag{ 1 }$$

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–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–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–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₁₂, 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₁₁ and g₂₂ are close but unequal: (a) the stability analysis is performed for g₁₁ = g₂₂ and their difference serves only to select appropriate instability modes; (b) the simulations are performed with g₂₂/g₁₁ = 1.01.

We restrict ourselves to repulsive interactions (g_ij > 0). Attractive interactions (g_ij < 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₁₂.

In section 4, we provide a relatively detailed account of the time evolution. For each trap strength λ, for values of g₁₂ not much larger than g₁₁, we observe simple collective modes. Above a threshold value of g₁₂, 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₁₂ 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.

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–37]:

$$\begin{eqnarray} \mathrm{i}\partial_t \psi_1 &=& \left( -\frac{1}{2} \partial^2_x +\frac{1}{2}\lambda^2x^2 +g_{11}|\psi_1|^2 +g_{12}|\psi_2|^2 \right) \psi_1, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \mathrm{i}\partial_t \psi_2 &=& \left( -\frac{1}{2}\partial^2_x +\frac{1}{2}\lambda^2x^2 +g_{12}|\psi_1|^2 +g_{22}|\psi_2|^2 \right) \psi_2. \end{eqnarray} \tag{ 3 }$$

Condensate wave functions ψ₁(x,t) and ψ₂(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_ij [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–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₁₁ = 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⁻⁵ in these units. For the purposes of this paper, 'tight' and 'shallow' mean, respectively, λ ≳ 10⁻² and λ ≲ 10⁻⁴ in our units. Note that the trap oscillator strength is λ^−1/2.

Of course, one can switch between different units via the transformation: x → x/l, t → t/l², λ → λl², g → gl and $\psi \rightarrow \psi \sqrt {l}$, 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₁₁, 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₁₁ = g₂₂ = g₁₂. The π/2 pulse may then be regarded as a sudden change (a quantum quench [39]) of the interaction parameters g₂₂ and g₁₂.

We use g₁₁ = g₂₂ for the stability analysis of section 3. For the explicit time evolution reported in section 4, we use g₁₁ and g₂₂ values close but unequal: g₁₁ = 1, g₂₂ = 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₁₁ and g₂₂ is somewhat larger (in the common case using ⁸⁷Rb hyperfine states $\left |{1}\right .\rangle =\left |{F=1,m_F=-1}\right .\rangle$ and $\left |{2}\right .\rangle =\left |{F=2,m_F=1}\right .\rangle$); however, our insights should be relevant to a broad regime of possible experiments. A full exploration of the regime of arbitrary differences (g₁₁ − g₂₂) 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].

We provide in this section a stability analysis for g₁₁ = g₂₂ that maps out the regions of λ–g₁₂ 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 ψ₀(x) the lowest-energy spatially symmetric stationary state of the case g₁₁ = g₂₂, with parameter g₁₂ set to its final value. (For large g₁₂, 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₁₁ + g₁₂. 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₁₁ ≠ g₂₂, because stationary states for such a case typically do not overlap completely in space. Instead, in our approach the difference between g₁₁ and g₂₂ 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₁₁ and g₂₂.

We linearize equations (2) and (3) around the reference stationary state ψ₀(x):

$$\begin{eqnarray} &&\psi_1(x, t) = \left[ \psi_{0}(x) + \delta\psi_1(x,t) \right] \exp(-\mathrm{i} \mu t),\nonumber\\\\ &&\psi_2(x, t) = \left[ \psi_{0}(x) + \delta\psi_2(x,t) \right] \exp(-\mathrm{i} \mu t),\nonumber \end{eqnarray} \tag{ 4 }$$

where μ is the chemical potential corresponding to the reference state. By keeping only terms of the first order in δψ₁(x,t) and δψ₂(x,t), we obtain a system of linear equations which can be cast in the form

$$\begin{equation} \partial_t^2 \left(\begin{array}{c} \delta\psi_1+ \delta\psi_1^* \\ \delta\psi_2 + \delta\psi_2^* \end{array} \right) +\mathcal{M} \left( \begin{array}{c} \delta\psi_1 + \delta\psi_1^* \\ \delta\psi_2 + \delta\psi_2^* \end{array} \right) = 0. \end{equation} \tag{ 5 }$$

Here $\mathcal {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 ψ₀(x) and expanding in the basis of harmonic trap (non-interacting) eigenstates.

Since we use g₁₁ = g₂₂ for the stability analysis, eigenmodes will have well-defined 'species parity', i.e. will all be either even [δψ₁(x,t) = δψ₂(x,t)] or odd [δψ₁(x,t) = −δψ₂(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₁₁ + g₁₂. Odd modes are more interesting—they describe out-of-phase 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 single-component 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 ω² = 3λ². 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.

Zoom In Zoom Out Reset image size

The two lowest-lying eigenmodes are odd modes encoding out-of-phase physics. For g₁₂ ≳ 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₁₂, thus leading to the onset of modulation instabilities. The antisymmetric mode becomes unstable at smaller value of g₁₂ (g^a₁₂ ≈ 1.6 for λ = 0.2) in comparison with the symmetric mode (g^s₁₂ ≈ 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₁₁ and g₂₂ naturally excites the second (spatially symmetric) mode. Thus, the second mode, occurring at larger g₁₂, 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⁻³ (also shown in figure 1) the lowest instability sets in for g^a₁₂ ≈ 1.02, while the next one appears at g^s₁₂ ≈ 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₁₂ and g^s₁₂) as a function of trap strength λ. For very shallow traps, the two transition lines merge as g^s₁₂ ≈ g^a₁₂ ≈ 1. The lower transition line (g^a₁₂) 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₁₂ value, given by the g^s₁₂ line.

Zoom In Zoom Out Reset image size

In section 4, we will see that spin pattern dynamics is indeed only generated when the inter-component repulsion g₁₂ exceeds the second instability line (g₁₂ > g^s₁₂), and that crossing the first instability (g^a₁₂ < g₁₂ < g^s₁₂) is not enough for pattern formation in a spatially symmetric trap.

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₁₁ = g₂₂ = g₁₂ = 1. The subsequent dynamics is performed with g₁₁ = 1, g₂₂ = 1.01, and several different values of g₁₂ 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

$$\begin{equation} w_{1,2}^2(t)=\int_{-\infty}^{\infty} x^2 |\psi_{1,2}(x,t)|^2\,\mathrm{d}x, \end{equation} \tag{ 6 }$$

while figure 4 shows density plots of the density difference (spin density), |ψ₁(x,t)|² − |ψ₂(x,t)|². 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For each λ the four values of g₁₂ from table 1 are used for figures 3 and 4. We have chosen g₁₂ values such that the first panel in each row is in the parameter region where there are no instabilities, the second one is in the region where the only instability is the antisymmetric one, and the third on each row is at g₁₂ values just above the second, relevant, instability. The fourth panel on each row is at higher g₁₂ values. The choice of g₁₂ 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₁₂ = 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.

Table 1. Parameters from the first five columns are used for the plots in figures 3 and 4. The instability values g^a₁₂ and g^s₁₂ (introduced in figures 1 and 2 and discussed in section 3) are also given for each trap strength.

λ	$\frac {g_{12}^{(1)}}{\sqrt {g_{11} g_{22}}}$	$\frac {g_{12}^{(2)}}{\sqrt {g_{11} g_{22}}}$	$\frac {g_{12}^{(3)}}{\sqrt {g_{11} g_{22}}}$	$\frac {g_{12}^{(4)}}{\sqrt {g_{11} g_{22}}}$	ga12	gs12
10−1	1.3	1.8	2	2.3	1.37	1.92
10−2	1.08	1.17	1.25	1.5	1.085	1.23
10−3	1.01	1.04	1.06	1.3	1.018	1.050
10−4	1.003	1.01	1.02	1.12	1.004	1.011
10−5	1	1.005	1.03	1.08	≈1	≈1

Broadly speaking, we note that there is only regular (collective-mode) dynamics in the second-column figures (g^a₁₂ < g₁₂ < g^s₁₂) 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₁₂) is the relevant one. The fourth panel on each row is at higher g₁₂ values, showing more rich dynamics.

In figure 3, we show the time dependence of the individual widths (w₁, w₂) and also of the total root-mean-square width, $w(t) = \sqrt {(w_1^2(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\pi /\sqrt {3}\lambda \approx 3.63/\lambda$. This follows from the frequency of this mode being almost constant near $\sqrt {3}\lambda$.

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₁₂ such that this mode has small real frequencies, this is excited as a regular 'spin' mode. For example, at λ = 10⁻³ and g₁₂ = 1.04, we observe an out-of-phase 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₁₂ 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₁₂ 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' |ψ₁(x,t)|² − |ψ₂(x,t)|². 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.

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₁₂, 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.

Zoom In Zoom Out Reset image size

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₁₂/g₁₁ − 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₁₂, we get more intricate dynamical patterns.

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₁₁ = g₂₂. The π/2 pulse of the experiments (in the cases where g₁₁ ≠ g₂₂) can be considered as turning on a nonzero (g₁₁ − g₂₂), 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₁₁ = g₂₂ to analyze dynamics with g₁₁ ≠ g₂₂, an obvious question is how the ratio g₂₂/g₁₁ affects the regime of applicability of this scheme. We expect that features of this (g₁₁ = g₂₂) stability analysis are useful for dynamical predictions as long as g₁₂/g₁₁ − 1 is roughly more than g₂₂/g₁₁ − 1. For example, for shallow traps (small λ), the instabilities occur at g₁₂/g₁₁ − 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₁₁ = g₂₂ = g₁₂, while the reference state is the lowest-energy spatially symmetric stationary state corresponding to the final value of g₁₂. 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₁₂, as opposed to the changes of g₁₂ 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₁₂ parameter space, while assuming that the intra-component interactions g₁₁ and g₂₂ are unequal but close in value. The regime of large difference (g₁₁ − g₂₂) 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–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 real-life 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.

IV acknowledges discussions with A Balaž and support from the Ministry of Education and Science of the Republic of Serbia, under project no. ON171017. Numerical simulations were run on the AEGIS e-Infrastructure, supported in part by FP7 projects EGI-InSPIRE, PRACE-1IP, PRACE-2IP and HP-SEE. NJvD acknowledges support from FOM and NWO.