Spin modulation instabilities and phase separation dynamics in trapped two-component Bose condensates
aa r X i v : . [ c ond - m a t . qu a n t - g a s ] S e p Spin modulation instabilities and phase separation dynamics in trapped two-component Bosecondensates
Ivana Vidanovi´c,
N. J. van Druten, and Masudul Haque Scientific Computing Laboratory, Institute of Physics Belgrade,University of Belgrade, Pregrevica 118, 11080 Belgrade, Serbia Institut f¨ur Theoretische Physik, Goethe-Universit¨at, 60438 Frankfurt/Main, Germany Max Planck Institute for the Physics of Complex Systems, N¨othnitzer Str. 38, 01187 Dresden, Germany Van der Waals-Zeeman Institute, University of Amsterdam,Science Park 904, 1098 XH Amsterdam, The Netherlands
In the study of trapped two-component Bose gases, a widely used dynamical protocol is to start from theground state of a one-component condensate and then switch half the atoms into another hyperfine state. Theslightly 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 trapstrength and by the inter- to intra-component interaction ratio. A stability analysis suited to the trapped situationprovides us with a framework to explain the various types of dynamics in different regimes.
PACS numbers: 67.85.-d, 67.85.De, 67.85.Jk, 03.75.Kk, 03.75.Nt
I. INTRODUCTION
Two-component Bose-Einstein condensates (BECs) are in-creasingly appreciated as a rich and versatile source of intri-cate non-equilibrium pattern dynamics phenomena. In addi-tion to experimental observations [1–13], pattern dynamics intwo-component BECs has attracted significant theoretical in-terest (see, e.g., [14–28] and citations in [14]).In a number of two-component BEC experiments reportedover more than a decade, a standard technique has been to startfrom the equilibrium state of a single-component BEC, e.g.,populating a single hyperfine state of Rb, and then usinga π/ pulse to switch half the atoms to a different hyperfinestate [1–9]. This results in a binary condensate where the twointra-species interactions ( g and g ) and one inter-speciesinteraction ( g ) are all slightly different from each other, butthe starting state is the ground state determined by g alone.Since it has been realized several times in several differentlaboratory setups, this is a paradigm non-equilibrium initialstate for binary condensate dynamics. A thorough and gen-eral analysis of the dynamics subsequent to such a π/ pulseis thus clearly important. In this article we present such ananalysis, clarifying the combined role of the inter-species in-teraction ( g ) and the strength λ of the trapping potential. Weprovide a stability analysis mapping out regions of the λ – g parameter space hosting different types of dynamics. Since itis now routine to monitor real-time dynamics in such experi-ments (e.g. [6]), we also directly analyze the real-time evolu-tion after a π/ pulse.It is widely known that the ground state of a uniform two-species BEC is phase separated or miscible depending onwhether or not the inter-species repulsion dominates over theself-repulsions of the two species, i.e., if g g < ( g ) (1)then the ground state is phase separated [15]. This criterionis also a key ingredient in understanding dynamical featuressuch as pattern dynamics in the density difference between the two species — such “spin patterns” emerge when the phaseseparation condition is satisfied. This can be understood asthe onset of a modulation instability [16–18], identified by theappearance of an unstable mode in the excitation spectrumaround a reference stationary state. For a homogeneous situa-tion, linear stability analysis shows that modulation instabilitysets in when the condition of Eq. (1) is satisfied [16–18] .The situation is different in the presence of a trapping po-tential. Phase separation in the ground state, as well as the ap-pearance of modulation instability when starting from a mixedstate, now requires larger inter-species repulsion [14, 19].This suggests that the region of parameter space where pat-tern dynamics occurs also depends on the trap. A trap is al-most always present in cold-atom experiments, and it is easyto imagine experiments where the trapping potential is not ex-tremely shallow but varies between tight and shallow limits.It is thus necessary to examine the relevance of Eq. (1) fortrapped binary BECs. To this end, we explore different trapstrengths spanning several orders of magnitude, and identifythe appropriate extensions of Eq. (1) for the type of spin dy-namics resulting from the π/ protocol described above.We focus on the effects of two parameters. First, we studyeffects of changing cross-species interaction g , thus gener-alizing Eq. (1) for trapped situations. Second, we explore therole of the relative strength of the trap with respect to the inter-actions. Our analysis, performed for a one-dimensional (1D)geometry, sheds light on the situation where g and g areclose but unequal: (a) the stability analysis is performed for g = g and their difference serves only to select appropri-ate instability modes; (b) the simulations are performed with g /g = 1 . .In Section II, we introduce the formalism and geometry. InSection III, we show results from a linear stability analysisfor a sequence of trap strengths, and identify and analyze rel-evant modulation instabilities. Through an analysis of unsta-ble modes, we present a classification of the parameter spaceinto dynamically distinct regions, in relation to the prototyp-ical initial state explained above. This may be regarded as adynamical “phase diagram”. A remarkable aspect is that the“phase transition” line most relevant to spin pattern dynam-ics does not arise from the first modulation instability (studiedin Ref. [14]). This first instability mode is antisymmetric inspace, and as a result is not naturally excited in a symmet-ric trap with symmetric initial conditions. Complex dynamics(not due to collective modes but rather due to modulation in-stability) is generated only when the first spatially symmetric mode becomes unstable, which occurs at a higher value of g .In Section IV we provide a relatively detailed account ofthe time evolution. For each trap strength λ , for values of g not much larger than g , we observe simple collectivemodes. Above a threshold value of g , the oscillation ampli-tude becomes sharply stronger, and at the same time the mo-tion becomes notably aperiodic, signaling that the dynamicsis more complex than a combination of a few modes. Dy-namical spin patterns start appearing at this stage and becomemore pronounced as g is increased further. The thresholdvalue at which the dynamics changes sharply corresponds tothe second modulation instability line rather than the first, aswe demonstrate through careful choice of parameters in eachregion of the phase diagram derived from stability analysis.Some further connections between the stability analysis anddynamical features, relating to the length scale of generatedpatterns, appear in Section V. In the concluding Section VIwe place our results in context and point out open questions. II. GEOMETRY AND FORMALISM
The relevant time-resolved experiments have been per-formed in both quasi-1D geometries (highly elongated trapswith strong radial trapping) [6] and in a 3D BEC of cylindri-cal symmetry with the radial variable playing analogous roleas the 1D coordinate [2, 5]. Since the basic phenomena arevery similar, we expect the same theoretical framework to de-scribe the essential features of each case. For definiteness, inthis work we show results for 1D geometry. We expect thegeneral picture emerging from this work to be qualitativelytrue also for other geometries exhibiting the same type of spindynamics.We describe the dynamics in the mean field framework atzero temperature, i.e., by two coupled Gross-Pitaevskii equa-tions [30–32]: i∂ t ψ = (cid:18) − ∂ x + λ x + g | ψ | + g | ψ | (cid:19) ψ , (2) i∂ t ψ = (cid:18) − ∂ x + λ x + g | ψ | + g | ψ | (cid:19) ψ . (3)Condensate wave functions ψ ( x, t ) and ψ ( x, t ) are normal-ized to unity, and λ is the strength of the harmonic trap. Fac-tors of particle number and radial trapping frequency are ab-sorbed as appropriate into the effective 1D interaction param-eters g ij [6, 32, 33]. We consider purely non-dissipative dy-namics, i.e., we do not attempt to model experimental lossrates with a phenomenological dissipative term as done in,e.g., Refs. [5–7].The equations above are in dimensionless form because wemeasure lengths in units of trap oscillator length and time in g −10123 ω g . . ω x δ ψ & δ ψ asymmetric-mode eigenvectorsymmetric-mode eigenvector g a g s λ = 10 −3 λ = 0.2 BREATHING DIPOLE
FIG. 1. (Color online.) Results from stability analysis. Squaredeigenvalues ω of the stability matrix M are plotted against g , for atight trap (top) and for a shallow trap (bottom left). The arrows showthe values of g for onset of the two instabilities, namely g a (onsetof spatially antisymmetric modulation instability) and g s (onset ofspatially symmetric instability). Typical eigenvectors correspondingto these two modes are shown in the panels on lower right. units of inverse trapping frequency, for a hypothetical trap ofunit strength ( λ = 1 ). The scale for trap strengths is itselffixed by imposing g = 1 . With this convention, small val-ues of λ correspond to a BEC in the Thomas-Fermi limit. Forcomparison, we note that the parameters of the experiment ofRef. [6] corresponds to λ of order − in these units. Ofcourse one can switch between different units via the transfor-mation: x → x/l , t → t/l , λ → λl , g → gl , and ψ → ψ √ l ,where l is an appropriately chosen scale.The initial state after a π/ pulse involves both componentsoccupying the ground state of a single-component system ofinteraction g , because the atoms were all in the first hyper-fine state before the pulse. We model this initial situation as atwo-component BEC with g = g = g . The π/ pulsemay then be regarded as a sudden change (a quantum quench [29]) of the interaction parameters g and g .We use g = g for the stability analysis of SectionIII. For the explicit time evolution reported in Section IV,we use g and g values close but unequal: g = 1 , g = 1 . . This choice of close values is convenient for illus-trating the structure of the phase diagram, especially for shal-low traps. In rubidium experiments the difference between g and g is somewhat larger (in the common case using Rb hyperfine states | i = | F = 1 , m F = − i and | i = | F = 2 , m F = 1 i ); however our insights should be relevantto a broad regime of possible experiments. A full explorationof the regime of arbitrary differences ( g − g ) remains anopen task, beyond the scope of the present manuscript.Numerical simulations presented in Section IV were per-formed using a semi-implicit Crank-Nicolson method [34,35]. III. STABILITY ANALYSIS AND DYNAMICAL “PHASEDIAGRAM”
We provide in this section a stability analysis for g = g that maps out the regions of λ - g parameter space whichsupport pattern formation instabilities.Ideally, one might like to perform a stability analysisaround the initial state. However, in contrast to the homo-geneous case [16], we are faced with the situation that theinitial state is not a stationary state of the final Hamiltonian.The choice of reference state is therefore a somewhat subtleaspect of the present analysis.We use as reference state ψ ( x ) the lowest-energy spa-tially symmetric stationary state of the case g = g , withparameter g set to its final value. (For large g , this isnot the ground state for these parameters, which is phase-separated.) This reference state has the advantage of lookingrelatively similar to our actual initial state (two componentstotally overlapping in space), and of being a stationary stateof the Hamiltonian for which we analyze linear stability. Ourreference state can be regarded as placing both components inthe single-component ground state for interaction g + g .We are not aware of a suitable stationary state even more sim-ilar to the actual initial state. We will see that our stabilityanalysis around this reference state will predict remarkablywell the main observed time-evolution features described inSection IV.Note that it is not natural to use g = g , because station-ary states for such a case typically do not overlap completelyin space. Instead, in our approach the difference between g and g will play the important role of selecting certain in-stability modes over others. For this reason, inferences fromthe present analysis apply only to small relative differencesbetween g and g .We linearize Eqs. (2) and (3) around the reference station-ary state ψ ( x ) : ψ ( x, t ) = [ ψ ( x ) + δψ ( x, t )] exp( − iµt ) ,ψ ( x, t ) = [ ψ ( x ) + δψ ( x, t )] exp( − iµt ) , (4)where µ is the chemical potential corresponding to the ref-erence state. By keeping only terms of the first order in δψ ( x, t ) and δψ ( x, t ) , we obtain a system of linear equa-tions which can be cast in the form: ∂ t (cid:18) δψ + δψ ∗ δψ + δψ ∗ (cid:19) + M (cid:18) δψ + δψ ∗ δψ + δψ ∗ (cid:19) = 0 . (5)Here M is a matrix differential operator which, upon dis-cretization or upon expansion in a set of orthogonal functions,becomes the so-called stability matrix (e.g., [22, 24]). We an-alyze below the eigenmodes of the stability matrix, which wehave obtained by numerically calculating the reference sta-tionary state ψ ( x ) and expanding in the basis of harmonictrap (non-interacting) eigenstates.Since we use g = g for the stability analysis, eigen-modes will have well-defined “species parity”, i.e. will allbe either even [ δψ ( x, t ) = δψ ( x, t ) ] or odd [ δψ ( x, t ) = − δψ ( x, t ) ] with respect to the interchange of species. Even λ g asymmetric instability, g symmetric instability, g x x | ψ | & | ψ | x x | ψ | & | ψ | < −−−−− −−−− −−−− −−−− −−−− −−−−− − − − − − − − −−− −−− < −−−−−− −−−−− −−−−− −−−−− −−−− −−−− −−−− −−−− −−− < −−−−− −−− −−−− − − − − − − − − −−− −−− −−− FIG. 2. (Color online.) The dynamical “phase diagram” showingthe critical values of g for the onset of the two types of modulationinstability versus the trap strength λ . The instability lines are shownas straight lines joining numerically determined values. The oscilla-tion schematics on the right (and corresponding arrows) indicate thatleft-right-left oscillation modes are persistent everywhere above the g a line, while in-out-in modes are persistent only above the higher g s line. The spatially symmetric instability ( g s line) is the one rele-vant for experimental situations with symmetric traps. Squares markvalues used in the dynamical simulations of Figs. 3 and 4 (Table 1). modes describe in-phase motion of the two components andsimply correspond to the excitation spectrum of a single-component BEC with interaction constant g + g . Oddmodes are more interesting — they describe out-of-phase mo-tion of two components and are therefore reflected in the spindynamics. Additionally, due to the spatial inversion symme-try x → − x , the solutions will also have well-defined spatialparity, and we can distinguish spatially symmetric and anti-symmetric modes.Typical eigenspectra are presented in Fig. 1. In the caseof a tight trap λ = 0 . , we notice two modes whose fre-quencies are nearly constant. These are even modes encod-ing single-component or in-phase physics. The lower one isthe dipole (Kohn) mode with frequency equal to the trap fre-quency λ . The second nearly-constant mode is the breath-ing mode, which for elongated traps takes value close to ω = 3 λ . The breathing mode (oscillations of cloud size)is visible in the plots of Fig. 3 (Section IV) as a fast oscilla-tion of the total condensate widths.The two lowest-lying eigenmodes are odd modes encod-ing out-of-phase physics. For g & , their frequencies aresignificantly below the breathing mode, and therefore lead torelatively slow oscillations in the spin density. This will alsobe visible in the real-time dynamics presented in Section IV(first two columns of Figs. 3 and 4). The forms of the corre-sponding eigenvectors are shown in the lower right of Fig. 1.The nature of the eigenvectors shows that the motion relatedto the lowest mode corresponds to the left-right oscillations ofthe two species, while the next odd mode corresponds to spa-tially symmetric spin motion. The frequencies of these twomodes become imaginary at certain values of g , thus lead-ing to the onset of modulation instabilities. The antisymmetricmode becomes unstable at smaller value of g ( g a ≈ . for λ = 0 . ) in comparison to the symmetric mode ( g s ≈ . for λ = 0 . ). In a spatially symmetric trap, there is no naturalmechanism for exciting the spatially antisymmetric mode. Onthe other hand, any difference between g and g naturallyexcites the second (spatially symmetric) mode. Thus, the sec-ond mode, occurring at larger g , is the relevant instabilityfor understanding the dynamics observed in experiments andexplored numerically in Section IV.We find similar excitation spectra for trap strengths λ span-ning several orders of magnitude. The spatially antisymmetricmode 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 Fig. 1) the lowestinstability sets in for g a ≈ . , while the next one appearsat g s ≈ . . The distinction between two instabilities be-comes ever smaller as we go toward a uniform system λ → ,where the phase-separation condition Eq. (1) becomes exact.Nevertheless, even for shallow traps, the issue is not purelyacademic as the precision in experimental measurement andcontrol of scattering lengths continues to improve [6, 36].In Fig. 2 (main panel), the results of the stability analysesare combined to present a dynamical “phase diagram”. Thetwo lines show the two instabilities ( g a and g s ) as a functionof trap strength λ . For very shallow traps, the two transitionlines merge as g s ≈ g a ≈ . The lower transition line ( g a )was previously introduced in Ref. [14]. However, for a trapand initial state with left-right spatial symmetry, this is not therelevant dynamical transition line, because the first even modeonly becomes unstable at some higher g value, given by the g s line.In the next Section, we will see that spin pattern dynamicsis indeed only generated when the inter-component repulsion g exceeds the second instability line ( g > g s ), and thatcrossing the first instability ( g a < g < g s ) is not enoughfor pattern formation in a spatially symmetric trap. IV. DYNAMICAL FEATURES ACROSS THE PARAMETERSPACE
In this Section we present and analyze the dynamicsobtained from direct numerical simulation of the Gross-Pitaevskii equations (2) and (3), after the system is initiallyprepared in the ground state of the situation g = g = g = 1 . The subsequent dynamics is performed with g =1 , g = 1 . , and several different values of g for each trapstrength λ .It is difficult to show the full richness of pattern dynam-ics through plots of a few quantities. We choose to show thedynamics through two types of plots (Figs. 3 and 4). Fig. 3shows the time dependence of the root mean square widths ofthe two components w , ( t ) = Z ∞−∞ x | ψ , ( x, t ) | dx, (6)while Fig. 4 shows density plots of the density difference (spindensity), | ψ ( x, t ) | − | ψ ( x, t ) | . In both figures, each row λ g (1)12 √ g g g (2)12 √ g g g (3)12 √ g g g (4)12 √ g g g a g s − − − − − ≈ ≈ g a and g s (intro-duced in Figs. 1 and 2 and discussed in Section III) are also given foreach trap strength. corresponds to a different trap strength ( λ ), and we approachthe shallow trap (Thomas-Fermi) limit going from top to bot-tom.For each λ the four values of g from Table I are used forFigs. 3 and 4. We have chosen g values such that the firstpanel in each row is in the parameter region where there areno instabilities, the second one is in the region where the onlyinstability is the antisymmetric one, and the third on each rowis at g values just above the second, relevant, instability. Thefourth panel on each row is at higher g values. The choiceof g values with respect to instability lines is clear in thetighter traps of the top three rows, as also shown by squares inFig. 2. For shallow traps (lower rows), the instability lines aretoo close together and too close to g = 1 , so making suchchoices is not meaningful. In the following, as we comparefeatures of the different columns, we implicitly exclude thelowest row (smallest λ ). This is also indicated by the fact thatthe schematic instability lines in Figs. 3 and 4 are not extendedto the lowest row.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. Thereis generally a sharp difference between the second and thirdfigure 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 Fig. 3, we show time-dependence of the individualwidths ( w , w ) and also of the total root mean square width, w ( t ) = p ( w ( t ) + w ( t )) / . Consistent with our observa-tion that spatially symmetric modes (and not the antisymmet-ric ones) are naturally excited in the current setup, the dy-namics shows signatures of the two most prominent spatiallysymmetric modes noted in Fig. 1. The breathing mode is theeasiest to notice and most ubiquitous — it shows up in almostevery parameter choice as oscillations in the total density (in-phase in the two components), with a typical period given by π/ √ λ ≈ . /λ . This follows from the frequency of thismode being almost constant near √ λ .We also see out-of-phase motion of the two components,associated with the lower spatially-symmetric mode in Fig.1, which has odd species parity. In the first two columns ofFig. 3, corresponding to smaller values of g such that thismode has small real frequencies, this is excited as a regular λ t λ t λ t λ t g increasing λ= − λ= − λ= − λ= − λ= − T i gh t e r t o s h a ll o w e r t r a p s λ / w t ( ) λ / w t ( ) λ / w t ( ) λ / w t ( ) λ / w t ( ) FIG. 3. (Color online.) Time evolution of root-mean-square widths after π/ pulse (interaction quench). First component width w ( t ) is shown as blue dashed line, second component width w ( t ) is shown as red solid line (gray solid without color), the total width w ( t ) = p ( w ( t ) + w ( t )) / is the black solid line intermediate between the other two. From top to bottom: tight to shallow traps. For each trapstrength, four values of g (indicated near top of each panel) from Table I are used. The two lines separating first and second column (reddashed) and second and third column (black solid) indicate the ‘positions’ of instability lines, from Figure 2. While the first two columns lookqualitatively the same and show regular oscillatory dynamics, in the third column we observe aperiodic motion of stronger amplitude that werelate to the onset of spin pattern dynamics. The spin dynamics is even more pronounced in the fourth column. ‘spin’ mode. For example, at λ = 10 − and g = 1 . ,we observe an out-of-phase oscillation with the period of ap-proximately ≈ , much slower than the breathing mode. Inaddition, the oscillation period of the out-of-phase motion isslower in the second than in the first column of each row, cor-responding to the decreasing frequency of the mode, as seenin stability analysis (Fig. 1). Once g becomes large enoughthat the instability threshold for this mode is crossed, the os-cillation amplitudes increase sharply and the width dynam-ics becomes strongly aperiodic and irregular (third column ofFig. 3). This signifies the onset of pattern dynamics, as op-posed to the excitation of a regular collective mode around astable state. Irregularity of the width dynamics at stronger g is even more apparent in the fourth column of Fig. 3.It is noteworthy that the spatially antisymmetric modes playno role and do not show up in these dynamical simulations.We see no signature of the Kohn mode. Nor do we see anysharp change associated with the instability of the antisym-metric mode, i.e., there is no sharp difference between the firsttwo columns of Fig. 3.In Fig. 4 we show the dynamics of the “spin density” | ψ ( x, t ) | − | ψ ( x, t ) | . The case of very shallow traps (lastrow), resembles the data in Refs. [17, 18]. As in Fig. 3, thefirst two columns show regular oscillations, correspondingto collective modes without instability. A sharp change oc-curs, not across the first instability line (between 1st and 2ndcolumn), but instead across the second instability line (2ndand 3rd columns), especially for tighter traps (top three rows)where comparison with instability lines is meaningful. Thesharp change can be noted through the color scales, which isdramatically different between second and third columns inthe upper rows. V. LENGTH SCALES OF PATTERNS
In homogeneous stability analysis, the length scale of pat-terns is inferred from the wavevector (momentum) at which aninstability first occurs. Since we perform our stability analysisspecifically for trapped systems, we do not have a momentumquantum number. Nevertheless, the eigenvectors of the un-stable modes contain information about the form of patterns λ t λ t λ t λ t λ / x λ / x λ / x λ / x λ / x FIG. 4. (Color online.) Spin dynamics subsequent to the π/ protocol, shown via the density difference λ − / (cid:0) | ψ ( x, t ) | − | ψ ( x, t ) | (cid:1) .Traps and g values are the same as in Fig. 3 and Table I: λ decreases from − to − from top to bottom and g values are indicatednear top of each panel. As in Fig. 3, the black solid line and the red dashed line indicate the instability lines from the “phase diagram” ofFig. 2. Note the sharp change of color-scale ranges between second and third columns in the upper rows, indicating that the dynamics changesdramatically only across the second instability line. generated in the dynamics of the trapped system. This is il-lustrated in Fig. 5, where the eigenvectors of the lowest unsta-ble even-parity modes are shown for several values of g , to-gether with the spin patterns generated in the non-equilibriumdynamics. There is a close match between the distance be-tween nodes of the eigenvectors (rough analog of ‘wavevec-tor’) and the length scales involved in the patterns.In Fig. 4 we see that the patterns contain more spatial struc-ture in shallow traps. The top two rows (tight traps) onlyshow in-out-in type of patterns. This can be understood fromthe idea that the interactions induce length scales (‘healinglengths’) in the problem, which are smaller for larger interac-tions, and which set the length scale of spatial structures. Fortight traps, the healing length set by the interactions is largeor comparable to the cloud size, so that only global dynamicalpatterns are generated. In such traps, generation of complexpatterns with many spatial oscillations would require muchhigher values of ( g /g − . For shallow traps, the healinglength becomes much smaller than the cloud size; as a resultone can have a multitude of dynamical spin structures in thesystem, of the type seen in experiments and prior simulations[5, 6, 18]. This heuristic explanation can be made quantitativeby counting the number of nodes appearing in the eigenmodes(as in Fig. 5). VI. CONCLUSIONS AND OPEN PROBLEMS
In this article we have analyzed a widely used dynamicalprotocol for two-component BECs, which involves startingfrom the ground state of one component and switching halfthe atoms to a different component through a π/ pulse. Wehave presented a stability analysis suitable to the trapped sit-uation, and also presented results from extensive dynamicalsimulations. Through an analysis of unstable modes, we havepresented a classification of the parameter space into a numberof dynamically distinct regions, in relation to the prototypicalinitial state. This may be regarded as a dynamical “phase dia-gram”.In the ‘stable’ regime of parameter space (no modulationinstabilities), our stability analysis explains the observed slowspin oscillations compared to the fast breathing mode oscilla-tions of the total density. We demonstrate that the important“phase transition” line for spatially symmetric situations rele-vant to most experiments is not the first instability (studied inRef. [14]), but a second transition line. The first instability isantisymmetric in space, and as a result is not naturally excitedin a symmetric trap.Our stability analysis is performed relative to a stationarystate of the situation g = g . The π/ pulse of the ex-periments (in the cases where g = g ) can be considered as g =1.02 g =1.06 g =1.12 -0.2 0 0.2-8 0 8 λ x-0.2 0 0.2-8 0 8 λ x-0.15 0 0.15-8 0 8 λ x δ ψ & δ ψ FIG. 5. (Color online.) Top: Eigenvectors of the most unstableeven eigenmodes, from the stability analysis of Section III, for λ =10 − , and g = 1 . , . , and . from left to right. Below eacheigenvector, the corresponding spin dynamics after the π/ protocol(parameters of Section IV) is shown through the time evolution of | ψ ( x, t ) | − | ψ ( x, t ) | . turning on a nonzero ( g − g ) , i.e., turning on ‘buoyancy’such that one component gains more energy by being in theinterior of the trap compared to the other. This helps to selectinstability modes which are symmetric in space.Since we have used a stability analysis with g = g toanalyze dynamics with g = g , an obvious question is howthe ratio g /g affects the regime of applicability of thisscheme. We expect that features of this ( g = g ) stabil-ity analysis are useful for dynamical predictions as long as g /g − is roughly more than g /g − . For example, forshallow traps (small λ ), the instabilities occur at g /g − values comparable to 0.01, which is why the placement of pa-rameters in the three dynamical regions of the ‘phase diagram’(Fig. 2) is not meaningful for the smallest λ values (lowestrows of Figs. 3 and 4).For the stability analysis we used a reference stationarystate which is of course not the initial state: the initial state isthe ground state for g = g = g , while the reference stateis the lowest-energy spatially symmetric stationary state cor-responding to the final value of g . The instability lines foundin this stability analysis would describe even better a situationwhere the dynamics is triggered by a small quench of g , asopposed to the changes of g that we consider here, whichcan be relatively large. We have looked at some examples ofthis type of dynamics and indeed find instabilities matchingthe stability analysis extremely well. However, although the initial state in the π/ dynamics is somewhat different fromthe reference state of our stability analysis, our results showthat this stability analysis does provide an excellent overallpicture of the dynamics generated by the π/ protocol.Our work opens up a number of questions deserving furtherstudy. First of all, we have thoroughly explored the λ – g pa-rameter space, while assuming that the intra-component inter-actions g and g are unequal but close in value. The regimeof large difference ( g − g ) clearly might have other inter-esting dynamical features which are yet to be explored.Second, in this work we have restricted ourselves to themean field regime. While the mean field description cap-tures well the richness of pattern formation phenomena (c.f.Refs. [5, 14, 17, 18, 20] in addition to present work), it maybe worth asking whether quantum effects beyond mean fieldmight have interesting consequences for the patter dynam-ics generated by a π/ pulse. For bosons in elongated traps,regimes other than mean field (such as Lieb-Liniger or Tonksregimes) may occur naturally in experiments [37–41]. Dy-namics subsequent to a π/ pulse in strongly interacting 1Dgases outside the mean field regime is an open area of investi-gation.Third, we have assumed a spatially symmetric trap and aninitial condition with spatial symmetry, and this plays a cru-cial 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 spatiallyantisymmetric excitations. The extent to which a small spatialasymmetry affects spin dynamics remains unexplored; in sucha case we would have some type of competition between twotypes of instabilities. Ref. [14] has studied dynamical effectsof fluctuations (noise), but the effects of thermal and quan-tum fluctuations is yet to be studied in the context of a π/ protocol.Finally, one could consider time evolution and spatiotem-poral patterns generated by a π/ pulse in the presence of anoptical lattice, described by the dynamics of a two-componentBose-Hubbard model. This is a situation easy to imagine real-izing experimentally. One could speculate complex interplaybetween spin dynamics and the spatial arrangement of Mottand superfluid regions. ACKNOWLEDGMENTS
I.V. acknowledges discussions with A. Balaˇz and supportby the Ministry of Education and Science of the Republic ofSerbia, under project No. ON171017. Numerical simulationswere run on the AEGrandomIS e-Infrastructure, supported inpart by FP7 projects EGI-InSPIRE, PRACE-1IP, PRACE-2IP,and HP-SEE. NJvD acknowledges support from FOM andNWO. [1] M. R. Matthews, D. S. Hall, D. S. Jin, J. R. Ensher, C. E. Wie-man, E. A. Cornell, F. Dalfovo, C. Minniti, and S. Stringari, Phys. Rev. Lett. , 243 (1998).[2] D. S. Hall, M. R. Matthews, J. R. Ensher, C. E. Wieman, and E. A. Cornell, Phys. Rev. Lett. , 1539 (1998).[3] H.-J. Miesner, D. M. Stamper-Kurn, J. Stenger, S. Inouye,A. P. Chikkatur, and W. Ketterle, Phys. Rev. Lett. , 2228(1999).[4] H. J. Lewandowski, D. M. Harber, D. L. Whitaker, andE. A. Cornell, Phys. Rev. Lett. , 070403 (2002)[5] K. M. Mertes, J. W. Merrill, R. Carretero-Gonz´alez,D. J. Frantzeskakis, P. G. Kevrekidis, and D. S. Hall, Phys. Rev.Lett. , 190402 (2007).[6] P. Wicke, S. Whitlock, and N. J. van-Druten, arXiv:1010.4545.[7] R. P. Anderson, C. Ticknor, A. I. Sidorov, and B. V. Hall, Phys.Rev. A , 023603 (2009)[8] P. Boehi, M. F. Riedel, J. Hoffrogge, J. Reichel, T. W. Haensch,and P. Treutlein, Nat. Phys. , 592 (2009).[9] M. Egorov, R. P. Anderson, V. Ivannikov, B. Opanchuk,P. Drummond, B. V. Hall, and A. I. Sidorov, Phys. Rev. A ,021605 (2011).[10] S. B. Papp, J. M. Pino, and C. E. Wieman, Phys. Rev. Lett. ,040402 (2008).[11] C. J. Myatt, E. A. Burt, R. W. Ghrist, E. A. Cornell, andC. E. Wieman, Phys. Rev. Lett. , 586 (1997).[12] G. Modugno, M. Modugno, F. Riboli, G. Roati, and M. Ingus-cio, Phys. Rev. Lett. , 190404 (2002).[13] C. Hamner, J. J. Chang, P. Engels, and M. A. Hoefer, Phys. Rev.Lett. , 065302 (2011).[14] R. Navarro, R. Carretero-Gonz´alez, and P. G. Kevrekidis, Phys.Rev. A , 023613 (2009).[15] Tin-Lun Ho and V. B. Shenoy, Phys. Rev. Lett. , 3276 (1996).[16] E. Timmermans, Phys. Rev. Lett. , 5718 (1998).[17] K. Kasamatsu and M. Tsubota, Phys. Rev. Lett. , 100402(2004)[18] K. Kasamatsu and M. Tsubota, Phys. Rev. A , 013617 (2006).[19] L. Wen, W. M. Liu, Y. Cai, J. M. Zhang, and J. Hu, Phys. Rev.A , 043602 (2012).[20] S. Ronen, J. L. Bohn, L. E. Halmo, and M. Edwards, Phys. Rev.A , 053613 (2008). [21] J. Sabbatini, W. H. Zurek, and M. J. Davis, Phys. Rev. Lett. ,230402 (2011).[22] C. K. Law, H. Pu, N. P. Bigelow, and J. H. Eberly, Phys. Rev.Lett. , 3105 (1997).[23] Th. Busch, J. I. Cirac, V. M. Perez-Garcia, and P. Zoller, Phys.Rev. A , 2978 (1997).[24] H. Pu and N. P. Bigelow, Phys. Rev. Lett. , 1134 (1998).[25] R. Graham and D. Walls, Phys. Rev. A , 484 (1998).[26] A. Sinatra, P. O. Fedichev, Y. Castin, J. Dalibard, andG. V. Shlyapnikov, Phys. Rev. Lett. , 251 (1999).[27] Y. Li, P. Treutlein, J. Reichel, and A. Sinatra, Eur. Phys. Journ.B , 365 (2009).[28] A. Balaˇz and A. Nicolin, Phys. Rev. A , 023613 (2012).[29] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore,Rev. Mod. Phys. , 863 (2011).[30] L. P. Pitaevskii, Sov. Phys. JETP , 451 (1961).[31] E. P. Gross, Nuovo Cimento , 454 (1961).[32] L. Salasnich, A. Parola, and L. Reatto, Phys. Rev. A , 043614(2002).[33] M. Olshanii, Phys. Rev. Lett. , 938 (1998).[34] P. Muruganandam and S. K. Adhikari, Comput. Phys. Commun. , 1888 (2010).[35] D. Vudragovi´c, I. Vidanovi´c, A. Balaˇz, P. Muruganandam, andS. K. Adhikari, Comput. Phys. Commun. , 2021 (2012).[36] M. Egorov, B. Opanchuk, P. Drummond, B. V. Hall, P. Han-naford, and A. I. Sidorov, arXiv:1204.1591.[37] D. S. Petrov, G. V. Shlyapnikov, and J. T. M. Walraven, Phys.Rev. Lett. , 3745 (2000).[38] J. N. Fuchs, D. M. Gangardt, T. Keilmann, and G. V. Shlyap-nikov, Phys. Rev. Lett. , 150402 (2005).[39] A. H. van Amerongen, J. J. van Es, P. Wicke,K. V. Kheruntsyan, and N. J. van Druten, Phys. Rev. Lett. ,090402 (2008).[40] T. Kinoshita, T. Wenger, and D. S. Weiss, Science , 1125(2004).[41] E. Haller M. Gustavsson, M. J.. Mark, J. G. Danzl, R. Hart,G. Pupillo, and H.-C. Ngerl, Science325