Probing ferromagnetic order in few-fermion correlated spin-flip dynamics

Magnetism constitutes a principal feature of a large class of materials and represents a macroscopic phenomenon of quantum origin [1–3]. In conductors the magnetic properties of the delocalized (itinerant) electrons are qualitatively understood in terms of the Stoner instability [4]. To verify and emulate the latter mechanism ultracold fermionic ensembles have been employed [5–7]. However, the nature of the interparticle interaction exhibited in three-dimensional ultracold gases hindered the study of itinerant ferromagnetism as the repulsive Fermi gas is metastable due to bosonic Feshbach molecule formation [8]. Utilizing fast interaction quenches, it has been shown that no ferromagnetic phase can be achieved as the decay into molecules is faster than the formation of ferromagnetic domains [6, 7]. Instead, recent pump-probe experiments [9] indicate that the formation rate of ferromagetic domains with a size comparable to the interatomic separation is larger than the corresponding molecular decay rate. Furthermore, ferromagnetic properties have been observed indirectly either by the spectroscopic study of strongly particle-imbalanced [10, 11] (supplemented by [12]) and particle-balanced [9] two-component Fermi mixtures or by employing a binary Fermi gas prepared in a magnetic domain wall structure [13]. The latter experimental evidence poses the question whether stable ferromagnetism can be observed in the absence of molecule formation.

A controllable setting that can shed light on such inquiries is the experimentally accessible few-fermion quasi-one dimensional (1D) gas [14]. Owing to its 1D character, a shallow two-body bound state for effectively repulsive interactions is absent and thus the molecule formation is suppressed. Moreover, the experimental [15, 16] and theoretical [17–21] study of the magnetic properties of few-fermion systems has led to the insight that for near zero and infinite interactions there is an approximate mapping of the 1D spin-1/2 Fermi gas to an effective spin-chain model [22–28]. Most importantly, these spinor systems possess experimentally accessible eigenstates of ferromagnetic nature⁴ , namely the interaction-independent spin-polarized states. Consequently, the study of the dynamical stability of these ferromagnetic states is essential for our understanding of the magnetic properties of 1D systems. The study of the mechanisms emanating in 1D systems might in turn provide important insights for higher dimensional settings.

We study the dynamical stability of the fully polarized 1D parabolically-confined few-fermion spin-1/2 gas under the effect of inhomogeneous Rabi coupling of the spin-states. This coupling scheme introduces a spatially dependent spin-flip transition amplitude and thus probes the stability of the initial state by breaking the spin-symmetries of the unperturbed system (see below). An argumentation based on the Hartree–Fock (HF) framework of the Stoner model testified within the time-dependent Hartree–Fock (TDHF) [29] showed that the spin-polarization of the Fermi gas is stable for interparticle repulsions that exceed the kinetic and spin-flip contributions [30, 31]. Inspecting the correlated spin-flip dynamics within the latter interaction regime (where TDHF predicts stable ferromagnetism) we reveal that the many-body (MB) spin–spin correlator exhibits ferromagnetic spin–spin correlations throughout the dynamics. Moreover and in contrast to the TDHF results the MB state of the ferromagnetically correlated gas shows an unstable polarization fluctuating between fully polarized and almost completely unpolarized. This outcome cannot be retrieved within the HF description and exposes the crucial role of correlations in the magnetic properties of spin-1/2 fermions even away from the strongly interacting regime. We show that the decay of the polarization and the emerging correlated spin-order can be understood by generalizing the spin-chain model of [28]. The coupling of the initial state to lower spin-S values is found to cause the dephasing of the collective Larmor precession of the spins. This dephasing dynamically leads the system to an almost equal superposition of the two ferromagnetic fully polarized states of opposite spin-orientation i.e. a NOON-like state [32, 33] with zero total polarization. For weaker and stronger interactions lying outside the above-mentioned ferromagnetically ordered regime the system undergoes a demagnetization dynamics which is identified and characterized. Our results generalize to other particle numbers within the few-body regime. The employed setup can be implemented in state-of-the-art ${}^{40}$K experiments and the corresponding findings can be probed by fluorescence imaging techniques. Additionally, we showcase that our findings can be generalized to a broader class of dynamical scenarios characterized by different initial states and Rabi-coupling potentials. We explicitly demonstrate the robustness of our results against common noise sources exhibited in such experiments by performing simulations of single-shot images.

The presentation of our results proceeds as follows. In section 2 we discuss the setup and the basic observables used for the interpretation of the spin-dynamics. In section 3 we present our results for the correlated spin-flip dynamics employing the Multi-Configuration Time-Dependent Method for Fermions (MCTDHF) [34–39] and interpret them in terms of two approximate methods for the case N = 3. Section 4 provides a generalization of our findings to the case of N = 5 fermions. In section 5 we also demonstrate that the observed dynamical phenomena persist for different initial states and inhomogeneous Rabi-couplings. A possible experimental probe of our predictions and its feasibility are discussed in section 6. In section 7 we summarize our results and provide an outlook. Appendix A addresses our numerical methodology based on MCTDHF. The numerical implementation of the single-shot simulations for spinor fermions is briefly discussed in appendix B. Finally, in appendix C we derive the effective spin-chain model for our system.

We consider an interacting system of N spin-1/2 fermions of mass m, confined in an 1D parabolic trap of frequency ω. The interparticle interaction emanating in such systems is well-described by the s-wave contact interaction of strength, g [40]. The latter can be manipulated by exploiting either Fano-Feshbach or confinement induced resonances [8]. The MB Hamiltonian that models such a system reads $\hat{\tilde{H}}={\hat{H}}_{0}+{\hat{H}}_{I}$, where the single-particle Hamiltonian ${\hat{H}}_{0}$ is

$$\begin{eqnarray}&&{\hat{H}}_{0}=\displaystyle \sum _{\alpha \in \{\uparrow ,\downarrow \}}\int {\rm{d}}x\,{\hat{\psi }}_{\alpha }^{\dagger }(x)\left(-\displaystyle \frac{{{\hslash }}^{2}}{2m}\displaystyle \frac{{{\rm{d}}}^{2}}{{\rm{d}}{x}^{2}}+\displaystyle \frac{1}{2}m{\omega }^{2}{x}^{2}\right){\hat{\psi }}_{\alpha }(x).\end{eqnarray} \tag{ 1 }$$

The corresponding interparticle interaction term is encoded by

$$\begin{eqnarray}&&{\hat{H}}_{I}=g\int {\rm{d}}x\,{\hat{\psi }}_{\downarrow }^{\dagger }(x){\hat{\psi }}_{\uparrow }^{\dagger }(x){\hat{\psi }}_{\uparrow }(x){\hat{\psi }}_{\downarrow }(x),\end{eqnarray} \tag{ 2 }$$

where ${\hat{\psi }}_{\alpha }(x)$ denotes the fermionic field-operator with spin $\alpha \in \{\uparrow ,\downarrow \}$. The Hamiltonian $\hat{\tilde{H}}$ exhibits several crucial spin-symmetries. It can be shown that $\hat{\tilde{H}}$ commutes with each component of the total spin-vector operator

$$\begin{eqnarray}&&\hat{{\boldsymbol{S}}}=\displaystyle \frac{{\hslash }}{2}\displaystyle \sum _{\alpha \alpha ^{\prime} }\int {\rm{d}}x\,{\hat{\psi }}_{\alpha }^{\dagger }(x){{\boldsymbol{\sigma }}}_{\alpha ,\alpha ^{\prime} }{\hat{\psi }}_{\alpha ^{\prime} }(x),\end{eqnarray} \tag{ 3 }$$

where ${\boldsymbol{\sigma }}$ denotes the Pauli vector. Additionally, it possesses a SU(2) symmetry stemming from its commutation with the total spin-magnitude operator

$$\begin{eqnarray}&&\hat{{S}^{2}}=\displaystyle \frac{3N{{\hslash }}^{2}}{4}+\mathop{\underbrace{\displaystyle \frac{{{\hslash }}^{2}}{4}{\displaystyle \sum }_{\alpha ,\alpha ^{\prime} ,\beta ,\beta ^{\prime} }\int {{\rm{d}}}^{2}x\,\left({{\boldsymbol{\sigma }}}_{\alpha ,\alpha ^{\prime} }\cdot {{\boldsymbol{\sigma }}}_{\beta ,\beta ^{\prime} }\right){\hat{\psi }}_{\alpha }^{\dagger }({x}_{1}){\hat{\psi }}_{\beta }^{\dagger }({x}_{2}){\hat{\psi }}_{\beta ^{\prime} }({x}_{2}){\hat{\psi }}_{\alpha ^{\prime} }({x}_{1})}}\limits_{={\displaystyle \sum }_{i\ne j}{\hat{{\boldsymbol{S}}}}_{i}\cdot {\hat{{\boldsymbol{S}}}}_{j}}.\end{eqnarray} \tag{ 4 }$$

These symmetries imply that the eigenvalues of ${\hat{S}}_{z}$ and ${\hat{S}}^{2}$ define good quantum numbers. Consequently, the ferromagnetic fully spin-polarized state, $| {{\rm{\Psi }}}_{F}\rangle ={\prod }_{i=0}^{N-1}\int {\rm{d}}x\,{\phi }_{i}(x){\hat{\psi }}_{\uparrow }^{\dagger }(x)| 0\rangle$, where ϕ_i(x) refers to the ith eigenfunction of the 1D harmonic oscillator, is the energetically lowest eigenstate of $\hat{\tilde{H}}$ (note that ${\hat{H}}_{I}| {{\rm{\Psi }}}_{F}\rangle =0$) with total spin eigenvalues $S={S}_{z}=\tfrac{N}{2}$.

To controllably probe the stability of such a ferromagnetic state $| {\rm{\Psi }}(0)\rangle =| {{\rm{\Psi }}}_{F}\rangle$ under coherent processes that break both ${\hat{S}}_{z}$ and ${\hat{S}}^{2}$ symmetries we employ an inhomogeneous Rabi coupling between the spin $\uparrow$ and $\downarrow$ states. Note here that similar Rabi-coupling techniques have been employed in several experiments e.g. see [41–44] involving binary bosonic mixtures. The resulting Hamiltonian of the total system reads

$$\begin{eqnarray}&&\hat{H}={\hat{H}}_{0}+{\hat{H}}_{S}+{\hat{H}}_{I},\end{eqnarray} \tag{ 5 }$$

where the externally-imposed Rabi coupling term is

