Probing Ferromagnetic Order in Few-Fermion Correlated Spin-Flip Dynamics

We unravel the dynamical stability of a fully polarized one-dimensional ultracold few-fermion spin-1/2 gas subjected to inhomogeneous driving of the itinerant spins. Despite the unstable character of the total spin-polarization the existence of an interaction regime is demonstrated where the spin-correlations lead to almost maximally aligned spins throughout the dynamics. The resulting ferromagnetic order emerges from the build up of superpositions of states of maximal total spin. They comprise a decaying spin-polarization and a dynamical evolution towards an almost completely unpolarized NOON-like state.

Introduction.-Magnetism constitutes a principal feature of a large class of materials and represents a macroscopic phenomenon of quantum origin [1][2][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][6][7]. However, the nature of the 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, ferromagnetic properties have been observed indirectly either by the spectroscopic study of repulsive Fermi polarons [9,10] (supplemented by [11]) or by employing a Fermi gas prepared in a magnetic domain wall structure [12]. 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 quasione dimensional gas [13]. Owing to its one-dimensional (1D) character, a two-body bound state for effectively repulsive interactions is absent and thus the molecule formation is suppressed. Moreover, the experimental [14,15] and theoretical [16][17][18][19][20] study of the magnetic properties of few-fermion systems has led to the insight that for near zero and infinite interactions there is an oneto-one correspondence between the spatial representation of the system and its description in terms of a spin-chain model [21][22][23][24]. Most importantly, these spinor systems possess experimentally accessible eigenstates of ferromagnetic nature [25], 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 for small spin-flip coupling and g approaches a XXZ model [24], which allows for spin-flip transitions between all states with maximal total spin S 2 .
systems. The latter in turn is important to advance also for higher dimensional settings.
In this Letter 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 (see below). An argumentation based on the Hartree-Fock (HF) framework of Stoner's instability testified within the time-dependent Hartree-Fock (TDHF) [26] showed that the spin-polarization of the Fermi gas is stable for interparticle repulsions [27] that exceed the kinetic and spin-flip contributions [see also Fig. 1(a)]. We demonstrate that within the latter interaction regime where TDHF predicts stable ferromagnetism the many-body (MB) spin-spin correlator exhibits ferromagnetic spin-correlations throughout the dynamics. Moreover and in contrast to the TDHF re-sults, the MB state of the ferromagnetically correlated gas exhibits an unstable polarization fluctuating between fully polarized and almost completely unpolarized. We show that the decay of the polarization can be understood by generalizing the spin-chain model of Ref. [24]. Virtual transitions to lower spin-S values are 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 [see also Fig. 1(b)] i.e. a NOON-like state. Our setup can be implemented in state-of-the-art 40 K experiments and the corresponding findings can be probed by fluorescence imaging techniques [31].
The contact interaction of effective interaction strength g between the distinct spin components is contained in .Ĥ exhibits the following symmetries. It commutes with each of the components ofŜ = 2 αα dxψ † α (x)σ α,α ψ α (x), with σ denoting the Pauli vector. Second, due to the SU(2) symmetryŜ 2 Therefore, the spin-polarized state, refers to the i-th eigenfunction of the 1D harmonic oscillator, is an eigenstate ofĤ and the ground state with total spin eigenvalues S = S z = N 2 . To probe the stability of the initial ferromagnetic state |Ψ(0) = |Ψ F under coherent processes that break botĥ S z andŜ 2 symmetries we employ an inhomogeneous Rabi coupling between the spin ↑ and ↓ states. The employed Hamiltonian readsĤ =Ĥ +Ĥ S wherê H S induces spin-flip transitions with a spatially dependent transition amplitude, modelled by a Gaussian of width w and intensity B 0 . We choose the values w = 2 and B 0 = 2.5/ √ 8π (in harmonic oscillator units) leading to a precession (Larmor) frequency, ω L ≈ 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 [24] occurring when the length scale of the modulation w is smaller than that of the trap l ω = /mω. To inspect the stability of the ferromagnetism we track two main observables, directly related with the system's broken symmetries. Namely, the normalized spin polarization magnitude P (1) S = 2 N | Ŝ | and the spin-spin correlator C S expresses the averaged spin-order (magnetization) and refers to the magnitude of the polarization. The correlator C (2) S probes the alignment of each two spins and serves as an indicator for the distinction of ferromagnetic C  S appear over regular time intervals for later times. Our results for this interaction interval are compliant with the demagnetization dynamics analyzed in [24,35] and we shall refer to this regime as the weak-g demagnetization regime. For intermediate interactions of either sign, i.e. 0.5 < |g| < 4, the decay and revival of P (1) S 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 (2) S indicates that the spins are close to be maximally aligned (e.g. C (2) S ≥ 0.85 for 2 < g < 4 and C (2) S ≥ 0.95 for g ∼ −2) throughout the evolution, signifying ferromagnetic spin-correlations. Therefore, ferromagnetism is unstable in this interaction interval as the polarization (P 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-correlations rather than the polarization. For g > 4 a suppression of the ferromagnetic spin-correlations occurs as the amplitude of the C (2) S oscillations increases for stronger g, see Fig. 2(b). For instance, at g ≈ 10, C  is also oscillatory taking values between unity and 1/3, with a significantly smaller oscillation frequency than C (2) S [see Fig. 2(a)]. In the following this interaction interval (g > 4) is referred to as the strong-g demagnetization regime.
At this point it is instructive to compare the MB results with the corresponding TDHF approximation presented in Fig. 2(c) and 2(d). For weak |g| (|g| < 0.5) the demagnetization dynamics is qualitatively captured by the TDHF approximation. However, upon increasing |g|, TDHF predicts no loss of P  Concluding, the ferromagnetically ordered regime exhibited in the MB case corresponds to a stable ferromagnetic one within the HF framework. In both cases the interaction regime is limited to intermediate values of g, 0.5 < |g| < 4 [see Fig. 2(b), 2(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 [36,37]. Instead, it is known that the states of maximum S [see Fig. 1(b)] and the MB antiferromagnetic ground states exhibit an avoided crossing in the Tonks-Girardeau limit [15][16][17][18][19][20] due to the breaking of the SU(2) symmetry. As we shall argue later on, the fluctuations of C (2) S in the strong-g demagnetization regime for the MB case can be attributed to this avoided crossing.
To uncover the main mechanisms responsible for the emergence of the ferromagnetically ordered regime (0.5 < |g| < 4) we employ an extended version of the spinchain model [31] presented in Ref. [24]. Our spinchain model incorporates additional effective magnetic field terms (for details see [31]) when compared to [24], that are necessary for the treatment of generic spin-1/2 fermion systems with a single conserved spin-component (hereŜ x ). Within this model the N -body wavefuntion is decomposed as of H 0 +Ĥ S . n = (n 1 , . . . , n N ) and α = (α 1 , . . . , α N ) denote the spatial and spin configuration respectively. The crucial approximation in this model is that all the interaction terms,Ĥ δ n I , that couple different spatial configurations n can be neglected. By settingĤ eff I ≡Ĥ I −Ĥ δ n I , the different spatial configurations n are decoupled. It can be shown that the time-evolution of each |Ψ n (t) can be described by an XXZ spin-chain consisting of N -spins [24]. The employed approximation limits the expected range of validity of the spin-chain model to small interaction values, g < ω 3 /m = 1, where the interaction energy is smaller than the energy spacing between the singleparticle eigenstates. Note here that for the particular values of B 0 , w, the dominant configuration n = (0, 1, 2) possesses approximately 99.72% of the contribution to |Ψ(0) .
The spin-dynamics within the spin-chain model consists of a weak-g demagnetization (|g| < 0.5) and a ferromagnetically ordered regime (0.5 < |g| < 4), complying with the MB case [compare [see Fig. 2(b) and 2(f)] between the two methods we observe that the spincorrelation dynamics is almost identical in the weak-g demagnetization regime. On the contrary, in the ferromagnetically ordered regime the ferromagnetic spin-correlations are overestimated by the spin-chain method [hardly visible in Fig. 2(b) and 2(f)]. For increasing interactions (g > 4) no strong-g demagnetization regime appears within the spin-chain model, because the coupling with the antiferromagnetic ground states is neglected as they correspond to a different spatial configuration. In conclusion, the MB polarization P (1) S dynamics within the ferromagnetically ordered regime is well-captured by the spin-chain model and emerges due to the couplings among the spin-states within the dominant n = (0, 1, 2) spatial configuration. In contrast to this, the (small) depletion of C (2) S stems from the neglected couplings to different spatial configurations.
To identify the underlying mechanisms of the MB spindynamics we invoke the spectrum of the spin-polarization F{P Fig. 3(a) and 3(b) for −2 < g < 10 and |g| < 1 respectively. 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 n = (0, 1, 2). For g = 0 a multitude of interactiondependent frequency branches emerge from each Larmor frequency. The three branches in the vicinity of Ω ≈ 0.85 [see also Fig. 3(b) for g ≈ 1] represent the dominant contributions within the ferromagnetically ordered regime. The remaining frequencies appearing in this regime [ Fig.  3(b)] are negligible as they correspond to eigenenergy differences between states of different S, which become energetically well-separated due to the interaction.
To further expose the modes participating in the P (1) S dynamics, we compare the MB polarization spectrum F{P (1) S }, with the eigenstates of the spin-chain model. The latter eigenstates can be labeled as |S, S x in the repulsive ( = +) and attractive ( = −) part of the ferromagnetically ordered regime. S x is conserved due to [Ŝ x ,Ĥ] = 0 and S is an approximate good quantum number due to the aforementioned energy separation. Figure  3(c) reveals that the energy differences between the spinchain eigenstates of S ≈ 3 2 , ∆E{| 3 2 , S x , | 3 2 , S x }, possess the main contribution to F{P (1) S } in the ferromagnetically ordered regime as they match well the dominant branches of the MB dynamics. The ferromagnetic spin-correlations emanate from the occupation of states with S ≈ 3/2 each characterized by a different value of S x . The frequency difference ∆Ω between the highest and lowest lying of the above-mentioned dominant branches [see also Fig. 3(a)] results in the decay dynamics of P ∆Ω is attributed to the energy shift of the | 3 2 , ± 1 2 states, induced by virtual transitions (due to the large energetic separation) to states of S ≈ 1 2 . We observe that for increasing |g| within the ferromagnetically ordered regime, ∆Ω decreases (or equivalently τ D increases) giving rise to the observed P   Fig. 3, enables us to further characterize the superpositions that emerge in the ferromagnetically ordered regime. To remove the effect of the collective Larmor precession, we introduce the precessing basis . Note that when the particles are collectively precessing with the frequency Ω L then | 3 2 , S z (t) = 3 2 |Ψ(t) | is constant in time. However, Fig. 4(a) presenting our MCTDHF results for a representative case (g = 2) within the ferromagnetically ordered regime testifies otherwise. We observe a low-frequency (∼ ∆Ω) population transfer from the state | 3 2 , S z (t) = 3 2 to the state | 3 2 , S z (t) = − 1 2 . For t > 350 the latter mechanism results in Ŝ z (t) ≈ 0, as | 3 2 , S z (t) = − 1 2 possesses approximatively a three times larger population as compared to | 3 2 , S z (t) = 3 2 . The nature of this superposition can be understood by transforming to the orthogonal precessing axis, y (t) [see Fig. 4(b)]. In this case the populations of | 3 2 , S y (t) = 3 2 and | 3 2 , S y (t) = − 3 2 are almost equal for t > 350, signifying the tendency to dynamically approach a NOON state, characterized by Ŝ y (t) ≈ 0, ) with a relative phase φ. These results combined with the conserved quantity S x = 0 explain the decay of the total magnetization P (1) S → 0 for t > 350. Resultingly, the spin-dynamics within the ferromagnetically ordered regime describes the dynamical evolution of a fullypolarized state to a superposition one consisting of two antiparallel-oriented fully-polarized states.
In the strong-g demagnetization regime [ Fig. 3(a)] the frequency branches that correspond to the | 3 2 , ± 1 2 + states deviate from Ω ≈ 0.85 as they couple to the antiferromagnetic ground states of S ≈ 1 2 and S x = ± 1 2 character. Due to the same mechanism additional branches also appear in F{P  Fig.  2(b)]. Finally let us note that all of the aforementioned results can be verified for larger particle numbers e.g. N = 5 [31].
Summary and Outlook.-By exploring the spin-flip dynamics of a few ultracold fermions we showed that the polarization of the confined Fermi gas cannot be stabilized by strong interactions. A stable correlation-induced ferromagnetic spin-order enstablishes itself in spite of the strongly fluctuating polarization. Our setup is experimentally accessible in 40 K few-atom experiments by employing a Raman coupling scheme of the two energetically lowest hyperfine states [31]. The observables P (1) S and C (1) S can be measured by employing Ramsay spectroscopy and fluorescence imaging [31]. A future interesting prospect is to examine the demagnetization dynamics of few-fermions when exposed to Rashba and Dresselhaus spin-orbit coupling [38,39] thereby providing a link to relevant condensed matter systems [40][41][42], also recently examined for thermal Fermi gases in the collisionless regime [43][44][45].
Inspecting F{P S } for N = 5 fermions, see Fig. S1(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 | 5 2 , S x states) possess a significant amplitude for |g| > 0.5. The frequency difference, ∆Ω, between the highest and lowest lying of the above five branches [see Fig. S1(g)] gives rise to the decay of P (1) S within the ferromagnetically ordered regime observed in Fig. S1(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 | 5 2 , S j (t) = e i 1 2Ŝ xΩLt | 5 2 , S j , Ω L = 1 5 ∆E{| 5 2 , 5 2 , | 5 2 , − 5 2 } and expand the MB wavefunction in terms of the latter. Figures S1(h) and S1(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.

II. Experimental Realization
Our setup can be realized using 40 K atoms under the influence of a Raman coupling of the two energetically lowest hyperfine states, while the observables P (1) S and C (2) S are accessible by fluorescence imaging. The effective Rabi coupling scheme, seeĤ S in Eq. (1), can be achieved by employing a resonant, δ = 0 twophoton 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 K atoms a broad s-wave Fano-Feshbach resonance between the hyperfine states | ↑ = | 2 S 1/2 ; F = Fluorescence imaging is commonly used to probe the state of the system in few-atom (N < 10) experiments [13]. Here a certain number of atoms is ejected from the trap and recaptured into a magneto-optical trap [S2]. Subsequently, the number of ejected particles can be inferred by measuring the intensity of the scattered light. We show that P S depend on the average and the variance of the magnitude of the spin polarization respectively. Because of the employed Raman scheme the Hamiltonian,Ĥ =Ĥ 0 +Ĥ S +Ĥ I (see main text) 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 spinvector precesses around the z spin-axis with frequency ω ↑↓ = 2π × 44.8 + 0.156G −1 (B − B FF ) MHz (where B refers to the bias magnetic field). ω ↑↓ corresponds to the energy offset between the | ↑ and | ↓ states of 40 K for magnetic fields in the vicinity of B FF . As a consequence only the spin-polarization along the z axis (i.e. population-imbalance in the occupation of the hyperfine states | ↑ , | ↓ ) 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 |Ψ(0) , namely all atoms reside in the | ↑ 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 |Ψ(0) to |Ψ(t) (in the interaction picture) under the influence ofĤ. At time t, the Raman coupling is suddenly switched off and the system is evolved for a dark time t dark . Within this time interval the reestablished symmetries of the Hamiltonian H 0 +Ĥ I (see also main text) prohibit any change to S and S 2 (this is equivalent to the statement δ = 0). 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 σ x , π 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 t σ 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-↓ are removed from the trap by applying a high-intensity resonant laser pulse at time t ex [15]. 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 = 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) [34], see [S3-S6] and section IV for details and evaluate P (1) S and C (2) S from the simulated experimental data. Note that ∆φ might be induced by variations in the optical path of the −σ y or σ x , π 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 δφ along the z axis that follows a Gaussian distribution of width ∆φ. Figures S2(a) and S2(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 (1) S , as well as, the correlation dynamics C (2) S . The uncertainty of the simulated single-shot results is of the order ∼ 0.05. Agreement between the MCTDHF data and the singleshot estimates is observed also for the ferromagnetically ordered [see Fig. S2(c) and S2(d)] and the strong-g demagnetization regime [see Fig. S2(e) and S2(f)]. Therefore, we conclude that the error in the phase φ is not prohibitive for accurate measurements of P 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. [S7-S12]) and recently also benchmarked in the field of ultracold atoms [33,34,S13].
The key idea of MCTDHF lies in the usage of a timedependent (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 | n(t) , with t-d weight coefficients A n (t) 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 |φ 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 Hamilto-nianĤ reduces to the determination of the A-vector coefficients and the SPFs, which in turn follow the variationally obtained MCTDHF equations of motion [32,33,34]. In the limiting case of M = N , the method reduces to the t-d Hartree-Fock 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 }.
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 | Ψ|Ψ − 1| < 10 −8 for the total wavefunction and | ϕ i |ϕ j − δ ij | < 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 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 [S3-S5]. Here the role of entanglement between particles in different spin states plays an important role. For example consider the procedure that the spin-↑ atoms are imaged before the spin-↓ atoms. Then, the total number of spin-↑ 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-↑ atoms the number of spin-↓ atoms is exactly known N imag ↓ = N − N 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 [S14, S15] since the spatial degree of freedom for each particle in the spin-↑ or spin-↓ state resides in the Fock space F ↑ , F ↓ respectively. The latter results in a total Fock space F S=1/2 = F ↑ ⊗ F ↓ . Then, the MB wavefunction can be expressed in the Schmidt decomposition form (herewith we omit the temporal dependence for simplicity) The coefficient λ k is referred to as the natural occupation of the species function k [S16]. Note that, |Ψ α k ∈ F α and as such the number of α-spin particles varies for different Schmidt modes, k. A state of the bipartite system [see Eq. (S3)] cannot be expressed as a direct product of two states from the two subsystem Fock spaces F α if at least two coefficients λ k are nonzero. In the latter case the system is referred to as entangled [S17]. 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. ρ α = Tr α [|Ψ Ψ|], where α refers to the spin state orthogonal to α and subsequently diagonalized resulting in its Schmidt rep- Then, the corresponding species wavefunction of the spin state α can be calculated by |Ψ α 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 [S6]. This generalization is based on the treatment of the vacuum state |0 α . 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 λk, wherẽ k is the Schmidt mode for which |Ψ α k = |0 α holds. If P rand < λk the imaging of the spin state α ends and the MB wavefunction is projected to |Ψ = |0 α ⊗ |Ψ α k . Then the simulation of the imaging of the α spin state is initiated. The MB wavefunction in this case is the species wavefunction |Ψ α k and as such the single-shot procedure reduces to the well-established single species case (see [S3-S5] and also the discussion below). For P rand > λk a particle in the spin state α is imaged. First, a random position is drawn according to the constraint ρ (1) α (x 1 ) = Ψ|ψ † α (x 1 )ψ α (x 1 )|Ψ > l 1 where l 1 refers to a random number within the interval [0, max{ρ (1) α (x)}]. Then we project the N -body wavefunction to the (N − 1)-body one by employing the operator 1 Nψ α (x 1 ), where N = Ψ|ψ † α (x 1 )ψ α (x 1 ) |Ψ 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 ρ ↓ (x) = Ψ|ψ ↓ (x)ψ ↓ (x)|Ψ change. This can be easily understood by employing again the Schmidt decomposition. Indeed after this first measurement the (N − 1)particle MB wavefunction reads where |Ψ Niψ α (x 1 ) |Ψ α k refer to the species wavefunction after the imaging and N k = denotes the corresponding normalization factor. Finally, the Schmidt coefficients readλ The above-mentioned procedure is repeated N imag ) is convoluted with a point spread function leading to a single-shot A α (x) for the spatial configuration of spin-α particles, wherex 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 imag α = N is zero, it can be easily shown that upon annihilating the last spin-α particle (provided that P rand > λ (−N imag α +1) k is chosen) the (N − N imag α )-particle MB wavefunction becomes After this last step the entanglement between the spin states has been destroyed and the single component wavefunction of the spin α atoms |Ψ (−N imag α ) corresponds to the second term on the right hand side of Eq. (S5).
In this way, it becomes evident that after the imaging of spin α particles the resulting wavefunction |Ψ (−N imag α ) = |Ψ α [see Eq. (S6) and (S5)] is a nonentangled (N − N imag α )-particle MB wavefunction and its corresponding single-shot procedure is the same as in the single species case [S3]. The latter is well-established (for details see [S3, S4]) and here it is only briefly outlined below. We first calculate ρ  α (x 1 ) > l 2 where l 2 is a random number in the interval [0, max[ρ (1) α (x)]]. Next, one particle located at position x 1 is annihilated and ρ The weight of each spatial configuration to the MB wavefunction |Ψ(t) is constant in time as the w n = Ψ(t)|P n |Ψ(t) are conserved. Therefore, the time evolution of the MB wavefunction within the spin-chain approximation reads w n e −iĤ n t |Ψ n (0) , where |Ψ n (0) = 1 √ w nP n |Ψ(0) is the normalized initial wavefunction for each XXZ spin-chain.
The generalization of the presented method compared to the one developed in Ref. [24] is the inclusion of the interaction-dependent local magnetic potential [see Eq. (S7) and (S9)], which vanish for a linear gradient as the one considered in Ref. [24]. However, in the present case such a term is important for obtaining the correct behaviour of the polarization magnitude P (1) S in the ferromagnetically ordered regime. Within our implementation we numerically diagonalize the one-body Hamiltonian,Ĥ 0 +Ĥ S , by employing the basis consisting of the 80 energetically lowest eigenstates of the harmonic oscillator and truncate the summation over n of Eq. (S10) by taking into account only the contributions of spatial mode configurations with |w n | > 10 −12 . This truncation results in 1520 and 38304 configurations for N = 3 and N = 5 fermions respectively.