Dynamics of spin and density fluctuations in strongly interacting few-body systems

The decoupling of spin and density dynamics is a remarkable feature of quantum one-dimensional many-body systems. In a few-body regime, however, little is known about this phenomenon. To address this problem, we study the time evolution of a small system of strongly interacting fermions after a sudden change in the trapping geometry. We show that, even at the few-body level, the excitation spectrum of this system presents separate signatures of spin and density dynamics. Moreover, we describe the effect of considering additional internal states with SU(N) symmetry, which ultimately leads to the vanishing of spin excitations in a completely balanced system.

symmetry after a sudden change -a quench -in the trapping potential. Recently, systems of cold atoms in optical lattices with SU(3) and SU(4) symmetries have been explored theoretically 43,44 . The correlations of SU(N) impenetrable systems have also been calculated for an increasing number of internal components and display interesting properties 45 .
We begin by presenting the formalism used to describe a system of strongly interacting atoms in a trap, where the Hamiltonian can be mapped into a spin chain with position-dependent exchange coefficients. We then describe the quench protocol, which essentially consists of changing the trap from a split well (where we assume a Gaussian barrier in the center of the system) to a simple harmonic well (see Fig. 1 for a sketch of this protocol). The ground-state configurations for these two systems are considerably different, and by changing the potential we can expect a non-trivial time evolution in the spin and density sectors. Initially, we describe the effect of the quench in the density sector and its consequences on the spin chain dynamics. By combining the dynamics of both sectors, we can extract the signatures of the separation between the density and spin oscillations in the system, showing how this effect can be observed in few-body ensembles, even as the number of internal components is increased. Moreover, we show that for a completely balanced system (where each atom is in a different internal state) the spin signature vanishes, and the excitation spectrum is analogous to that of a gas of impenetrable bosons. Multi-component cold atomic ensembles with strong interactions are currently within experimental reach 46 and often exhibit exotic dynamical effects, such as edge states 47,48 . In these systems, the internal states of the atoms can be manipulated with laser pulses, and the behavior of each component can be measured with precision. Studying cold atoms with different internal symmetries in a highly controllable environment can lead to insight on particle physics models and even shed light on the validity of unification schemes 49 . Parts of the work presented here have been published, with modifications, in ref. 50 .

Results
System description. Our goal is to describe the dynamics of a strongly interacting few-body system with internal ("pseudospin") degrees of freedom. We focus our description on a fermionic system, but the formalism is equally valid for bosons with the correct adaptations to the Hamiltonian. We consider initially an SU(2) system, where the internal degrees of freedom are described by |↑〉 and |↓〉. Later we will generalize our approach to systems with higher symmetries. We denote the particle numbers in each species by N ↑ and N ↓ , the total number of particles thus being given by N = N ↑ + N ↓ . For simplicity, we adopt the notation N ↑ + N ↓ (e.g. 3 + 1 for a system with three particles of species |↑〉 and one of species |↓〉). Experimentally, two-component fermionic systems can be realized by preparing trapped ensembles of 6 Li atoms in the two lowest hyperfine states 51,52 .
The Hamiltonian for the system under consideration is given by Sketch of the quench protocol adopted in this work. An imbalanced (N ↑ = 3, N ↓ = 1) strongly interacting two-component system of atoms is initialized in the ground state of a harmonic potential modified by a Gaussian bump centered at x = 0. For t > 0 this perturbation is suddenly switched off and the time evolution of the density for each component is registered. The majority (minority) density ρ ↑(↓) (x) is shown in green (blue). The black curves, which are rescaled for clarity, show the shape of the trapping potential at t = 0 and t > 0. (b,c) Time evolution of the spin densities, where length and time are given in units of the harmonic potential (see text for details). Yellow and blue colors indicate higher and lower densities, respectively. The slices shown in (a) correspond to t/T 0 = 0 and t/T 0 = 5 in the figures on the right. This figure has been published, with modifications, in ref. 50 .
www.nature.com/scientificreports www.nature.com/scientificreports/ is a trapping potential. The remaining term accounts for the contact interactions, where atoms in different internal states interact with strength g (interactions between atoms in the same internal state are forbidden due to the Pauli principle). Since we are dealing with atoms of the same element in different internal states, we consider all masses equal. In our calculations, we use the length, energy and time units of the harmonic trap, that is ω = l m / 0 0  , ε 0 = ω 0 and T 0 = 1/ω 0 , respectively. The interaction strength g is given in units of 2  /ml. We also assume, for simplicity, A case of particular interest in one-dimensional systems is the limit of strong interactions (g 1  ), where the Hamiltonian can be mapped to that of a spin chain [32][33][34][35]53,54 . In the case of fermions, Eq. (1) can be written, up to linear order in 1/g, as 36 where P i,i+1 is the permutation operator, which exchanges two neighboring atoms of different components. When all interactions between atoms in different internal states are the same, the system exhibits an SU(N) symmetry 55 . In the following sections, we focus only on strong repulsive interactions, and consider a fixed value of g in our calculations. For the formalism considered here, it has been shown that the static properties and spatial distributions are well-described already for g ~ 10 32 . In the strongly attractive limit, the model described here corresponds to the fermionic analog 56 of the Super-Tonks-Girardeau gas 57 .
The exchange coefficients α i are solely determined by the geometry of the trapping potential through the wave function for a system of spinless fermions, which we label Φ(x 1 , ..., x N ). Obtaining the eigenstates of Eq. (2), we can calculate the spatial distributions of each atomic component in the trap as , ↑ ↓ m i , is the probability of finding (↑, ↓) spins at site i and ρ i (x) represents the spatial distribution of each individual particle in the trap. In Section 3 we provide details on how to calculate the exchange coefficients and the spatial densities.