$$\begin{eqnarray}&&{\hat{H}}_{S}={B}_{0}\displaystyle \sum _{\alpha ,\alpha ^{\prime} }\int {\rm{d}}{x}\,{{\rm{e}}}^{-\tfrac{{x}^{2}}{2{w}^{2}}}{\hat{\psi }}_{\alpha }^{\dagger }(x){\sigma }_{\alpha ,\alpha ^{\prime} }^{x}{\hat{\psi }}_{\alpha ^{\prime} }(x).\end{eqnarray} \tag{ 6 }$$

In particular, ${\hat{H}}_{S}$ induces spin-flip transitions with a spatially dependent transition amplitude, modelled by a Gaussian of width w and intensity B₀. This coupling scheme can be realized in ultracold experiments by optical Raman dressing of the two lowest hyperfine levels of ${}^{40}{\rm{K}}$ (see also section 4). We choose the values w = 2 and ${B}_{0}=2.5/\sqrt{8\pi }$ (in harmonic oscillator units) leading to an average precession (Larmor) frequency, ${\omega }_{L}\approx 0.85$, for the spins which is lower than all collective mode frequencies (e.g. ω_R ≈ 2 for the breathing mode). This choice enables us to avoid spin segregation phenomena [28] occurring when the length scale of the modulation w is smaller than that of the trap ${l}_{\omega }=\sqrt{{\hslash }/m\omega }$.

Our goal is to inspect the stability of the ferromagnetism when the ${\hat{H}}_{S}$ term, equation (6), is abruptly switched on at t = 0. To achieve this we track two main observables, directly related with the system's broken symmetries. Namely, the normalized spin polarization magnitude ${P}_{S}^{(1)}=\tfrac{2}{{\hslash }N}| \langle \hat{{\boldsymbol{S}}}\rangle |$ and the spin–spin correlator ${C}_{S}^{(2)}=\tfrac{4\langle {\hat{S}}^{2}\rangle -3N}{{{\hslash }}^{2}N(N-1)}=\tfrac{{\sum }_{i\ne j}\langle {\hat{{\boldsymbol{S}}}}_{i}\cdot {\hat{{\boldsymbol{S}}}}_{j}\rangle }{{{\hslash }}^{2}N(N-1)}$. P⁽¹⁾_S expresses the averaged spin-order (magnetization) and refers to the magnitude of the polarization. Due to its one-body character P⁽¹⁾_S does not probe the correlations that might emerge in the system despite being affected by them. For this purpose we employ the spin–spin correlator, C⁽²⁾_S, which probes the alignment of each two spins and serves as an indicator for the distinction of ferromagnetic C⁽²⁾_S ≈ 1, antiferromagnetic ${C}_{S}^{(2)}$ ≈ −1 and paramagnetic C⁽²⁾_S ≈ 0 spin–spin correlations.

3.1. Many-body correlated spin-flip dynamics

Figures 1(a) and (b) present our MB results for the paradigmatic case of N = 3 fermions obtained via MCTDHF, that enables us to capture all interparticle correlation effects [34–39]. The MCTDHF method is a variational, numerically exact, ab initio method for solving the time-dependent MB Schrödinger equation which includes all correlation effects. It is based on expanding the MB wavefunction in terms of a time-dependent variationally optimized basis. In this way it enables us to efficiently truncate the MB Hilbert space relevant for our system by using a computationally feasible basis size. The MCTDHF method exhibits increased numerical efficiency when compared to an expansion relying on a time-independent basis since the number of basis states can be significantly reduced. A detailed discussion on the capabilities and the MB wavefunction ansatz of the above-mentioned method is presented in appendix A.

Zoom In Zoom Out Reset image size

For weak repulsive or attractive interactions, $| g| \lt 0.5$, a rapid demagnetization (see the decaying behaviour of P_S⁽¹⁾) is observed, accompanied by a loss of the spin alignment, C⁽²⁾_S, at a time scale ∼20. Partial revivals of both P_S⁽¹⁾ and C⁽²⁾_S appear over regular time intervals for later times. Our results for this interaction interval are compliant with the spin-dynamics analysed in [28, 45] and we shall refer to this regime as the weak-g demagnetization regime. Indeed within this regime each of the particles precesses with a different Larmor frequency leading to the loss of the polarization magnitude P⁽¹⁾_S after a few precession cycles. For intermediate interactions of either sign, i.e. $0.5\lt | g| \lt 4$, the decay and revival of ${P}_{S}^{(1)}$ also occurs but at a drastically increased time-scale (which increases further for larger $| g|$) when compared to the weak-g demagnetization regime. In contrast, C⁽²⁾_S indicates that the spins are close to be maximally aligned (e.g. C⁽²⁾_S ≥ 0.85 for 2 < g < 4 and C⁽²⁾_S ≥ 0.95 for g ∼ −2 in figure 1(b)) throughout the evolution, signifying ferromagnetic spin-correlations. Therefore, ferromagnetism is unstable in this interaction interval as the polarization (P⁽¹⁾_S) of the ensemble features large fluctuations despite the ferromagnetic order captured by C⁽²⁾_S which is almost perfect. Hence, we refer to this regime as ferromagnetically ordered. In particular, it involves a different spin-order than ferromagnetism, as its order is inferred by the ferromagnetic spin–spin correlations rather than the polarization. For g > 4 a suppression of the ferromagnetic spin–spin correlations occurs as the amplitude of the C⁽²⁾_S oscillations increases for stronger g, see figure 1(b). For instance, at g ≈ 10, C⁽²⁾_S fluctuates between the values 0.5 and unity. P⁽¹⁾_S is also oscillatory taking values between unity and 1/3, with a significantly smaller oscillation frequency than ${C}_{S}^{(2)}$ (see figure 1(a)). In the following this interaction interval (g > 4) is referred to as the strong-g demagnetization regime.

3.2. Spin-flip dynamics within TDHF

To demonstrate the crucial role of correlations within the MB dynamics we compare the above MB findings with the TDHF approximation presented in figures 1(c) and (d). For weak $| g|$ ($| g| \lt 0.5$) the demagnetization dynamics is qualitatively captured by the TDHF approximation. However, upon increasing $| g|$, (in particular for $0.5\lt | g| \lt 4$) TDHF predicts no loss of ${P}_{S}^{(1)}$ in contrast to the MB case (compare figures 1(a) and (c)), while a similar spin-correlation dynamics is observed compare figures 1(b) and 1(d). This behaviour of the TDHF can be interpreted in terms of the Stoner model [4, 30], see figure 2(a). Indeed, within HF the interaction energy of contact interacting spin-1/2 fermions is proportional to the density overlap between the two spin-components [4, 30]. Therefore, for large enough g the system initialized in a spin-polarized state characterized by zero interaction energy, cannot access states with a single (or more) spin-flips due to their large interaction energy. Thus, each of the spins has to preccess with the same frequency resulting in the constant polarization magnitude, P⁽¹⁾_S. For strong g > 4, Rabi oscillations between the ferromagnetic initial state $| {{\rm{\Psi }}}_{F}\rangle$ (characterized by P⁽¹⁾_S = 1, C⁽²⁾_S = 1) and the paired HF ground states $| {{\rm{\Psi }}}_{\alpha }^{\mathrm{HF}}\rangle \sim \int {{\rm{d}}}^{3}x{\phi }_{0}({x}_{1}){\phi }_{0}({x}_{2}){\phi }_{1}({x}_{3}){\psi }_{\uparrow }^{\dagger }({x}_{1}){\psi }_{\downarrow }^{\dagger }({x}_{2}){\psi }_{\alpha }^{\dagger }({x}_{3})| 0\rangle$ (referring to P⁽¹⁾_S = 1/3, C⁽²⁾_S = − 1/2) take place. This indicates that the interparticle repulsion between the paired fermions exactly balances the energy benefit of their pairing in the same state and corresponds to the Stoner instability of the ground state. It is important to note here that these Rabi oscillations are absent in the MB case, see also figures 1(c) and (d) for g ≈ 5.

Zoom In Zoom Out Reset image size

Concluding, the ferromagnetically ordered regime exhibited in the MB case corresponds to a stable ferromagnetic one within the HF framework. This observation exposes the correlated nature of the ferromagnetically ordered regime. In both cases the interaction regime is limited to intermediate values of g and in particular to $0.5\lt | g| \lt 4$ (see figures 1(b), (d)). However, for g > 4 the mechanism that breaks the spin-order differs. In the HF case the ground state Stoner instability takes place which is forbidden for any finite repulsive interaction in the MB case [46, 47]. Instead, it is known that the unpaired states of maximum S (see figure 1(b)) and the MB antiferromagnetic ground states exhibit a crossing in the Tonks-Girardeau limit [16–21]. In our case due to the breaking of the SU(2) symmetry the above consist an avoided crossing [27, 48, 49] which is approached for increasing g > 4. As we shall argue in section 3.4, the fluctuations of C⁽²⁾_S in the strong-g demagnetization regime for the MB case can be attributed to this avoided crossing.

3.3. Effective spin-chain model

3.4. Analysis of the microscopic mechanisms

To identify the underlying mechanisms of the MB spin-dynamics we invoke the spectrum of the spin-polarization, namely ${ \mathcal F }\{{P}_{S}^{(1)}\}({\rm{\Omega }})=| \int {\rm{d}}t\,{{\rm{e}}}^{{\rm{i}}{\rm{\Omega }}t}{P}_{S}^{(1)}(t){| }^{2}$, presented in figures 3(a) and (b) for −2 < g < 10 and $| g| \lt 1$ respectively. Recall that each branch in the spectrum of P_S⁽¹⁾ corresponds to an energy difference between two eigenenergies of $\hat{H}$, (see equation (5)). For g = 0 three distinct Larmor frequencies⁵ occur that correspond to the three energy differences among the occupied single-particle eigenstates in the spatial configuration $\vec{n}=(0,1,2)$. For g $\ne$ 0 a multitude of interaction-dependent frequency branches emerge from each Larmor frequency. The failure of TDHF to capture even on the qualitative level the spin dynamics even for low g is evident in ${ \mathcal F }\{{P}_{S}^{(1)}\}({\rm{\Omega }})$. Indeed, the TDHF captures only one frequency per particle for $| g| \lt 0.1$ (see figure 3(c)) in contrast to the multitude of interaction-dependent frequency branches emerging from each Larmor frequency in the correlated case (see figure 3(b)). Thus we can conclude that the build up of correlations in the MB case even for very small g completely invalidates the HF picture for the spin dynamics. Such correlations are of particular importance in the ferromagnetically ordered regime. In this case, three dominant branches appear in the vicinity of Ω ≈ 0.85 (see also figure 3(b) for g ≈ 1) within the correlated case that lead to the beating dynamics of P⁽¹⁾_S, observed in figures 1(a), (b).

Zoom In Zoom Out Reset image size