Dynamics.
We now describe the procedure that induces the dynamical evolution of the system, which consists of a sudden change of the trapping potential. Our initial choice of V(x) is given by a harmonic trap with an additional Gaussian bump in the center, as shown in Fig. 1(a). For t > 0, the bump is suddenly turned off. The time evolution of the spinless fermion wave function can then be described in terms of the evolution of the single particle orbitals under the same quench protocol 19,58 . We are thus able to construct Φ(x 1 , ..., x N , t) for all times t, which in turn determines the time-dependence of the exchange coefficients in Eq. (2). As a general rule, we can assume that the exchange coefficients are proportional to the overlap between the single particle distributions; additionally, since the trapping potentials are symmetric at all times, we have α i (t) = α N−i (t). In Fig. 2(a) we show the time evolution of the spatial densities obtained from Φ(x 1 , ..., x N , t). At this point, we are not considering the spin sector, so the densities shown correspond to the total ("charge") density Time evolution of (a) the total density ρ c (x, t) and (b) the exchange coefficients α 1 (t) (green) and α 2 (t) (yellow), for the case of N = 4. Since the trapping potential is spatially symmetric at all times, we have α 3 (t) = α 1 (t). In (b), the black dashed curves show the results obtained with an analytical fit of the data with a Fourier expansion. In (b), we assume g = 20. This figure has been published, with modifications, in ref. 50 . www.nature.com/scientificreports www.nature.com/scientificreports/ which is normalized to the total number of particles. Its dynamical behavior is what should be expected for the coherent density oscillations of a Tonks-Girardeau gas 58 . In panel (b), we see the behavior of α 1 (t)/g and α 2 (t)/g after the sudden change in the potential. We observe that the periodicity of the dynamics in the motion of the infinitely repulsive gas is reflected in the dynamics of the exchange coefficients. Particularly, since the system is approximately split in two parts at t = 0, we have that α(t = 0) ~ 0 (there is nearly no overlap between the initial densities at the center of the system). When the densities in the central region become larger, we see that the numerical value of α 2 (t) surpasses that of α 1 (t). This indicates that the spin correlations of the system should change between these two points, as we will see next.
Given the established time-periodicity of Hamiltonian (2), we can analyze the dynamics of the spin sector using Floquet theory (see Section 3 for details). While most studies performed in this context must deal with the issue of thermalization due to the external driving 59 , in our case the time-dependence of the spin chain originates directly from the dynamics of the density sector generated by the sudden change of the trapping potential. To describe the time evolution of the system in these terms, we first find an analytical fit of the exchange coefficients in terms of Fourier modes. In Fig. 2(b), this approximation is shown as the black dashed curves. The full time evolution of the spin sector is performed through the numerical integration of the Schrödinger equation with the Crank-Nicolson method 60 . Our quantities of interest are the dynamical densities ρ c (x, t), ρ ↑ (x, t) and ρ ↓ (x, t), as well as the squared width of the distribution for each spin density, defined as . While we choose to focus on the time evolution of parity-symmetric operators, other quantities (like, for instance, 〈x ↑,↓ (t)〉) could also be explored. Such cases, however, would require a different choice of quench protocol or initial spin state to break the symmetry across the system.

SU(2).
We choose initially a fermionic system with SU(2) symmetry. In this case, the permutation operator is , which allows us to write the Hamiltonian as where E 0 has the same meaning as in Eq. (2). We notice that this Hamiltonian reproduces the fermionic cases presented, for instance, in refs 33,37 . We focus mostly on the 3 + 1 case, which can be interpreted as a few-body Fermi background in the presence of an impurity 61 . As we will show next, in this simple setup it is already possible to describe the decoupling of spin and density dynamics. At the end of this section, we include some results for a 4 + 2 combination, which points to the possibility of realizing the protocol described here in larger systems.
The ground state of Hamiltonian (5) with repulsive interactions (g > 0) has antiferromagnetic correlations and can be described (in the 3 + 1 case) by aside from a normalization factor. Here, we have assumed a homogeneous potential, such that the exchange coefficients are identical and equal to 1. However, it can be shown that similar results hold for a harmonic trap. To observe how the ground state correlations change with the choice of α i , we define the operator P edge = |〈↓↑↑↑|gs〉| 2 + |〈↑↑↑↓|gs〉| 2 , which provides information regarding the position of the impurity (the ↓ atom) in the system. In Fig. 3 we show the values of P edge for different choices of α 1 and α 2 , where again we assume a parity-symmetric potential.
We can readily see that, for α 1 = α 2 , we have a constant result of P (2 2 ) edge 1 4 = − . Above the diagonal (α 2 > α 1 ) we have the region that includes, for instance, the coefficients obtained from a harmonic trap. In this case, the antiferromagnetic correlations are even more prevalent. If, however, , we obtain a larger probability of finding the impurity at the edges. These ground state correlations are obtained when considering a potential such as the double-well assumed for t = 0. In fact, we can plot the trajectory of the exchange coefficients after the quench in the trapping potential described above. This is shown by the dashed curve in Fig. 3. The white dot denotes the values at t = 0; the curve is then traveled back and forth periodically as t increases. The fact that this trajectory crosses over the diagonal indicates that the sudden quench in the potential should induce major changes in the spin correlations of the system.
To quantify this effect we now calculate the time evolution of the spin density and of the squared width for each component. It is important to notice that, for t > 0, the external potential is a simple harmonic trap. In Fig. 1(b,c), we show the time evolution of ρ ↑ (x, t) and ρ ↓ (x, t), respectively. We notice that, while the underlying dynamics seen in Fig. 2 is still present, we now have an additional oscillation mode. Specifically, after the sudden change in the potential, we observe a tendency of the majority atoms to spread to the edges, while the impurity localizes towards the center.
In Fig. 4(a) we show the time evolution of the squared width for the density of each component, over a larger time interval. This can be interpreted as induced breathing modes for the background and the impurity. Additionally, we show the dynamical behavior of the total density (Eq. (4)). Besides corroborating the results found in Fig. 1(b,c), these curves show how the fluctuations in the density and spin sectors are captured as two oscillations modes in the dynamics of the spin densities for each individual component. In Fig. 4(b,c) we show the Fourier transform of the width oscillations, defined as www.nature.com/scientificreports www.nature.com/scientificreports/ spin and density excitations appear as two separate peaks, the lower frequency corresponding to the spin dynamics.
Here, we can see that the dynamics of the minority component is strongly dominated by the spin excitations. On the other hand, the majority component has a more balanced distribution of oscillations in the density and spin sectors. In Fig. 4(b,c), we additionally include the theoretical predictions for the density and spin oscillations (as black and gray dashed lines, respectively). The first is obtained by calculating the oscillation frequency of the single particle spatial orbitals and corresponds to the expected result for the Tonks-Girardeau or the spinless fermion gas, namely ω = 2ω 0 . The second is extracted by calculating the gaps in the Floquet quasienergy spectrum of the time-periodic spin chain. We point out that the coherence in the oscillations of each component over long times is due to the assumption that the spatial sector is described by a spinless fermion wave function, which in a real system is only approximately true for very strong interactions. Depending on the real value of the parameter g, we could expect to find excitation peaks which are slightly different than the ones predicted here. In our formalism, modifying g means changing the energy gaps in Eq. (5), which translates into a frequency shift in the spin excitation peaks (with the frequency decreasing as the interaction strength is increased).
In Fig. 4(d-f), we extend our results to the 4 + 2 case. The time-dependence of the exchange coefficients for N = 6 is obtained by calculating these coefficients for one or two periods and then fitting a function by expanding  www.nature.com/scientificreports www.nature.com/scientificreports/ the result in Fourier modes, as done in Fig. 2 for the N = 4 case. While the general behavior (relative magnitude and position of the excitation peaks) are maintained with respect to the N = 4 system, we find additional peaks in the low-frequency side of the spectrum, which originates from having a larger number of particles in each component. Nevertheless, these excitations can also be captured by analyzing the Floquet quasienergy gaps of the time-dependent spin chain. These results indicate that the protocol we employ is also suited for larger systems, provided that the imbalance in spin the populations is respected. In the next sections, we show how increasing the number of internal components will affect the behavior of these quantities.