The origin of the above-mentioned frequencies can be exposed by comparing ${ \mathcal F }\{{P}_{S}^{(1)}\}$, with the energy differences of the eigenstates of the spin-chain model. As anticipated by our discussion in section 3.3 the eigenstates of the configuration $\vec{n}=(0,1,2)$ are expected to well capture the P⁽¹⁾_S dynamics. The eigenspectrum of the spin-chain Hamiltonian ${\hat{H}}_{\vec{n}}={\hat{H}}_{(\mathrm{0,1,2})}$ is presented in figure 3(d). The spin-chain eigenstates within the weak-g demagnetization regime ($| g| \lt 0.5$) are ordered in terms of increasing S_x due to the Zeeman effect induced by the effective magnetic field along the x axis that ${\hat{H}}_{S}$, equation (6), introduces. Within this weak-g demagnetization regime an avoided crossing between the distinct eigenstates of the same ${S}_{x}=\pm \tfrac{1}{2}$ occurs which can be attributed to the breaking of the ${\hat{S}}^{2}$ symmetry by the inhomogeneity of the ${\hat{H}}_{S}$ term. The states of the spin-chain model can be labelled as $| S,{S}_{x}{\rangle }^{{\ell }}$ in the repulsive (${\ell }=+$) and attractive (ℓ = −) part of the ferromagnetically ordered regime ($0.5\lt | g| \lt 4$). This is possible because S is an approximate good quantum number within this interaction range. This property can be identified by examining the expectation value of $\langle {\hat{S}}^{2}\rangle \equiv {}^{{\ell }}\langle S,{S}_{x}| {\hat{S}}^{2}| S,{S}_{x}{\rangle }^{{\ell }}$. Figure 3(e) presents this expectation value for the states with ${S}_{x}=\tfrac{1}{2}$ and varying interaction strength. It becomes evident that for increasing $| g|$ the $\langle {\hat{S}}^{2}\rangle$ of one of these states [i.e. ${\left|\tfrac{3}{2},\tfrac{1}{2}\right\rangle }^{+}$ and ${\left|\tfrac{3}{2},\tfrac{1}{2}\right\rangle }^{-}$ in the case of g > 0 and g < 0 respectively] approaches the value $\langle {\hat{S}}^{2} \rangle =\tfrac{3}{2}\left(\tfrac{3}{2}+1\right)$, indicating that $S\approx \tfrac{3}{2}$. Moreover, the states $\{| \tfrac{3}{2},\tfrac{1}{2}{\rangle }^{-},| \tfrac{3}{2},\tfrac{1}{2}{\rangle }^{\pm }\}$ and $\{| \tfrac{3}{2},\tfrac{1}{2}{\rangle }^{+},| \tfrac{3}{2},\tfrac{1}{2}{\rangle }^{\pm }\}$ in the case of g > 0 and g < 0 respectively tend to $\langle {\hat{S}}^{2}\rangle =\tfrac{1}{2}(\tfrac{1}{2}+1)$. We remark here that the states with ${S}_{x}=-\tfrac{1}{2}$ exhibit a similar behaviour as the aforementioned ${S}_{x}=\tfrac{1}{2}$ case and the eigenstates $| \tfrac{3}{2},\pm \tfrac{3}{2}{\rangle }^{\pm }$ correspond to fully polarized states along the x axis with $\langle {\hat{S}}^{2} \rangle =\tfrac{3}{2}\left(\tfrac{3}{2}+1\right)$ for all g (not shown here for brevity). By employing the spin-chain eigenspectrum, figure 3(d), we can identify all the energy branches appearing in figure 3(b) with the corresponding eigenenergy differences of the spin-chain model. Most importantly, figure 3(f) presents this identification within the ferromagnetically ordered regime ($0.5\lt | g| \lt 4$). As it can be seen the energy differences between the spin-chain eigenstates with $S\approx \tfrac{3}{2},{\rm{\Delta }}E\{| \tfrac{3}{2},{S}_{x}{\rangle }^{{\ell }},| \tfrac{3}{2},{S}_{x}^{{\prime} }{\rangle }^{{\ell }^{\prime} }\}$, possess the main contribution to ${ \mathcal F }\{{P}_{S}^{(1)}\}$ in the ferromagnetically ordered regime as they match well with the dominant energy branches appearing in the MB dynamics. Based on this identification we can draw several conclusions regarding the order exhibited within the ferromagnetically ordered regime. First, the ferromagnetic spin-correlations emanate from the predominant occupation of MB eigenstates with $S\approx 3/2$ each characterized by a different value of S_x. The frequency difference ΔΩ between the highest (${\rm{\Delta }}E\{| \tfrac{3}{2},\tfrac{3}{2}{ \rangle }^{\pm },| \tfrac{3}{2},\tfrac{1}{2}{ \rangle }^{\ell }$) and lowest lying (${\rm{\Delta }}E\{| \tfrac{3}{2},-\tfrac{1}{2}{ \rangle }^{\ell },| \tfrac{3}{2},-\tfrac{3}{2}{ \rangle }^{\pm }$) of the above-mentioned dominant branches (see also figure 3(d)) results in the dephasing of the initial superposition. In terms of the spin-$\uparrow$ and spin-$\downarrow$ states this means that the system performs spin-flips between the different states with $S=\tfrac{3}{2}$ and varying S_z leading to the decay of P⁽¹⁾_S, see figure 2(b). The corresponding timescale is ${\tau }_{D}\approx \tfrac{\pi }{{\rm{\Delta }}{\rm{\Omega }}}$. ${\rm{\Delta }}{\rm{\Omega }}$ is attributed to the energy shift of the $| \tfrac{3}{2},\pm \tfrac{1}{2}{\rangle }^{{\ell }}$ states from being equidistantly spaced which is induced by an avoided crossing at g = 0, as in figure 3(d). The increase of $| g|$ within the ferromagnetically ordered regime, leads to a decrease of ΔΩ (or equivalently increase of τ_D) giving rise to the observed P_S⁽¹⁾ dynamics (figures 1(a) and (e)).

Regarding the strong-g demagnetization regime (see figure 3(a) for g > 4) the frequency branches that correspond to the $| \tfrac{3}{2},\pm \tfrac{1}{2}{\rangle }^{+}$ states deviate from Ω ≈ 0.85 as they couple to the antiferromagnetic ground states of $S\approx \tfrac{1}{2}$ and ${S}_{x}=\pm \tfrac{1}{2}$ character. Due to the same mechanism additional branches also appear in ${ \mathcal F }\{{P}_{S}^{(1)}\}$ (see the corresponding arrows in figure 3(a)) that contribute to the complex P⁽¹⁾_S dynamics exhibited in this interaction regime (see figure 1(a)) and results in the oscillatory patterns of ${C}_{S}^{(2)}$ (see figure 1(b), and also our discussion in section 3.2).

3.5. Characterization of entanglement in the ferromagnetically ordered regime

The separation of the energy scales ${\rm{\Delta }}{\rm{\Omega }}\sim 0.05$ (see figure 3(c)) and the average Larmor frequency, ${{\rm{\Omega }}}_{L}=\tfrac{1}{3{\hslash }}{\rm{\Delta }}E\{| \tfrac{3}{2},\tfrac{3}{2}\rangle ,| \tfrac{3}{2},-\tfrac{3}{2}\rangle \}\sim 0.85$ observed in figure 3, enables us to further characterize the superpositions that emerge in the ferromagnetically ordered regime. To this end we introduce the precessing bases $| \tfrac{3}{2},{S}_{j^{\prime} (t)}\rangle ={{\rm{e}}}^{{\rm{i}}\tfrac{1}{2}{\hat{S}}_{x}{{\rm{\Omega }}}_{L}t}| \tfrac{3}{2},{S}_{j}\rangle$ ($j\in \{x,y,z\}$) and project the MB wavefunction obtained within MCTDHF to these basis states. Note that if all the particles were collectively precessing with the frequency Ω_L then $| \langle \tfrac{3}{2},{S}_{z^{\prime} (t)}=\tfrac{3}{2}| {\rm{\Psi }}(t)\rangle |$ would be constant in time. However, as we have already enstablished in section 3.2, this is not the case. Figure 4(a) presents the results of this projection for a representative case (g = 2) within the ferromagnetically ordered regime. We observe a low-frequency (∼ΔΩ) population transfer from the state $| \tfrac{3}{2},{S}_{z^{\prime} (t)}=\tfrac{3}{2}\rangle$ to the state $| \tfrac{3}{2},{S}_{z^{\prime} (t)}=-\tfrac{1}{2}\rangle$. For t > 350 the latter mechanism results in $\langle {\hat{S}}_{z^{\prime} (t)}\rangle \approx 0$, as $| \tfrac{3}{2},{S}_{z^{\prime} (t)}=-\tfrac{1}{2}\rangle$ possesses approximatively a three times larger population as compared to $| \tfrac{3}{2},{S}_{z^{\prime} (t)}=\tfrac{3}{2}\rangle$. The nature of this superposition can be understood by transforming to the orthogonal precessing axis, y'(t) (see figure 4(b)). In this case the populations of $| \tfrac{3}{2},{S}_{y^{\prime} (t)}=\tfrac{3}{2}\rangle$ and $| \tfrac{3}{2},{S}_{y^{\prime} (t)}=-\tfrac{3}{2}\rangle$ are almost equal for t > 350, signifying the tendency to dynamically approach a NOON state, characterized by $\langle {\hat{S}}_{y^{\prime} (t)}\rangle \approx 0$, i.e. $| {\rm{\Psi }}({t}_{0})\rangle \,\sim \tfrac{1}{\sqrt{2}}(| \tfrac{3}{2},{S}_{y^{\prime} (t)}=\tfrac{3}{2}\rangle +{{\rm{e}}}^{{\rm{i}}\phi }| \tfrac{3}{2},{S}_{y^{\prime} (t)}=-\tfrac{3}{2}\rangle )$ with a relative phase ϕ. These results combined with the conserved quantity $\langle {\hat{S}}_{x}\rangle =0$ explain the decay of the total magnetization P⁽¹⁾_S → 0 for t > 350. Accordingly, the spin-dynamics within the ferromagnetically ordered regime describes the dynamical evolution of a fully polarized state to a superposition one consisting of two antiparallel-oriented fully polarized states.

Zoom In Zoom Out Reset image size

We next demonstrate the robustness of our main findings for the case of larger particle numbers by examining a system consisting of N = 5 fermions.

Figures 5(a)–(f) present P⁽¹⁾_S and C⁽²⁾_S for N = 5 fermions within the three different approaches employed above, i.e. the MCTDHF, TDHF and spin-chain approach. A similar spin-dynamics as for the N = 3 case is observed for both quantities, but in different interaction regimes, caused by the increase of N. The ferromagnetically ordered regime occurs in the range $0.5\lt | g| \lt 2$, where both P⁽¹⁾_S and ${C}_{S}^{(2)}$ possess a value proximal to unity within the TDHF approach (see figures 5(c) and (d)). In the same range the MB treatment provided by MCTDHF reveals that P⁽¹⁾_S is decaying (see figure 5(a)), a feature which is also well captured by the spin-chain method (see figure 5(e)). The only additional structures that emerge in the MB spin-flip dynamics when compared to the N = 3 case are very narrow interaction windows where C⁽²⁾_S gets significantly depleted from unity (see figure 5(b), $g\approx 2.5$). These regions can be attributed to avoided crossings between the different spin-states of the dominant spatial configuration $\vec{n}=(0,1,2,3,4)$ with states characterized by spatial configurations with double occupations of single-particle spatial mode(s) [e.g. $\vec{n}^{\prime} =(0,1,2,2,3)$].

Zoom In Zoom Out Reset image size

Inspecting ${ \mathcal F }\{{P}_{S}^{(1)}\}$ for N = 5 fermions, see figure 5(g), similar microscopic mechanisms to the N = 3 case can be observed in both the weak-g and the ferromagnetically ordered regime. Despite the fact that more states are involved, the main features essentially remain the same. The weak-g demagnetization regime originates from the multitude of branches emerging from the five available non-interacting Larmor frequencies. However, only five (which can be identified as the energy differences between the $| \tfrac{5}{2},{S}_{x}{ \rangle }^{\ell }$ states) possess a significant amplitude for $| g| \gt 0.5$. The frequency difference, ΔΩ, between the highest and lowest lying of the above five branches (see figure 5(g)) gives rise to the decay of P_S⁽¹⁾ within the ferromagnetically ordered regime observed in figure 5(a).

Finally, we show that even the superpositions emerging in the dynamics are of the same character as for N = 3 fermions. To reveal this we construct the precessing basis analogously to the case N = 3, namely $| \tfrac{5}{2},{S}_{j}^{{\prime} }(t)\rangle ={{\rm{e}}}^{{\rm{i}}\tfrac{1}{2}{\hat{S}}_{x}{{\rm{\Omega }}}_{L}t}| \tfrac{5}{2},{S}_{j}\rangle$, ${{\rm{\Omega }}}_{L}=\tfrac{1}{5{\hslash }}{\rm{\Delta }}E\{| \tfrac{5}{2},\tfrac{5}{2}\rangle ,| \tfrac{5}{2},-\tfrac{5}{2}\rangle \}$ and expand the MB wavefunction in terms of the latter. Figures 5(h) and (i) present the results of this expansion for the axes z'(t) and y'(t) respectively and for g = 1 within the ferromagnetically ordered regime. Figure 5(h) demonstrates that the collective precession of the spins characterized by $| \langle \tfrac{5}{2},{S}_{z^{\prime} (t)}=\tfrac{5}{2}| {\rm{\Psi }}(t)\rangle {| }^{2}=1$ gets quickly dephased. At later times t ≈ 350 the dephasing of the collective Larmor precession leads to the formation of a NOON-like state characterized by $| \langle \tfrac{5}{2},{S}_{y^{\prime} (t)}=\tfrac{5}{2}| {\rm{\Psi }}(t)\rangle {| }^{2}=| \langle \tfrac{5}{2},{S}_{y^{\prime} (t)}=-\tfrac{5}{2}| {\rm{\Psi }}(t)\rangle {| }^{2}\approx 1/2$ (see figure 5(i)), compliant with our N = 3 results.

Below we demonstrate that the above identified ferromagnetic properties are not restricted to the previously examined out-of-equilibrium scenario. Indeed, we will show that the 1D spin-1/2 Fermi gas exhibits a similar spin-dynamics for different initial states characterized by $S\approx \tfrac{N}{2}$ within the ferromagnetic-like regime, which, furthermore, does not depend on the exact form of the Rabi coupling potential⁶ . The special feature of the specific dynamical protocol investigated in the previous sections is that it can be readily implemented in state-of-the-art experiments (see also section 6).

It can be shown that any initial state which is dominated by the states of the $\vec{n}=(0,\,\ldots ,\,N-1)$ spatial configuration with $S=\tfrac{N}{2}$ possesses a similar spin–spin correlation dynamics to the ferromagnetic one, $| {{\rm{\Psi }}}_{F}\rangle ={\prod }_{i=0}^{2}\int {\rm{d}}x\,{\phi }_{i}(x){\hat{\psi }}_{\uparrow }^{\dagger }(x)| 0\rangle$. The reason is that the time-evolution of the expectation value of ${\hat{S}}^{2}$ reads

$$\begin{eqnarray} &&\langle {\rm{\Psi }}(t)| {\hat{S}}^{2}| {\rm{\Psi }}(t) \rangle =\displaystyle \sum _{{S}_{x}=-\tfrac{N}{2}}^{\tfrac{N}{2}}{\left| \langle {\rm{\Psi }}(t)\left|\displaystyle \frac{3}{2};{S}_{x}\right. \rangle | \right.}^{2}\left.\Space{0ex}{0.50ex}{0ex}\displaystyle \frac{3}{2};{S}_{x}\right|{{\rm{e}}}^{\tfrac{{\rm{i}}t}{\hslash }\hat{H}}{\hat{S}}^{2}{{\rm{e}}}^{-\tfrac{{\rm{i}}t}{\hslash }\hat{H}}\left|\displaystyle \frac{3}{2};{S}_{x}\right. \rangle +{ \mathcal O }\left(1-\displaystyle \sum _{{S}_{x}=-\tfrac{N}{2}}^{\tfrac{N}{2}}{\left|{\rm{\Psi }}(t)\left|\displaystyle \frac{3}{2};{S}_{x}\right. \rangle | \right.}^{2}\right).\end{eqnarray} \tag{ 7 }$$

This equation stems from the fact that the quantum numbers $S=\tfrac{N}{2}$ and ${S}_{x}\in \{-\tfrac{N}{2},\,\ldots ,\,\tfrac{N}{2}\}$ uniquely identify a single N-body state of the $\vec{n}=(0,\,\ldots ,\,N-1)$ configuration. Additionally, the ${\hat{S}}_{x}$ operator commutes with both ${\hat{S}}^{2}$ and $\hat{H}$ cancelling all cross terms that would appear in the first term of equation (7). Therefore, within the ferromagnetically ordered regime where $\langle {\rm{\Psi }}(t)| {\hat{S}}^{2}| {\rm{\Psi }}(t)\rangle \approx \tfrac{N}{2}\left(\tfrac{N}{2}+1\right)$, for $| {\rm{\Psi }}(0)\rangle =| {{\rm{\Psi }}}_{F}\rangle$, all of the contributing expectation values need to satisfy $\langle \tfrac{N}{2};{S}_{x}| {{\rm{e}}}^{\tfrac{{\rm{i}}{t}}{{\hslash }}\hat{H}}{\hat{S}}^{2}{{\rm{e}}}^{-\tfrac{{\rm{i}}{t}}{{\hslash }}\hat{H}}| \tfrac{N}{2};{S}_{x}\rangle \approx \tfrac{N}{2}\left(\tfrac{N}{2}+1\right)$, since the probabilities are positive definite, i.e. $| \langle {\rm{\Psi }}(t)| \tfrac{3}{2};{S}_{x}\rangle {| }^{2}\gt 0$. This implies that for all initial states which satisfy ${\sum }_{{S}_{x}=-\tfrac{N}{2}}^{\tfrac{N}{2}}| \langle {\rm{\Psi }}(0)| \tfrac{3}{2};{S}_{x}\rangle {| }^{2}\approx 1$ the correlation dynamics within the ferromagnetically ordered regime is stable.

To provide concrete numerical evidence supporting the above-mentioned theoretical argument, we present in figure 6, the time-evolution of the polarization, P⁽¹⁾_S, and the spin–spin correlator, C⁽²⁾_S, utilizing the Hamiltonian of equation (5), in the case of the initial state [16]

$$\begin{eqnarray}\begin{array}{c}| {{\rm{\Psi }}}_{S=3/2,{Sz}=1/2} \rangle =\int \displaystyle \frac{{{\rm{d}}}^{3}x}{\sqrt{3!}}\left|\begin{array}{ccc}{\phi }_{0}({x}_{1}) & {\phi }_{0}({x}_{2}) & {\phi }_{0}({x}_{3})\\ {\phi }_{1}({x}_{1}) & {\phi }_{1}({x}_{2}) & {\phi }_{1}({x}_{3})\\ {\phi }_{2}({x}_{1}) & {\phi }_{2}({x}_{2}) & {\phi }_{2}({x}_{3})\end{array}\right|{\hat{\psi }}_{\uparrow }^{\dagger }({x}_{1}){\hat{\psi }}_{\uparrow }^{\dagger }({x}_{2}){\hat{\psi }}_{\downarrow }^{\dagger }({x}_{3})| 0 \rangle .\end{array}\end{eqnarray} \tag{ 8 }$$

As it can be clearly seen the correlation dynamics (see figure 6(b)) is almost identical to the one observed in section 3.1 for the fully spin-polarized initial state, $| {{\rm{\Psi }}}_{F}\rangle ={\prod }_{i=0}^{2}\int {\rm{d}}x\,{\phi }_{i}(x){\hat{\psi }}_{\uparrow }^{\dagger }(x)| 0\rangle$ (see also figure 1(b)). Most importantly, the ferromagnetically ordered regime appears in the $0.5\lt | g| \lt 4$ interaction range, characterized by ${C}_{S}^{2}\approx 1$, while losses of spin-alignment (i.e. C²_S < 1) are found outside of this interaction regime. The polarization dynamics (see figure 6(a)), however, shows different patterns from the dynamics obtained with $| {{\rm{\Psi }}}_{F}\rangle$, since the initial polarization in the present case is ${P}_{S}^{(2)}=\tfrac{1}{3}$ rather than unity. Nevertheless, within the ferromagnetically ordered regime we observe large fluctuations of the polarization while the spin–spin correlator is almost constant, similarly to the dynamics that the system with $| {\rm{\Psi }}(0) \rangle =| {{\rm{\Psi }}}_{F} \rangle$ follows (see figures 1(a), (b)).

Zoom In Zoom Out Reset image size

According to our previous discussion (see section 2) the Rabi-coupling between the spin-$\uparrow$ and the spin-$\downarrow$ states is assumed to be weak and the characteristic length of its modulation is larger than the length scale of the trap, $\sqrt{\tfrac{{\hslash }}{m\omega }}$. Due to these assumptions it is reasonable to approximate the Rabi-coupling potential by its Taylor series. We can, thus, demonstrate that our results generalize to all Rabi-coupling potentials with a non-vanishing second-order derivative⁷ by showing that a similar dynamics as in section 3.1 can be obtained for the parabolic spin-coupling potential

$$\begin{eqnarray}&&{\hat{H}}_{S}^{{\prime} }=\displaystyle \frac{1}{2}m{({\rm{\Delta }}\omega )}^{2}\int {\rm{d}}x\,{x}^{2}\displaystyle \sum _{\alpha ,\beta }{\hat{\psi }}_{\alpha }^{\dagger }(x){\sigma }_{\alpha \beta }^{x}{\hat{\psi }}_{\beta }(x).\end{eqnarray} \tag{ 9 }$$

For our simulations we employ ${\rm{\Delta }}\omega =\sqrt{0.1}$, while the system is initialized in the fully polarized state, $| {{\rm{\Psi }}}_{F}\rangle$, though as argued above a similar dynamics takes place when the system is initialized in the $| {{\rm{\Psi }}}_{S=3/2,{Sz}=1/2}\rangle$ state (results not shown here for brevity). As it can be seen, the behaviour of the system in terms of the spin-polarization, P⁽¹⁾_S (see figure 6(c)) and the spin–spin correlations, C⁽²⁾_S (see figure 6(d)) is almost identical to the case of ${\hat{H}}_{S}$, equation (6) (see also figures 1(a) and (b)), with deviations occurring only within the g ≈ 0 region (compare figures 6(c) and (b)). We also note that in the case of strong spin-dependent potentials, where the exact shape of the Rabi-coupling potential might play an important role, spin segregation phenomena are induced [28, 45]. These compete with the ferromagnetic order identified here, as the overlap of the spin-densities provides an upper bound for $\langle {\hat{S}}^{2}\rangle$ and therefore such investigations lie beyond the scope of this work.

Our setup can be realized using ⁴⁰K atoms under the influence of a Raman coupling of the two energetically lowest hyperfine states, while the observables P⁽¹⁾_S and C⁽²⁾_S are accessible by fluorescence imaging. Below we propose a possible experimental realization in order to probe our findings. The robustness of the suggested measurement scheme is demonstrated by comparing our MCTDHF results with simulated sets of single-shot images that contain additional noise emulating this way the noise sources inherent in a corresponding experiment [14–16].

The effective Rabi coupling scheme, see ${\hat{H}}_{S}$ equation (6), can be achieved by employing a two-photon resonant, δ = 0, Raman transition via two Gaussian focussed laser beams. To incorporate non-negligible interatomic interactions one needs to apply a bias magnetic field close to the point of an s-wave broad Fano-Feshbach resonance [8]. For ${}^{40}{\rm{K}}$ atoms a broad s-wave Fano-Feshbach resonance between the hyperfine states $| \uparrow \rangle =| {}^{2}{S}_{1/2};F=\tfrac{9}{2},{m}_{F}=-\tfrac{9}{2}\rangle$ and $| \downarrow \rangle =| {}^{2}{S}_{1/2};F=\tfrac{9}{2},{m}_{F}=-\tfrac{7}{2}\rangle$, is located at the magnetic field strength B_FF = 202.10 G [50].

Fluorescence imaging is commonly used to probe the state of the system in few-atom (N < 10) experiments [14]. Here a certain number of atoms is ejected from the trap and recaptured into a magneto-optical trap [51]. Subsequently, the number of ejected particles can be inferred by measuring the intensity of the scattered light. We show that P⁽¹⁾_S and C⁽²⁾_S can be experimentally detected using fluorescence imaging. P⁽¹⁾_S and C⁽²⁾_S depend on the average and the variance of the magnitude of the spin polarization respectively. Because of the employed Raman scheme the Hamiltonian (6) is implemented in the interaction picture. This implies that in the Schrödinger picture and in the absence of the Raman fields the orientation of the spin-vector precesses around the z spin-axis with frequency ${\omega }_{\uparrow \downarrow }=2\pi \times \left[44.8+0.156{{\rm{G}}}^{-1}(B-{B}_{\mathrm{FF}})\right]$ MHz (where B refers to the bias magnetic field). ${\omega }_{\uparrow \downarrow }$ corresponds to the energy offset between the $| \uparrow \rangle$ and $| \downarrow \rangle$ states of ${}^{40}{\rm{K}}$ for magnetic fields in the vicinity of ${B}_{\mathrm{FF}}$. As a consequence only the spin-polarization along the z axis (i.e. population-imbalance in the occupation of the hyperfine states $| \uparrow \rangle ,| \downarrow \rangle$) in spin–space can be directly probed. To measure the spin-state in such atomic systems Ramsay spectroscopy is employed to coherently rotate the rotating x or y axes in the interaction picture to the stationary z axis, which is common for both pictures. A Ramsay spectroscopy sequence (described in the interaction picture) is utilized. Initially all of the atoms are prepared in the N-body state $| {\rm{\Psi }}(0)\rangle$, namely all atoms reside in the $| \uparrow \rangle$ hyperfine state. At time t = 0 the inhomogeneous Raman coupling of the hyperfine states is suddenly switched on and the fermions are exposed to it for time t. By the end of this process, the MB wavefunction has evolved from $| {\rm{\Psi }}(0)\rangle$ to $| {\rm{\Psi }}(t)\rangle$ (in the interaction picture) under the influence of $\hat{H}$. At time t, the Raman coupling is suddenly switched off and the system evolves for a dark time t_dark. Within this time interval the reestablished symmetries of the Hamiltonian ${\hat{H}}_{0}+{\hat{H}}_{I}$ (see equations (3), (4)) prohibit any change to ${\boldsymbol{S}}$ and S². To measure the S_x or S_y components we need to rotate the desired spin component to the z axis by applying a −σ_y or ${\sigma }_{x},\tfrac{\pi }{2}$-pulse respectively by means of spatially homogeneous two-photon-resonant optical Raman fields with the appropriate phase shift, ϕ, from the inhomogeneous one (for S_z no $\tfrac{\pi }{2}$-pulse is used and the sequence continues directly with the next step). This sequence stops the precession dynamics of the desired spin component in the Schrödinger picture as it is mapped to the stationary z axis. In the following, all the spin-$\downarrow$ are removed from the trap by applying a high-intensity resonant laser pulse at time t_ex [16]. The surviving atoms are loaded into the magneto-optical trap (at t = t_meas) and counted to provide a measurement for the spin polarization ${S}_{i}=\tfrac{2{N}_{\uparrow }-N}{2}$ along the selected axis i ∈ { x, y, z }.

As a proof-of-principle of the above-mentioned imaging procedure we simulate single experimental measurements, where we take into account a random error in the phase Δ ϕ. We employ a generalized version of the recent single-shot implementation offered by Multi-Layer Multi-Configuration Hartree method for atomic Mixtures (ML-MCTDHX) [36] (see [52–55] and appendix B for details) and evaluate P_S⁽¹⁾ and C⁽²⁾_S from the simulated experimental data. Note that Δ ϕ might be induced by variations in the optical path of the −σ_y or σ_x, $\tfrac{\pi }{2}$-pulse Raman beams. To incorporate this source of error we simulate experimental measurements for each of the x, y and z components of the spin vector and incorporate a random rotation by $\delta \phi$ along the z axis that follows a Gaussian distribution of width Δ ϕ.

Figures 7(a) and (b) offer a comparison of our MCTDHF data with the simulated experimental estimates based on 600 single-shot realizations containing an error Δ ϕ = π/12. We observe that despite the latter error the single-shot results follow closely the MCTDHF data and reproduce both the spin polarization P_S⁽¹⁾, as well as, the correlation dynamics C⁽²⁾_S. The uncertainty of the simulated single-shot results is of the order ∼0.05. Agreement between the MCTDHF data and the single-shot estimates is observed also for the ferromagnetically ordered (see figures 7(c) and (d)) and the strong-g demagnetization regime (see figures 7(e) and (f)). Therefore, we conclude that the error in the phase ϕ is not prohibitive for accurate measurements of ${P}_{S}^{(1)}$ and C⁽²⁾_S as long as it is kept sufficiently smaller than $\tfrac{\pi }{2}$.

Zoom In Zoom Out Reset image size

We have explored the spin-flip dynamics of few ultracold fermions subject to spatially inhomogeneous external driving of the spins. In particular, we showed that in this case the polarization of the confined Fermi gas cannot be stabilized for any interaction strength. A result that lies in contrast to the picture of ferromagnetism provided by the celebrated Stoner model. Most importantly, a stable correlation-induced ferromagnetic spin-order emerges in spite of the strongly fluctuating polarization for moderate interactions. We have characterized the emerging spin-order by comparing ab initio simulations with an effective spin-chain model in the few-body case. The influence of correlations and the emergence of entangled NOON-like states during the dynamical evolution of the system has been explicitly demonstrated. In the weak and strong interaction limit the behaviour of the system is characterized by a significant depletion of the spin–spin correlator which can be related to the corresponding avoided crossings appearing in the eigenspectrum. Our setup is experimentally accessible in ⁴⁰K few-atom experiments by employing a Raman coupling scheme of the two energetically lowest hyperfine states. The observables P_S⁽¹⁾ and C_S⁽¹⁾ can be measured by employing Ramsay spectroscopy and fluorescence imaging.

It is known that the properties of itinerant magnetism vary depending on the confining potential [56]. Studying the stability of ferromagnetism in the case of a double well potential or an optical lattice can yield further insights into the magnetic properties exhibited in 1D systems. Notice, also, that the spin-chain model presented here is easily extendable to higher dimensional settings. The investigation on whether a similar order occurs in higher dimensional settings also provides an intriguing perspective for future study. Another interesting prospect is to examine the demagnetization dynamics of few-fermions when exposed to Rashba and Dresselhaus spin–orbit coupling [57, 58]. This might establish a link to relevant condensed matter systems where such demagnetization mechanisms are well-studied [59–61]. Such dynamics have recently been examined in the case of thermal Fermi gases in the collisionless regime [62–64].

This work is supported by the Cluster of Excellence 'Advanced Imaging of Matter' of the Deutsche Forschungsgemeinschaft (DFG) - EXC 2056 - project ID 390715994. SIM and PS gratefully acknowledge financial support by the Deutsche Forschungsgemeinschaft (DFG) in the framework of the SFB 925 'Light induced dynamics and control of correlated quantum systems'.

Our approach to solve the MB Schrödinger equation $\left({\rm{i}}{\hslash }{\partial }_{t}-\hat{H}\right)| {\rm{\Psi }}(t)\rangle =0$ relies on the ML-MCTDHX [36]. In particular we employ a reduction of the ML-MCTDHX method for spin-1/2 fermions being referred in the following as the spinor-variant of MCTDHF. MCTDHF has been applied extensively for the treatment of fermions with or without spin-degrees of freedom, in a large class of condensed matter, atomic and molecular physics scenarios (see e.g. [65–70]) and recently also applied in the field of ultracold atoms [35–39, 71]. The key idea of MCTDHF lies in the usage of a time-dependent (t-d) and variationally optimized MB basis set, which allows for the optimal truncation of the total Hilbert space. The ansatz for the MB wavefunction is taken as a linear combination of t-d Slater determinants $| \vec{n}(t)\rangle$, with t-d weight coefficients ${A}_{\vec{n}}(t)$

$$\begin{eqnarray}&&| {\rm{\Psi }}(t)\rangle =\displaystyle \sum _{\vec{n}}{A}_{\vec{n}}(t)| \vec{n}(t)\rangle .\end{eqnarray} \tag{ A1 }$$

Each t-d Slater determinant is expanded in terms of M t-d variationally optimized single-particle functions (SPFs) $\left|{\phi }_{l}(t)\right\rangle ,l\,=\,1,2,\,\ldots ,\,M$ with occupation numbers $\vec{n}=({n}_{1},\,\ldots ,\,{n}_{M})$. The SPFs are subsequently expanded in a primitive basis $\{\left|k,s\right\rangle \}$, being the tensor product of a discrete variable representation (DVR) basis for the spatial degrees of freedom $\{| k\rangle \}$ of dimension M_p and the two-dimensional spin basis $\{| \uparrow \rangle ,| \downarrow \rangle \}$

$$\begin{eqnarray}&&| {\phi }_{j}(t)\rangle =\displaystyle \sum _{k=1}^{{M}_{p}}\displaystyle \sum _{\alpha =\{\uparrow ,\downarrow \}}{C}_{k\alpha }^{j}(t)| k\rangle | \alpha \rangle .\end{eqnarray} \tag{ A2 }$$

C^j_{k α}(t) refer to the corresponding t-d expansion coefficients. Note here that each t-d SPF is a general spinor wavefunction of the form $| {\phi }_{j}(t)\rangle =\int {\rm{d}}x\,[{\chi }_{j}^{\uparrow }(x){\hat{\psi }}_{\uparrow }^{\dagger }(x)+{\chi }_{j}^{\downarrow }(x){\hat{\psi }}_{\downarrow }^{\dagger }(x)]| 0\rangle$ and hence the employed method is termed as the spinor-variant of MCTDHF. The time-evolution of the N-body wavefunction under the effect of the Hamiltonian $\hat{H}$ reduces to the determination of the A-vector coefficients and the SPFs, which in turn follow the variationally obtained MCTDHF equations of motion [34–36]. In the limiting case of M = N, the method reduces to the TDHF method, while for the case of M = 2M_p, it is equivalent to a full configuration interaction approach (commonly referred to as 'exact diagonalization' in the literature) within the basis $\{| k,s\rangle \}$.

For our implementation we have used a harmonic oscillator DVR, which results after a unitary transformation of the commonly employed basis of harmonic oscillator eigenfunctions, as a primitive basis for the spatial part of the SPFs. To study the dynamics, we propagate the wavefunction by utilizing the appropriate Hamiltonian within the MCTDHF equations of motion. To verify the accuracy of the numerical integration, we impose the following overlap criteria $\left|\langle {\rm{\Psi }}| {\rm{\Psi }}\rangle -1\right|\lt {10}^{-8}$ for the total wavefunction and $\left|\langle {\varphi }_{i}| {\varphi }_{j}\rangle -{\delta }_{{ij}}\right|\lt {10}^{-9}$ for the SPFs. To infer about convergence, we increase the number of SPFs and DVR basis states such that the observables of interest (P_S⁽¹⁾, C_s⁽²⁾) do not change within a given relative accuracy which is in our case 10⁻⁴. More specifically, we have used M_p = 60, M = 26 and M_p = 80, M = 20 for the N = 3 and the N = 5 case respectively. Note that a full configuration interaction treatment of the above-mentioned systems in the employed primitive bases would require 280 840 number-states for N = 3 and 820 384 032 ones for N = 5.

The single-shot simulation procedure relies on a sampling of the MB probability distribution, being available within the ML-MCTDHX framework. In a spinor Fermi gas the single-shot procedure is altered significantly when compared to the single component case [52–54]. Here the role of entanglement between particles in different spin states plays an important role. For example consider the procedure that the spin-$\uparrow$ atoms are imaged before the spin-$\downarrow$ atoms. Then, the total number of spin-$\uparrow$ atoms N_↑^imag that will be imaged is not a priori known due to the breaking of the S_z symmetry. However, after imaging all of the spin-$\uparrow$ atoms the number of spin-$\downarrow$ atoms is exactly known ${N}_{\downarrow }^{\mathrm{imag}}=N-{N}_{\uparrow }^{\mathrm{imag}}$ since the total number of atoms N is definite.

To capture the entanglement between the different spin states the MB wavefunction obtained by ML-MCTDHX should be expressed such that the entanglement between the spin states is evident. The spin-1/2 Fermi gas under consideration is a bipartite system [72, 73] since the spatial degree of freedom for each particle in the spin-$\uparrow$ or spin-$\downarrow$ state resides in the Fock space ${{ \mathcal F }}^{\uparrow },{{ \mathcal F }}^{\downarrow }$ respectively. The latter results in a total Fock space ${{ \mathcal F }}^{S=1/2}={{ \mathcal F }}^{\uparrow }\otimes {{ \mathcal F }}^{\downarrow }$. Then, the MB wavefunction can be expressed in the Schmidt decomposition form (herewith we omit the temporal dependence for simplicity)

$$\begin{eqnarray}&&| {\rm{\Psi }}\rangle =\displaystyle \sum _{k=1}^{K}\sqrt{{\lambda }_{k}}| {{\rm{\Psi }}}_{k}^{\uparrow }\rangle | {{\rm{\Psi }}}_{k}^{\downarrow }\rangle .\end{eqnarray} \tag{ B1 }$$

The coefficient λ_k is referred to as the natural occupation of the species function k⁸ . Note that, $| {{\rm{\Psi }}}_{k}^{\alpha }\rangle \in {{ \mathcal F }}^{\alpha }$ and as such the number of α-spin particles varies for different Schmidt modes, k. A state of the bipartite system (see equation (B1)) cannot be expressed as a direct product of two states from the two subsystem Fock spaces ${{ \mathcal F }}^{\alpha }$ if at least two coefficients ${\lambda }_{k}$ are nonzero. In the latter case the system is referred to as entangled [74]. The Schmidt decomposition of the MB wavefunction is obtained as follows. The reduced density matrix for one of the spin states, let it be α, is evaluated i.e. ${\rho }^{\alpha }={\mathrm{Tr}}_{\alpha ^{\prime} }[| {\rm{\Psi }}\rangle \langle {\rm{\Psi }}| ]$, where α' refers to the spin state orthogonal to α and subsequently diagonalized resulting in its Schmidt representation ${\rho }^{\alpha }={\sum }_{k=1}^{M}{\lambda }_{k}| {{\rm{\Psi }}}_{k}^{\alpha }\rangle \langle {{\rm{\Psi }}}_{k}^{\alpha }|$. Then, the corresponding species wavefunction of the spin state α' can be calculated by $| {{\rm{\Psi }}}_{k}^{\alpha ^{\prime} }\rangle =\tfrac{1}{\sqrt{{\lambda }_{k}}}\langle {{\rm{\Psi }}}_{k}^{\alpha }| {\rm{\Psi }}\rangle$.

The single-shot process in spinor gases represents a generalization of the single-shot process for a mixture with a definite number of atoms in each species [55]. This generalization is based on the treatment of the vacuum state $| {0}^{\alpha }\rangle$. Before each step of the single-shot process the existence of particles in the imaged spin state is checked. To perform the latter a random number in the interval P_rand ∈ [0, 1] is compared with ${\lambda }_{\tilde{k}}$, where $\tilde{k}$ is the Schmidt mode for which $| {{\rm{\Psi }}}_{\tilde{k}}^{\alpha }\rangle =| {0}^{\alpha }\rangle$ holds. If ${P}_{\mathrm{rand}}\lt {\lambda }_{\tilde{k}}$ the imaging of the spin state α ends and the MB wavefunction is projected to $| \tilde{{\rm{\Psi }}}\rangle =| {0}^{\alpha }\rangle \otimes | {{\rm{\Psi }}}_{\tilde{k}}^{\alpha ^{\prime} }\rangle$. Then the simulation of the imaging of the α' spin state is initiated. The MB wavefunction in this case is the species wavefunction $| {{\rm{\Psi }}}_{\tilde{k}}^{\alpha ^{\prime} }\rangle$ and as such the single-shot procedure reduces to the well-established single species case (see [52–54] and also the discussion below). For ${P}_{\mathrm{rand}}\gt {\lambda }_{\tilde{k}}$ a particle in the spin state α is imaged. First, a random position is drawn according to the constraint ${\rho }_{\alpha }^{(1)}({x}_{1}^{{\prime} })=\langle {\rm{\Psi }}| {\hat{\psi }}_{\alpha }^{\dagger }({x}_{1}^{{\prime} }){\hat{\psi }}_{\alpha }({x}_{1}^{{\prime} })| {\rm{\Psi }}\rangle \gt {l}_{1}$ where l₁ refers to a random number within the interval [0, max{ρ⁽¹⁾_α(x)}]. Then we project the N-body wavefunction to the (N − 1)-body one by employing the operator $\tfrac{1}{{ \mathcal N }}{\hat{\psi }}_{\alpha }({x}_{1}^{{\prime} })$, where ${ \mathcal N }=\sqrt{ \langle {\rm{\Psi }}| {\hat{\psi }}_{\alpha }^{\dagger }({x}_{1}^{{\rm{{\prime} }}}){\hat{\psi }}_{\alpha }({x}_{1}^{{\rm{{\prime} }}})| {\rm{\Psi }} \rangle }$ is a normalization factor. The latter process directly affects the Schmidt coefficients λ_k's (entanglement weights) and thus despite the fact that the spin-α' atoms have not been imaged yet, both ${\rho }_{\uparrow }^{(1)}(x)=\langle {\rm{\Psi }}| {\hat{\psi }}_{\uparrow }^{\dagger }(x){\hat{\psi }}_{\uparrow }(x)| {\rm{\Psi }}\rangle$ and ${\rho }_{\downarrow }^{(1)}(x)=\langle {\rm{\Psi }}| {\hat{\psi }}_{\downarrow }^{\dagger }(x){\hat{\psi }}_{\downarrow }(x)| {\rm{\Psi }}\rangle$ change. This can be easily understood by employing again the Schmidt decomposition. Indeed after this first measurement the (N − 1)-particle MB wavefunction reads

$$\begin{eqnarray}&&| {\tilde{{\rm{\Psi }}}}^{(-1)}\rangle =\displaystyle \sum _{k}\sqrt{{\tilde{\lambda }}_{k}^{(-1)}}| {\tilde{{\rm{\Psi }}}}_{k}^{\alpha (-1)}\rangle | {{\rm{\Psi }}}_{k}^{\alpha ^{\prime} }\rangle ,\end{eqnarray} \tag{ B2 }$$

where $| {\widetilde{{\rm{\Psi }}}}_{k}^{\alpha (-1)} \rangle =\tfrac{1}{{N}_{k}}{\hat{\psi }}_{\alpha }({x}_{1}^{{\rm{{\prime} }}})| {{\rm{\Psi }}}_{k}^{\alpha } \rangle$ refer to the species wavefunction after the imaging and ${N}_{k}=\sqrt{\langle {{\rm{\Psi }}}_{k}^{\alpha }| {\hat{\psi }}_{\alpha }^{\dagger }({x}_{1}^{{\prime} }){\hat{\psi }}_{\alpha }({x}_{1}^{{\prime} })| {{\rm{\Psi }}}_{k}^{\alpha }\rangle }$ denotes the corresponding normalization factor. Finally, the Schmidt coefficients read ${\tilde{\lambda }}_{k}^{(-1)}={\lambda }_{k}{N}_{k}/{\sum }_{m}{\lambda }_{m}{N}_{m}^{2}$. The above-mentioned procedure is repeated N^imag_α times until the condition $| | \langle {0}^{\alpha }| {{\rm{\Psi }}}^{(-{N}_{\alpha }^{\mathrm{imag}})}\rangle | | =1$ is reached or if a random number satisfying ${P}_{\mathrm{rand}}\lt {\lambda }_{\tilde{k}}^{(-{N}_{\alpha }^{\mathrm{imag}})}$ is selected. The resulting distribution of positions (${x}_{1}^{{\prime} }$,..., ${x}_{{N}_{\alpha }^{\mathrm{imag}}}^{{\prime} }$) is convoluted with a point spread function leading to a single-shot ${{ \mathcal A }}^{\alpha }(\tilde{x})$ for the spatial configuration of spin-α particles, where $\tilde{x}$ refers to the spatial coordinates within the image. It is worth mentioning at this point that, in the special case for which the probability of ${N}_{\alpha }^{\mathrm{imag}}=N$ is zero, it can be easily shown that upon annihilating the last spin-α particle (provided that ${P}_{\mathrm{rand}}\gt {\lambda }_{\tilde{k}}^{(-{N}_{\alpha }^{\mathrm{imag}}+1)}$ is chosen) the $(N-{N}_{\alpha }^{\mathrm{imag}})$-particle MB wavefunction becomes

$$\begin{eqnarray}&&| {\tilde{{\rm{\Psi }}}}^{(-{N}_{\alpha }^{\mathrm{imag}})}\rangle =| {0}^{\alpha }\rangle \otimes \displaystyle \sum _{k\ne \tilde{k}}\displaystyle \frac{\sqrt{{\tilde{\lambda }}_{k}^{(-{N}_{\alpha }^{\mathrm{imag}}+1)}}\langle {x}_{{N}_{\alpha }^{\mathrm{imag}}}^{{\prime} }| {{\rm{\Phi }}}_{k}^{\alpha }\rangle }{{\displaystyle \sum }_{m\ne \tilde{k}}\sqrt{{\tilde{\lambda }}_{k}^{(-{N}_{\alpha }^{\mathrm{imag}}+1)}| \langle {x}_{{N}_{\alpha }^{\mathrm{imag}}}^{{\prime} }| {{\rm{\Phi }}}_{k}^{\alpha }\rangle {| }^{2}}}| {{\rm{\Psi }}}_{k}^{\alpha ^{\prime} }\rangle .\end{eqnarray} \tag{ B3 }$$

After this last step the entanglement between the spin states has been destroyed and the single component wavefunction of the spin $\alpha ^{\prime}$ atoms $| {{\rm{\Psi }}}^{(-{N}_{\alpha }^{\mathrm{imag}})}\rangle$ corresponds to the second term on the right hand side of equation (B3).

In this way, it becomes evident that after the imaging of spin α particles the resulting wavefunction $| {{\rm{\Psi }}}^{(-{N}_{\alpha }^{\mathrm{imag}})}\rangle =| {{\rm{\Psi }}}^{\alpha ^{\prime} }\rangle$ (see equation (B3)) is a non-entangled (N − N_α^imag)-particle MB wavefunction and its corresponding single-shot procedure is the same as in the single species case [52]. The latter is well-established (for details see [52, 53]) and here it is only briefly outlined below. We first calculate ${\rho }_{\alpha ^{\prime} }^{(1)}(x)$ from the MB wavefunction $| {{\rm{\Psi }}}^{(-{N}_{\alpha }^{\mathrm{imag}})}\rangle =| {{\rm{\Psi }}}^{\alpha ^{\prime} }\rangle$. Then, a random position ${x}_{1}^{{\rm{{\prime} }}{\rm{{\prime} }}}$ is drawn obeying ρ⁽¹⁾_α'( ${x}_{1}^{{\rm{{\prime} }}{\rm{{\prime} }}}$ ) > l₂ where l₂ is a random number in the interval $[0,\max [{\rho }_{\alpha ^{\prime} }^{(1)}(x)]$. Next, one particle located at position ${x}_{1}^{{\rm{{\prime} }}{\rm{{\prime} }}}$ is annihilated and ρ^{(1) '}_α'(x) is calculated from $| {{\rm{\Psi }}}^{(-{N}_{\alpha }^{\mathrm{imag}}-1)}\rangle ={\left[{\rho }_{\alpha ^{\prime} }^{(1)}({x}_{1}^{{\prime\prime} })\right]}^{-1/2}{\hat{\psi }}_{\alpha ^{\prime} }({x}_{1}^{{\prime\prime} })| {{\rm{\Psi }}}^{(-{N}_{\alpha }^{\mathrm{imag}})}\rangle$. To proceed, a new random position ${x}_{2}^{{\rm{{\prime} }}{\rm{{\prime} }}}$ is drawn from ρ^{(1) '}_α'(x; t_im). Following this procedure for $N-{N}_{\alpha }^{\mathrm{imag}}$ steps we obtain the distribution of positions (${x}_{1}^{{\rm{{\prime} }}{\rm{{\prime} }}}$, ..., ${x}_{N-{N}_{\alpha }^{{\rm{imag}}}}^{{\rm{{\prime} }}{\rm{{\prime} }}}$) which is then convoluted with a point spread function resulting in a single-shot image ${{ \mathcal A }}^{\alpha ^{\prime} }(\tilde{x}^{\prime} | {{ \mathcal A }}^{\alpha }(\tilde{x}))$.