SU(3).
We now consider the case of a three-component strongly interacting fermionic gas with SU(3) symmetry. These systems are particularly interesting due to their connections to the quark model in the framework of quantum chromodynamics. We label the three internal states as |↑〉, |→〉 and |↓〉. While the Hamiltonian can still be described by Eq. (2), the permutation operator is now given by where λ is the vector composed by the eight generators of the SU(3) group, namely the Gell-Mann matrices. A system described by Eq. (7) can be mapped into the Lai-Sutherland model 62,63 which is a particular case of the spin-1 bilinear biquadratic model 64,65 . We keep the number of particles fixed as N = 4, with N ↑ = 2, N → = 1 and N ↓ = 1. The quench protocol and the quantities considered are the same as in the previous section.
In Fig. 5(a), we show 〈x(t) 2 〉 for each component and for the total density. Since we have two minority particles, each interacting with the remaining atoms with interaction strength g, the results for each of these components are identical. The excitation peaks seen in Fig. 5(b-d) reveal the contributions of the density and spin oscillations to the dynamics of each component (naturally, since we have N → = N ↓ = 1 plots (c) and (d) show identical results). Still, we can see that the majority component has a larger contribution to the total density excitations. The minority cases, however, show a slight increase in these frequencies as compared to the two-component case, with the spin oscillations remaining dominant.
An interesting perspective when dealing with multicomponent strongly interacting gases is a case where interactions are slightly imbalanced and a particular symmetry is broken. Here, we analyze the three-component case with broken SU(3) symmetry. It is useful, in this context, to rewrite the SU(3) permutation operator in terms of raising and lowering operators. These are defined as T ± = (λ 1 ± iλ 2 )/2, V ± = (λ 4 ± iλ 5 )/2 and U ± = (λ 6 ± iλ 7 )/2, where once again λ i are the Gell-Mann matrices. We still focus the particular case of N = 4, with N ↑ = 2, N → = 1 and N ↓ = 1. www.nature.com/scientificreports www.nature.com/scientificreports/ Below we rewrite the permutation operator with these modifications, including an additional symmetry-breaking parameter 1/η. The inclusion of η above means we are explicitly breaking the symmetry of the system by changing the energy contribution of turning |↑〉 into |→〉 and vice-versa. In Fig. 5(e-h) we show the result of breaking the SU(3) symmetry (by making η = 0.5) on the dynamics. While the effects in the |↑〉 and |→〉 components are subtle -an increase in the spin excitation peak as seen in panels (f) and (g) -in |↓〉 it is more drastic, with the spin contributions being distributed over different low frequencies. The remaining components still preserve isolated peaks for spin oscillations. These results point to the possibility of measuring the independence of spin and density dynamics even in a context where internal symmetries are not perfectly preserved. Naturally, different outcomes for the spin excitations can be expected by choosing a different value of η, or by breaking the symmetry in a different interaction channel.