The spin-chain Hamiltonian builds upon the spin dependent eigenstates, $| {\chi }_{n}^{\alpha }\rangle$, of the non-interacting Hamiltonian ${\hat{H}}_{0}+{\hat{H}}_{S}$, where $\alpha ,n$ denote the spin and spatial modes. To simplify the notation below, we perform a rotation in spin-space by employing the unitary operator $\hat{U}={{\rm{e}}}^{{\rm{i}}\tfrac{\pi }{4}{\hat{S}}_{y}}$ such that the x axis (see main text) is mapped to the z axis and thus the spin-modes correspond to $\alpha \in \{\uparrow ,\downarrow \}$. The nth spatial mode is considered as singly occupied if either $| {\chi }_{n}^{\uparrow }\rangle$ or $| {\chi }_{n}^{\downarrow }\rangle$ is occupied, doubly occupied if both are occupied and unoccupied if neither is occupied. Then the spatial configurations are defined by $\vec{n}=({n}_{1},\,\ldots ,\,{n}_{N})$ where n_i refers to the occupied spatial modes. There are 2^N − 2D (D denotes the number of double occupations) distinct states that correspond to the same spatial configuration, $\vec{n}$, corresponding to the different available spin-configurations $\vec{\alpha }=({\alpha }_{1},\,\ldots ,\,{\alpha }_{N})$. Consequently, a basis state of the N-body system, $| \vec{n};\vec{\alpha }\rangle ={\hat{c}}_{{n}_{1}{\alpha }_{1}}^{\dagger }...{\hat{c}}_{{n}_{N}{\alpha }_{N}}^{\dagger }| 0\rangle$, is completely defined by its spin and spatial configurations $\vec{\alpha }$ and $\vec{n}$ respectively.

To derive the effective Hamiltonian of the spin-chain model, ${\hat{H}}^{\mathrm{eff}}$, we neglect all terms that couple states of different spatial configurations. The non-interacting Hamiltonian ${\hat{H}}_{0}+{\hat{H}}_{S}$ is diagonal on the basis states $| \vec{n},\vec{\alpha }\rangle$ and thus its exact form is incorporated in the effective spin-chain Hamiltonian, ${\hat{H}}^{\mathrm{eff}}$. However, the same is not true for the interaction term ${\hat{H}}_{I}$. According to the above-mentioned approximation, the general form of the effective interaction term, ${\hat{H}}_{I}^{\mathrm{eff}}$, contains all the terms in ${\hat{H}}_{I}$ that preserve the spatial configuration of the state they act on. There are only two terms in ${\hat{H}}_{I}$ that possess the latter property and are linearly independent, namely the ${\hat{H}}_{I}^{\mathrm{shift}}$ and ${\hat{H}}_{I}^{\mathrm{exc}}$ terms. ${\hat{H}}_{I}^{\mathrm{shift}}=g{\sum }_{m,n=0}^{{m}_{\max }}{u}_{{mm}}^{{nn}}{\hat{c}}_{n\uparrow }^{\dagger }{\hat{c}}_{m\downarrow }^{\dagger }{\hat{c}}_{m\downarrow }{\hat{c}}_{n\uparrow }$ accounts for the energy shift of the single-particle modes due to interaction, where ${u}_{{mp}}^{{nq}}=\int {\rm{d}}{x}\,{[{\chi }_{n}^{\uparrow }(x){\chi }_{m}^{\downarrow }(x)]}^{* }{\chi }_{p}^{\downarrow }(x){\chi }_{q}^{\uparrow }(x)$ denote the corresponding interaction integrals. ${\hat{H}}_{I}^{\mathrm{exc}}=g{\sum }_{m=0}^{{m}_{\max }}{\sum }_{n\ne m}{u}_{{nm}}^{{mn}}{\hat{c}}_{m\uparrow }^{\dagger }{\hat{c}}_{n\downarrow }^{\dagger }{\hat{c}}_{m\downarrow }{\hat{c}}_{n\uparrow }$ allows for the exchange of the single-particle modes after a collision event. Therefore, the effective Hamiltonian reads ${\hat{H}}^{\mathrm{eff}}={\hat{H}}_{0}+{\hat{H}}_{S}+{\hat{H}}_{I}^{\mathrm{eff}}$, where ${\hat{H}}_{I}^{\mathrm{eff}}={\hat{H}}_{I}^{\mathrm{shift}}+{\hat{H}}_{I}^{\mathrm{exc}}$.