SU(4).
We now examine the effect of applying our formalism to the case where the number of particles N matches the number of internal components. To that end, we consider the SU(4) fermionic gas with N = 4 and internal states labeled as |↑〉, |↗〉, |↘〉 and |↓〉. The number of particles in each state is thus given by (the so-called 1 + 1 + 1 + 1 infinitely repulsive system with different masses is known to have interesting properties, which were described in 66 ). We rewrite the permutation operator for the SU(4) system as where now λ represents the vector spanning the 15 SU(4) generators 67 . In the following, we focus on describing the results only for the |↑〉 and |↗〉 components.
It becomes clear that the behavior 〈x 2 (t)〉 as a function of time is the same for both components shown in Fig. 6 (this also holds for the other two components not shown). This is expected since the number of internal components matches the total number of particles in the system. Moreover, the frequency spectrum shows that the only contributions in the oscillations stem from the density excitations, as opposed to the previous cases. This allows us to interpret the dynamics of the SU(N) system with strong interactions as the one expected for a gas of impenetrable bosons, as long as the number of particles matches the number of internal components. This conclusion is in agreement with the observation that breathing mode frequencies of the SU(N) Fermi gas with strong repulsion approaches that of a Tonks-Girardeau gas 46 .
In the models considered here, a vanishing spin signal in the spectrum can additionally be obtained by taking a balanced system with a lower number of internal components (e.g. a 2 + 2 SU(2) system). This can be explained as a result of the symmetric perturbation to the potential that determines the initial state of the system. Turning off the barrier in this particular case has no effect on the ground state of the spin chain, which remains unchanged www.nature.com/scientificreports www.nature.com/scientificreports/ even as the density oscillations take place. However, for a matching number of particles and internal components, the results described in this section are the only possible outcomes for a system initialized in the ground state.

Discussion
We have presented an analysis of the dynamics of density and spin fluctuations in trapped few-body systems with SU(N) symmetry. In the limit of strong interactions, the Hamiltonian of this system can be mapped into a spin chain. This mapping allows for a straightforward generalization to systems with more internal components. In addition, it takes into account the geometry of the trapping potential into the set of exchange coefficients of the spin chain.
The dynamics of the system is obtained after a quench in the trapping potential, where a Gaussian barrier in the center of the harmonic trap is suddenly switched off. This simple protocol is particularly interesting from an experimental point of view, since it requires only minor modifications to the potential, without resorting to spin-selective traps. The sudden change induces the motion of the spatial degrees of freedom, which in turn are reflected in a time-dependence of the exchange coefficients of the spin chain. It is important to point out that, since the system is initialized in the ground state of the spin chain, the motion observed in this sector is only possible due to the quench in the potential. By monitoring the time evolution of the breathing modes given by the oscillations in 〈x 2 (t)〉, we describe the excitation spectrum of SU(2), SU(3) and SU(4) systems. Moreover, in the particular case where the number of internal components exactly matches the number of particles, we see that the spin excitations are completely washed out, and the only contributions are due to density oscillations that agree with those of a spinless Fermi gas. While the theoretical scheme presented here may in principle be hard to generalize to a large number of particles (specially due to the difficulties in calculating the dynamical exchange coefficients over long periods of time), a numerical many-body approach employing the same dynamical protocol may be able to do so. Aside details regarding the true strength of the interactions between non-identical particles, we expect a similar behavior to arise in this context, provided that the same essential ingredients are maintained.
Our results indicate that the decoupling of spin and density dynamics, rather than being exclusively a bulk effect in many-body ensembles, can occur in few-body systems under fairly simple conditions. We point out that multicomponent fermionic systems are the object of current experimental investigation, and the scheme presented in this work could serve as a starting point for experiments with many-body ensembles. The formalism we explored can also be used to predict the behavior of quantum gases with different atomic species (e.g. a bosonic mixture), or generalized to more involved quench protocols, simply by mapping the time evolution of the spatial orbitals into the exchange coefficients of the spin chain under consideration.  where the spinless fermion wave function Φ(x 1 , ..., x i , ..., x N ) is built as the Slater determinant of the N lowest occupied orbitals of the potential V(x). The corresponding energy of this state is given by E 0 in Eq. (2), which is calculated as the sum of the energies of each occupied level. This wave function can be simply viewed as the antisymmetrized version of the Tonks-Girardeau wave function for impenetrable bosons 68 . Methods for obtaining the exchange coefficients in different trapping geometries are available and can efficiently calculate α i for systems with up to N ~ 30 69,70 .

Methods
The spin densities given by Eq. (3) require calculating the spatial distribution of each individual particle in the trap, which is given by The expressions above involve multidimensional integrals of the spinless fermion wave function Φ(x 1 , ..., x N ), which for large systems can become increasingly hard to calculate. Moreover, if we consider time-dependent orbitals φ j (x, t), these expressions must be evaluated for all times after the sudden change in the trapping potential. Therefore, we explore the determinant form of the spinless fermion wave function Φ(x 1 , ..., x N ) to rewrite these equations in a shape which is suited to calculations in a dynamical context. Equivalent forms of these expressions have been presented in 53,69,70 . The individual one-body densities can be written as , and the subscript jk on the right side of Eq. (13) indicates that the j-th row and k-th column are removed. Although we omit the time in these expressions, we consider the single-particle orbitals φ i (x) to be time-dependent. In Eq. (3), we calculate the spin densities by combining the single-particle spatial distributions with the probability of finding a given spin component at each site. This last quantity can be calculated, for a given spin state ψ Details of the quench protocol. The trapping potential at t = 0 is given by t x s 0 0 2 2 0 ( / ) 2 where V 0 determines the height of the Gaussian peak and s sets its width. The system is therefore separated in an effective double-well by taking ω 0 = 1, V 0 = 25 and s = 0.1. The initial spinless fermion wave function Φ(x 1 , ..., x N , t = 0) is constructed with the single particle orbitals obtained by numerical diagonalization, using the N s = 35 lowest energy states of the harmonic oscillator. We note that, since the Gaussian peak is large compared to the individual densities of the orbitals, the ground state is quasi-degenerate (the two lowest energy states have nearly the same distribution, with opposite parity).
For t > 0, we make V 0 = 0 and the spinless fermion wave function is constructed by considering the time-evolved single particles orbitals φ j (x, 0) according to Floquet quasienergy spectrum for a time dependent hamiltonian. To obtain the frequency of spin oscillations for the time-dependent spin chain, we focus on finding the Floquet modes for the periodic Hamiltonian and calculating the gaps in the Floquet quasienergy spectrum 71 . Generally speaking, we are interested in finding the time evolution operator at any time t by solving the following differential equation for U(t, t 0 ): By obtaining and diagonalizing U(T, 0) (where T is the period of the Hamiltonian) we obtain a set of quasienergies ε n and Floquet modes at t = 0, which we write as |u n (0)〉. The time evolution of the system is then given by n n n n where |u n (t)〉 denotes the time-evolved Floquet modes and c n = 〈u n (0)|ψ(0)〉. Once we have calculated ε n , we can easily obtain the energy gaps that contribute to the time evolution of the system by considering the dominant contributions given by c n . In a small system, this is a fairly straightforward process, since the initial state projects only onto a few Floquet modes.