To cast ${\hat{H}}^{\mathrm{eff}}$ in the spin-chain form we define the spin operators for each spatial mode ${\hat{\sigma }}_{n}^{0}={\hat{c}}_{n\uparrow }^{\dagger }{\hat{c}}_{n\uparrow }+{\hat{c}}_{n\downarrow }^{\dagger }{\hat{c}}_{n\downarrow }$, ${\hat{\sigma }}_{n}^{z}={\hat{c}}_{n\uparrow }^{\dagger }{\hat{c}}_{n\uparrow }-{\hat{c}}_{n\downarrow }^{\dagger }{\hat{c}}_{n\downarrow }$, ${\hat{\sigma }}_{n}^{+}={\hat{c}}_{n\uparrow }^{\dagger }{\hat{c}}_{n\downarrow }$ and ${\hat{\sigma }}_{n}^{-}={\hat{c}}_{n\downarrow }^{\dagger }{\hat{c}}_{n\uparrow }$. The effective MB Hamiltonian ${\hat{H}}^{\mathrm{eff}}$ conserves the spatial modes for each spatial configuration $\vec{n}=({n}_{1},\,\ldots ,\,{n}_{N})$ as it commutes with the projection operators

$$\begin{eqnarray}&&{\hat{P}}_{\vec{n}}=\displaystyle \sum _{{\alpha }_{1}=\{\uparrow ,\downarrow \}}...\displaystyle \sum _{{\alpha }_{N}=\{\uparrow ,\downarrow \}}{\hat{c}}_{{n}_{1}{\alpha }_{1}}^{\dagger }...{\hat{c}}_{{n}_{N}{\alpha }_{N}}^{\dagger }| 0\rangle \langle 0| {\hat{c}}_{{n}_{N}{\alpha }_{N}}...{\hat{c}}_{{n}_{1}{\alpha }_{1}}.\end{eqnarray} \tag{ C1 }$$

Employing this projection operator, we can derive the spin-chain Hamiltonian, ${\hat{H}}_{\vec{n}}={\hat{P}}_{\vec{n}}{\hat{H}}^{\mathrm{eff}}{\hat{P}}_{\vec{n}}$, for each configuration $\vec{n}$ with no double occupations (i.e. ${n}_{i}\ne {n}_{j},\forall \,i,j$), corresponding to the N-spin XXZ spin-chain

$$\begin{eqnarray}&&{\hat{H}}_{\vec{n}}={\varepsilon }_{\vec{n}}(g)+\displaystyle \sum _{i=1}^{N}{h}_{{n}_{i}}(g){\hat{\sigma }}_{{n}_{i}}^{z}-g\displaystyle \sum _{i=1}^{N-1}\displaystyle \sum _{j=i+1}^{N}\left[\left({J}_{{n}_{i}{n}_{j}}^{\perp }{\hat{\sigma }}_{{n}_{i}}^{+}{\hat{\sigma }}_{{n}_{j}}^{-}+{J}_{{n}_{i}{n}_{j}}^{\perp * }{\hat{\sigma }}_{{n}_{i}}^{-}{\hat{\sigma }}_{{n}_{j}}^{+}\right)+{J}_{{n}_{i}{n}_{j}}^{z}{\hat{\sigma }}_{{n}_{i}}^{z}{\hat{\sigma }}_{{n}_{j}}^{z}\right].\end{eqnarray} \tag{ C2 }$$

The spin–spin interactions are given by the overlap integrals ${J}_{{nm}}^{\perp }={u}_{{nm}}^{{mn}}$ and ${J}_{{nm}}^{z}=\tfrac{1}{4}({u}_{{mm}}^{{nn}}+{u}_{{nn}}^{{mm}})$. The interaction-dependent energy shift ${\varepsilon }_{\vec{n}}(g)$ and the local magnetic field ${h}_{{n}_{i}}(g)$ read

$$\begin{eqnarray}\begin{array}{rcl}{\varepsilon }_{\vec{n}}(g) & = & \displaystyle \sum _{i=1}^{N}\displaystyle \frac{{E}_{{n}_{i}}^{\uparrow }+{E}_{{n}_{i}}^{\downarrow }}{2}+g\displaystyle \sum _{i=1}^{N-1}\displaystyle \sum _{j=i+1}^{N}{J}_{{n}_{i}{n}_{j}}^{z},\\ {h}_{{n}_{i}}(g) & = & \displaystyle \frac{{E}_{{n}_{i}}^{\uparrow }-{E}_{{n}_{i}}^{\downarrow }}{2}+g\displaystyle \sum _{j=1,\,j\ne i}^{N}\displaystyle \frac{{u}_{{n}_{j}{n}_{j}}^{{n}_{i}{n}_{i}}-{u}_{{n}_{i}{n}_{i}}^{{n}_{j}{n}_{j}}}{4}.\end{array}\end{eqnarray} \tag{ C3 }$$

The configurations with D double occupations have to be treated separately because the creation operator of a double occupancy ${\hat{c}}_{n\uparrow }^{\dagger }{\hat{c}}_{n\downarrow }^{\dagger }$ possesses a non-trivial commutation relation with the ${\hat{\sigma }}_{n}^{0}$ one, $[{\hat{\sigma }}_{m}^{0},{\hat{c}}_{n\uparrow }^{\dagger }{\hat{c}}_{n\downarrow }^{\dagger }]=2{\delta }_{{nm}}{\hat{c}}_{n\uparrow }^{\dagger }{\hat{c}}_{n\downarrow }^{\dagger }$. In this case, it turns out that the projected Hamiltonian, ${\hat{H}}_{\vec{n}}$, is expressed in terms of a (N − 2D)-spin XXZ Hamiltonian, ${\hat{H}}_{\vec{m}}^{(2D)}$, with $\vec{m}$ being a $(N-2D)$-particle configuration composed of the singly occupied states of $\vec{n}$. The ${\hat{H}}_{\vec{n}}$ and ${\hat{H}}_{\vec{m}}^{(2D)}$ are related via the creation operator of all the double occupations ${\hat{{ \mathcal M }}}_{\vec{k}}^{\dagger }\equiv {\prod }_{j=0}^{D-1}{\hat{c}}_{{k}_{j}\downarrow }^{\dagger }{\hat{c}}_{{k}_{j}\uparrow }^{\dagger }$, where $\vec{k}$ is the vector of doubly occupied modes in $\vec{n}$, as ${\hat{H}}_{\vec{n}}={\hat{{ \mathcal M }}}_{\vec{k}}^{\dagger }{\hat{H}}_{\vec{m}}^{(2D)}{\hat{{ \mathcal M }}}_{\vec{k}}$. ${\hat{H}}_{\vec{m}}^{(2D)}$ has exactly the same form as equation (C2) but the energy shift, ${\varepsilon }_{\vec{m}}^{(2D)}(g)$, and local magnetic field, ${h}_{{m}_{i}}^{(2D)}(g)$, possess additional contributions when compared to the ones in equation (C3). Namely,

$$\begin{eqnarray}\begin{array}{rcl}{\varepsilon }_{\vec{m}}^{(2D)}(g) & = & {\varepsilon }_{\vec{m}}(g)+\displaystyle \sum _{i=1}^{D}2{\varepsilon }_{{k}_{i}}^{0}+g\left[\displaystyle \sum _{i=1}^{D}2{J}_{{k}_{i}{k}_{i}}^{z}+\displaystyle \sum _{i=1}^{D-1}\displaystyle \sum _{j=i+1}^{D}4{J}_{{k}_{i}{k}_{j}}^{z}+\displaystyle \sum _{i=1}^{D}\displaystyle \sum _{j=1}^{N-2D}{J}_{{k}_{i}{m}_{j}}^{z}\right],\\ {h}_{{m}_{i}}^{(2D)}(g) & = & {h}_{{m}_{i}}(g)+g\displaystyle \sum _{j=1}^{D}\displaystyle \frac{{u}_{{k}_{j}{k}_{j}}^{{m}_{i}{m}_{i}}-{u}_{{m}_{i}{m}_{i}}^{{k}_{j}{k}_{j}}}{2}.\end{array}\end{eqnarray} \tag{ C4 }$$

The weight of each spatial configuration to the MB wavefunction $| {\rm{\Psi }}(t)\rangle$ is constant in time as the ${w}_{\vec{n}}=\langle {\rm{\Psi }}(t)| {\hat{P}}_{\vec{n}}| {\rm{\Psi }}(t)\rangle$ are conserved. Therefore, the time evolution of the MB wavefunction within the spin-chain approximation reads

$$\begin{eqnarray}&&| {\rm{\Psi }}(t)\rangle =\displaystyle \sum _{\vec{n}}\sqrt{{w}_{\vec{n}}}\,{e}^{-i{\hat{H}}_{\vec{n}}t}| {{\rm{\Psi }}}_{\vec{n}}(0)\rangle ,\end{eqnarray} \tag{ C5 }$$

where $| {{\rm{\Psi }}}_{\vec{n}}(0)\rangle =\tfrac{1}{\sqrt{{w}_{\vec{n}}}}{\hat{P}}_{\vec{n}}| {\rm{\Psi }}(0)\rangle$ is the normalized initial wavefunction for each XXZ spin-chain.

The generalization of the presented method compared to the one developed in [28] is the inclusion of the interaction-dependent local magnetic potential (see equations (C2) and (C4)), which vanishes for a linear gradient as the one considered in [28] in the present case such a term is important for obtaining the correct behaviour of the polarization magnitude ${P}_{S}^{(1)}$ in the ferromagnetically ordered regime. Within our implementation we numerically diagonalize the one-body Hamiltonian, ${\hat{H}}_{0}+{\hat{H}}_{S}$, by employing the basis consisting of the 80 energetically lowest eigenstates of the harmonic oscillator and truncate the summation over $\vec{n}$ of equation (C5) by taking into account only the contributions of the spatial mode configurations with $| {w}_{\vec{n}}| \gt {10}^{-12}$. This truncation results in 1520 and 38304 configurations for N = 3 and N = 5 fermions respectively.