Interplay of phase separation and itinerant magnetism for correlated few fermions in a double-well

We explore the stability of the phase separation phenomenon in few-fermion spin-$1/2$ systems confined in a double-well potential. It is shown that within the SU(2) symmetric case, where the total spin is conserved, the phase separation cannot be fully stabilized. An interaction regime characterized by metastable phase separation emerges for intermediate interactions which is inherently related with ferromagnetic spin-spin correlations emanating within each of the wells. The breaking of the SU(2) symmetry crucially affects the stability properties of the system as the phase separated state can be stabilized even for weak magnetic potential gradients. Our results imply an intricate relation between the phenomena of phase separation and ferromagnetism that lies beyond the view of the Stoner instability.


I. INTRODUCTION
Understanding the properties of itinerant magnetism has been a long-standing problem in condensed matter physics [1,2]. Its importance stretches beyond this field of study since it impacts the behaviour of a large class of quantum systems encountered e.g. in atomic physics [3,4]. The emergence of ferromagnetism in systems of spatially delocalized short-range repulsively interacting spinor fermions has been historically qualitatively understood in the framework of the Stoner instability [5]. Within this framework ferromagnetism is related to the phase separation of the different spin components and the formation of ferromagnetic domains [6][7][8]. Ultracold atoms provide a fertile platform to investigate such quantum many-body (MB) phenomena due to their exceptional tunability [3]. Indeed, several experiments utilizing ensembles of ultracold fermions have attempted to implement and study the Stoner instability [9][10][11][12][13] but their results have been somewhat inconclusive [14,15].
The phase separation of Fermi systems has been studied in the case of strong attractive interactions [16][17][18] where the phenomenon of spin-segregation for weak attraction or repulsion has been identified [19,20]. However, only recently experiments attempted to address the relation between ferromagnetism and phase separation in the case of a repulsively interacting Fermi-gas [11,13]. For instance, it has been demonstrated [11] that an artificially prepared phase separated state becomes metastable for strong repulsions which in turn implies the presence of a ferromagnetic instability. Accordingly, by employing pump-probe spectroscopy the emergence of short-range two-body anti-correlations in the repulsive Fermi-gas supporting some sort of ferromagnetic order has been revealed [13], while the possibility of macroscopic phase separation has been ruled out. These experimental evidences indicate that the relation between phase separation and magnetism might be more intricate and involved than it appears within the framework of the Stoner instability manifested within the Hartree-Fock theory. Nevertheless, competing processes such as the Feshbach molecule formation [21] and its possible enhancement by coherent processes [14] have hindered the experimental progress in this direction. As a consequence a complete understanding on how and via which mechanism phase separation and ferromagnetism are related remains still elusive.
Here we propose that one-dimensional (1D) few-body systems offer an ideal platform to provide insight into these fundamental questions. Besides the suppression of the above-mentioned competing processes which render the magnetic properties of 1D spin-1/2 fermions experimentally addressable [22,23], the corresponding theoretical understanding of these properties is also advanced. Indeed, the availability of numerically-exact methods [24][25][26] and the development of powerful spin-chain models [27][28][29][30][31][32][33][34] allows for the accurate modeling of the magnetic properties emerging in 1D systems in the cases of strong [27][28][29][30][31][32] and weak [33,34] interactions. Regarding the occurrence of phase separation previous studies revealed the role of the breaking of the SU(2) symmetry, associated with the conservation of the total spin of the system. Moreover, manifestations of the interplay between the magnetic properties and the phase separation have also been reported [33][34][35][36][37][38]. Below, we provide some characteristic examples. It has been demonstrated [35,36] that contrary to mean-field treatments phase separation does not occur during the interaction-quench dynamics of an SU(2) symmetric system. However, the ground state of a system with weakly broken SU(2) symmetry is known [37,38] to be phase separated in the case of infinite repulsion. In contrast, it has been shown that a parabolically confined initially spin-polarized Fermi-gas in the case of weak interactions prefers a state of largely miscible spin components even when perturbed by a spin dependent potential which weakly breaks the SU(2) symmetry [33]. In particular, for sufficiently weak spin-dependent potentials a ferromagnetic order despite the miscible character of the Fermi-gas has been established [34]. However, a systematic study that clarifies the relation between the phase separation and the magnetic properties of 1D fermions unifying, also, the above results is currently absent. Furthermore, the comparison of the underlying mechanisms provided by such a unification with the expectations of the Stoner instability might provide invaluable insights into the study of magnetic phenomena emanating in more complex systems.
Here we attempt to bridge this apparent gap in the literature by studying the stability of the phase separated state during the correlated dynamics of fermionic ensembles confined in a double-well (DW). The employed DW confinement allows for the experimental implementation of the phase separated initial state [11]. This initial state is allowed to evolve for different values of the interaction strength and the degree of the dynamical phase separation between the spin components is monitored. To capture the correlated out-of-equilibrium dynamics of this spinor fermion system we resort to the multilayer multiconfiguration time-dependent Hartree method for atomic mixtures (ML-MCTDHX) [26]. Focussing on an SU(2) invariant system and following the above-mentioned procedure we find that for weak interactions the phase separation is unstable. While for increasing repulsion an interaction regime where the phase separated state becomes metastable is unveiled. To identify the emergence of this metastable state and its relation with the magnetic properties of the system we invoke an effective tight-binding model. The metastability of the phase separated state is shown to be inherently connected with the appearance of a quasi-degenerate manifold of eigenstates characterized by intra-well ferromagnetic correlations of both wells but a varying total spin. The occurrence of this manifold is attributed to the ferromagnetic Hund exchange interactions [39][40][41] emanating within each well of the DW setup. Moreover, the low-frequency tunneling dynamics that leads to the decay of the metastable initial state provides a manifestation of the antiferromagnetic Anderson kinetic exchange interactions [42]. These interactions act between the wells and result in the lifting of the degeneracy among states exhibiting intra-well ferromagnetic correlations.
For larger interactions, the interband coupling introduced by cradle-like processes [43][44][45] is shown to result in a fastly decaying dynamics of the phase separation, thus limiting the interaction regime where this metastability of the initial state is exhibited. The breaking of the SU(2) symmetry is found to substantially affect the dynamics of the system. Indeed, the initial phase separated state of the system can be stabilized by applying a linear magnetic potential gradient to the system. This stabilization is much more prevalent in the case of intermediate interactions due to the occurrence of quasi-degenerate eigenstates with different total spin. Our results demonstrate the relation of the phase separation to the stability of the intra-well ferromagnetic order. Indeed, the interplay of the Anderson and Hund exchange interactions is found to dictate the behaviour of the system in terms of these two above phenomena implying that their relation is more intricate than what is qualitatively expected in view of the Stoner instability.
This paper is structured as follows. In section II we introduce our setup and discuss its inherent spin symmetries. Section III presents the MB dynamics of our system and showcases the important features of the related eigenspectrum. An effective tight-binding model of our system is introduced in section IV which is subsequently utilized to expose the magnetic properties of the system during the dynamics. In section V we study the dynamics in the case of a broken SU(2) symmetry. Finally, in section VI we conclude and provide future perspectives. In Appendices A and B we generalize our results for more particles and different barrier heights respectively. Appendix C provides the derivation of the Anderson effective kinetic exchange interaction for our DW setup and Appendix D describes the employed numerical approach, namely the ML-MCTDHX method.

A. Hamiltonian and Symmetries
We consider an interacting system consisting of N spin-1/2 fermions of mass m being confined in an 1D DW trap. The latter is composed by a harmonic oscillator with frequency ω and a Gaussian barrier. Such a system is described by the MB HamiltonianĤ =Ĥ 0 +Ĥ I , whereĤ 0 andĤ I correspond to its non-interacting and interacting parts respectively. The Hamiltonian,Ĥ, expressed in harmonic oscillator units ( = m = ω = 1), readŝ whereψ α (x) denotes the fermionic field operator with spin-α ∈ {↑, ↓}. V 0 and w refer to the height and width of the Gaussian barrier respectively. In the ultracold regime, g describes the effective 1D s-wave contact interaction strength between anti-aligned spins. This effective interaction strength, g, is known to be related with the transverse confinement length and the 3D s-wave scattering length [46]. The above imply that the interaction strength is experimentally tunable with the aid of confinement-induced and Fano-Feshbach resonances [21]. The Hamiltonian of Eq. (1) is invariant under rotations in spin-space as it commutes with the totalŜ z ,Ŝ ± =Ŝ x ± iŜ y spin operators. The corresponding individual spin operators,Ŝ k , are defined aŝ with σ k , k ∈ {x, y, z}, referring to the corresponding Pauli matrix. The system additionally possesses an SU(2) symmetry sinceĤ [Eq. (1)] commutes with the total spin operator,Ŝ 2 =Ŝ +Ŝ− +Ŝ z (Ŝ z − 1). As we shall demonstrate later on, this symmetry has a deep impact on the eigenspectrum of the system. The behaviour of the single-particle HamiltonianĤ 0 for varying V 0 and w is well-known [47,48] and depicted in Fig. 1 (a). For V 0 = 0 the harmonic oscillator potential is retrieved and the single-particle spectrum consists of equidistant states. As V 0 is increased, gradually all the eigenenergies, starting with the energetically two lowest ones, form quasi-degenerate pairs of different parity states (herewith called bands). Employing linear combinations of the two eigenstates forming the band, b, it is possible to construct the so-called Wannier states, φ b s (x), which are localized either in the left, s = L or the right well, s = R [49]. The single-particle dynamics of a system initialized in such a Wannier state is rather simple as the particle tunnels from each well to the other during the evolution with a frequency given by the energy difference, 2t b , between the two quasi-degenerate states which form the corresponding band.

B. Initial State Characterization
The purpose of this work is to examine whether a phase separated state can be stabilized in the presence of interactions and reveal its relation to the (ferro)magnetic properties of the system. A promising candidate for such an investigation is the initial state where N ↑ = N 2 spin-↑ and N ↓ = N 2 spin-↓ fermions are localized in the left and right wells respectively [see Fig. 1 s (x) denotes the Wannier state corresponding to the s ∈ {L, R} well and band b. Herein, we intend to address the dynamics of a system initialized in the state described by Eq. (3), especially focussing on the stability properties of the phase separation. Evidently, in the non-interacting case each one of the constituting particles will perform its individual tunneling oscillation with a frequency 2t b and, consequently, the phase separation imprinted in the initial state will be periodically lost and recovered during the time-evolution. However, in the case of g = 0 the individual tunneling channels of each of the particles couple due to the interparticle interaction. The interaction between the spin components is particularly important since the system accesses via tunneling, states possessing a substantial density overlap for anti-aligned spins yielding an interaction energy Fig. 1(c). Employing a mean-field argumentation one arrives at the conclusion that the tunneling among the wells slows down and eventually terminates as the repulsion increases. This is due to the large interaction energy of a spin-↑ and a spin-↓ atom occupying the same well when compared to the interaction energy contained in |Ψ(0) where the spin components are phase separated. However, the interparticle interaction possibly induces two-(or more) body correlations crucially affecting the dynamics of the system. As we shall demonstrate later on this is indeed the case and the dynamics for g = 0 is more involved than what is expected by the above-mentioned mean-field argumentation.

C. Magnetization Imbalance and MB Eigenstate Categorization in Terms of Bands
To monitor the degree of phase separation between the spin components during the dynamics of the system we employ the experimentally accessible measure [11] Here ρ is the spin-dependent, α ∈ {↑, ↓}, one-body density. Notice that both the Hamiltonian, Eq. (1), and the initial state, Eq. (3), are invariant under the transformation x → −x, | ↑ → | ↓ and | ↓ → | ↑ , implying that M ↑ + M ↓ = 0 is conserved during the dynamics. The quantity M takes its extreme values M = 1 and M = −1 when the particles within each of the wells are fully-polarized, a situation equivalent to a perfect phase separation. The sign of M in this case depends on whether the spin-↑ particles reside in the left (M = 1), as is the case for |Ψ(0) , or right (M = −1) well. In the case that M = 0 the spin-↑ and spin-↓ particles are distributed over both wells showing that the spin components are miscible. Since M = 0 corresponds to states magnetized along the x spatial-axis [see also Eq.(3)], M will be herewith referred to as magnetization imbalance.
Furthermore, let us note that for large barrier heights and weak or intermediate interactions, we expect that the band-gaps between the non-interacting bands constitute the largest energy scale of the system, see Fig. 1(a). As a consequence, the energetic characterization of the MB eigenstates in terms of non-interacting bands will be of great importance in the following. We assign each eigenstate of the non-interacting N -body system, |Ψ g=0 to an energetic class by employing the vector n B = (n 0 B , n 1 B , . . . ). This vector contains the occupation numbers of each of the non-interacting bands, , withn b sα being the number operator that counts the number of spin-α particles residing in the Wannier state φ b s (x). Accordingly, each eigenstate of the interacting system, g = 0, will be assigned to an energy class, n B , if it constitutes a superposition of non-interacting eigenstates of this particular class. For instance the initial state, |Ψ(0) , belongs to the n B = (2, 2, 0, . . . ) class for N = 4, see also Fig. 1(b). Indeed, the initial state for N = 4 contains two fermions in the 0th band (n 0 B = 2) and two additional ones residing in the 1st excited band (n 1 B = 2).

III. MANY-BODY EIGENSPECTRUM AND CORRELATED DYNAMICS
In this section we examine the eigenspectrum of the full MB HamiltonianĤ [see Eq. (1)] in the case of N = 4 fermions. Then we analyze the correlated dynamics of such a system initialized in the state |Ψ(0) [Eq. (3)] and subsequently left to evolve withinĤ. This investigation permits us to identify the emergent phase separation behavior between the spin components for varying interaction strength. To track the correlated dynamics of this system we employ ML-MCTDHX [26] and, in particular, its reduction for spin-1/2 fermions (for more details see Appendix C). ML-MCTDHX is an ab initio variational method that takes all correlations into account enabling us to reveal their influence into the static properties and in particular the dynamics of MB systems. We generalize our results to the N > 4 case in Appendix A.

A. Many-Body Eigenspectrum
The eigenspectrum ofĤ [Eq. (1)], for N ↑ = N ↓ = 2 fermions and varying g, is presented in Fig. 2(a), in the case of a relatively deep (V 0 = 8, w = 0.5) DW potential. The overlap of the MB interacting eigenstates, |Ψ i with the initial state, |Ψ(0) is indicated by the different colours in Fig. 2(a). Based on the eigenspectrum we can identify four different interaction regimes, indicated by A, B, C and B in Fig. 2(a), where the overlap of the initial state, |Ψ(0) with the MB eigenstates |Ψ i ofĤ exhibits a qualitatively different behaviour. In addition, by expanding each eigenstate |Ψ i in the number states of the Wannier basis, φ b s (x; t), (not shown here for brevity) we are able to identify its energetic class, n B (see section II C) which is important for identifying the interband processes emanating in the eigenspectrum and dynamics.
For weak interactions, g < 0.5 within the interaction regime A we observe that multiple eigenstates [the ones with E > 13.42 are hardly visible in Fig. 2(a)] contribute to the initial state. We remark that these states belong to the energy class n B = (2, 2, 0, . . . ) according to the energy categorization given in section II C. The energies of the eigenstates with E > 13.42 increases for increasing g, while their overlap with the initial state decreases, see 2(a) for g ≈ 0.2] emerge but overall the MB eigenspectrum is only slightly modified. These narrow avoided crossings result from the coupling of states belonging to the energy classes n B = (3, 0, 1, 0, . . . ) and n B = (2, 2, 0, . . . ) by a weak two-particle interband transfer process. Entering the interaction regime B, 0.5 < g < 2.5, we observe that the three eigenstates ofĤ possessing the dominant overlap with |Ψ(0) , are quasi-degenerate. In terms of increasing energetic order we refer to these eigenstates as |α , |β and |γ , see also the inset of Fig. 2(a). The existence of the quasi-degenerate predominantly occupied eigenstates within the B and also B (4.5 ≤ g < 5) interaction regimes implies that the time-scales of the dynamical evolution, which are associated with the energy differences of these quasi-degenerate states, are rather large. Therefore, these interaction regimes are very promising for studying the dynamical stability of the phase separation exhibited by the initial state |Ψ(0) . Note that the physical reasoning behind the emergence of this quasi-degenerate eigenstate manifold will be the main focus of section IV. At g ≈ 3.5 the three aforementioned quasi-degenerate eigenstates show a wide avoided crossing [see the dashed circle in Fig. 2(a)] with the eigenstates of the n B = (3, 1, 0, . . . ) energy class within the interaction regime C, 2.5 < g < 4. As we shall explicate later on, this interband avoided crossing is the fermionic analogue of the so-called cradle mode that has been identified in the interaction quench dynamics of spinless lattice trapped bosonic ensembles [43][44][45]. For larger repulsions, g > 4.5, the quasi-degeneracy of the predominantly occupied eigenstates reappears giving rise to the B interaction regime. The eigenspectrum for these interactions (g > 4.5) possesses a similar structure to the one observed within the interaction regime B. Note that the Tonks-Girardeau limit of our system is approached for g > 5 (not shown here for brevity). It is known [32] that in this case the eigenspectrum features an avoided crossing between the aforementioned quasi-degenerate states and the ones belonging to the energetically lowest class n B = (4, 0, . . . ). The eigenspectrum in this case can be theoretically described by using standard spin-chain techniques [32]. We remark that the state |γ possesses an interaction independent eigenenergy, associated with its fully antisymmetric character under particle exchange.

B. Correlated Dynamics
To inspect the stability of the phase separation encoded in |Ψ(0) for different interaction strengths we let the system initialized in the state |Ψ(0) of Eq.  ρ (1) α (x; t), see Fig. 3(a 1 ) and 3(a 2 ). Here, the dominant process is single particle tunneling. More precisely, the particles occupying the first excited band tunnel between the wells with much higher frequency than the particles occupying the lowest band. This can be identified by comparing the rate of tunneling of the two humped density structure (b = 1 band) appearing in Fig. 3(a 1 ) with the tunneling of the density residing near the center of the well (b = 0 band) within the time-interval 0 < t < 100. Additionally, an interaction induced dephasing, due to the involvement of the multitude of eigenstates identified in Fig. 2(a), is evident as |Ψ(0) does not completely revive during the time evolution.
Further monitoring the dynamical evolution of the system we observe that the phase separated state, |Ψ(0) [Eq.
(3)], is a long-lived metastable state within the interaction regime B. Indeed, for 0.5 < g < 2.5 we can infer the decay of the phase separation, imprinted in the magnetization imbalance, M and its subsequent revival [ Fig. 2(b)]. This process is relatively fast for weak interactions within the interaction regime B. For instance, the phase separated state |Ψ(0) decays to a miscible state with M = 0 at t ≈ 500 for g = 1. Qualitatively similar dynamics occurs but it is shown to be significantly slower for 1 < g < 4, e.g. at g = 2, M = 0 is reached for t ≈ 2000, while the life-time of |Ψ(0) exceeds t = 4000 for g ≈ 2.5 [ Fig. 2(b)]. The metastability of |Ψ(0) is accordingly well-justified since its life-times are much larger than the inverse of the characteristic tunneling rate of the ground, π/(2t 0 ) ≈ 42, and the first excited band, π/(2t 1 ) ≈ 190. To shed light into the dynamical evolution of |Ψ(0) we also inspect the one-body densities of the spin components, ρ are delocalized over both wells and they are almost perfectly overlapping which is in accordance to the value M = 0 [see Fig. 2(b)]. Note here that the absence of any signature of phase separation within each of the wells justifies the use of M as a measure of phase separation. Subsequently, the density of each component accumulates in the opposite well, than it was residing initially, but a small density portion remains in the initially populated well. Finally, at t ≈ 2700 an almost perfect revival of |Ψ(0) occurs. For larger evolution times, the above-mentioned dynamics is repeated in a periodic manner. Regarding the underlying tunneling mechanisms, the evolution of ρ (1) α (x; t) is indicative of a low-frequency two-body correlated tunneling dynamics for both spin components, as the entire density of two spin-aligned fermions seems to tunnel among the wells without being deformed. In addition, a contribution stemming from a single-particle tunneling process is also visible in Fig. 3(b 1 ) and 3(b 2 ), notice for instance the dynamics of the faint two-humped structure for t ≈ 600, t ≈ 1200 and t = 1800. In the following section it will be shown that the occurrence of the interaction regime B can be explained by examining the spin-order exhibited in the system.
For strong interactions, g > 2.5, the eigenstates belonging to the energy class n B = (3, 1, 0, . . . ) cross with the predominantly occupied eigenstates of the class n B = (2, 2, 0, . . . ) as shown in Fig. 2(a) at g ≈ 3.5. The states of the two energy classes exhibit two avoided crossings (indicated in Fig. 2(a) by the dashed circle) due to the interband interaction-induced coupling which is a manifestation of the cradle mode [43][44][45]. This resonant behaviour is directly imprinted on M , which shows a strong dependence of the lifetime of |Ψ(0) on the value of g, see Fig. 2(b) at g ≈ 2.8 and g ≈ 3.5. The spin-dependent one-body densities also show a tunneling behavior similar to the weakly interacting case, compare Fig. 3(c 1 ) and 3(a 1 ). The cradle mode is manifested as a dipole-like oscillation within each well. To explicitly demonstrate its existence we invoke the total one-body density fluctuations [43][44][45] defined as Indeed, δρ (1) (x; t) reveals dipole-like oscillations within both wells [see for instance Fig. 3(d) around t ≈ 400 i.e. the encircled region] and a beating dynamics for the intensity of the cradle mode. This beating can be understood by inspecting the eigenspectrum of the system [ Fig. 2(a)], where two almost perfectly overlapping cradle resonances can be identified at g ≈ 3.5, yielding two cradle frequencies of comparable magnitude. Notice that the cradle mode exhibited in our system is slightly different from its bosonic counterpart [43][44][45] as it does not involve overbarrier transport between the different wells but rather a direct interband population transfer within a particular well. The absence of overbarrier transport can be identified in Fig. 3(d) as the density fluctuations in the spatial region of the barrier, x ≈ 0 are vanishing.

IV. INTERPRETATION OF THE MAGNETIC PROPERTIES AND THE EFFECTIVE TIGHT-BINDING MODEL
Having appreciated the magnetic properties of the system within the fully-correlated ML-MCTDHX approach, we next proceed by constructing a reduced effective model. This model as we shall discuss below facilitates the qualitative interpretation of the correlated MB dynamics. In particular, the qualitative understanding of the underlying magnetic properties of the system via the effective model enables the identification of the decay mechanisms of the phase separation in a straightforward and intuitive way, allowing also, for comparisons with previous studies.

A. The Effective Tight-Binding Model
As already mentioned in section II, the band-gaps constitute the largest energy-scale of the system for both weak and intermediate interactions. It is therefore, well-justified to assume that a corresponding tight-binding model might sufficiently capture the observed dynamics. Within such a tight-binding model the Wannier states, φ b s (x), with s ∈ {L, R}, constitute the basis states of the MB Hamiltonian. The non-interacting Hamiltonian readŝ is the average energy of the non-interacting eigenstates forming the band, b. Also,â b † sα (â b sα ) is the operator that creates (annihilates) a spin-α particle in the Wannier sα . The exact form of the interaction term,Ĥ I , involves all matrix elements between the different Wannier states and it is, thus, quite complicated in appearance. Within the lowest-band approximation the Fermi-Hubbard model circumvents this issue by considering only on-site interactions and neglecting all density-induced tunneling effects [53]. It constitutes a valid approximation for large V 0 , where the underlying Wannier basis-states are well-localized to the corresponding wells. Additionally, g should define a sufficiently smaller energy scale than the band gap, ensuring that no significant interaction-induced interband tunneling, such as the cradle mode, occurs. Fermi-Hubbard models have been very successful in describing various effects emanating in a variety of settings where DW or lattice potentials are involved [54,55].
Therefore, it is tempting to approximate the exact interaction term,Ĥ I , by the following effective onê where to the inter and intraband on-site interactions respectively. However, as it can be easily verified the last term of Eq. (6) breaks the SU(2) symmetry of H [Eq. (1)], since it does not commute withŜ 2 . In order to avoid this artificial symmetry breaking one needs, also, to include into the effective Hamiltonian the term The termĤ exc I , which is present in the exactĤ I of Eq. (1), incorporates the effect where two fermions in different bands but on the same well can exchange their spin due to their mutual interaction. Models that extend the Hubbard model in a similar manner to the above-mentioned have been employed in the context of the metal-insulator transition emanating in d-electron systems, for a review see [56].
Including all of the above-mentioned terms into an effective tight-binding Hamiltonian results in the following multi-band tJU modelĤ sβ , k ∈ {x, y, z}, s ∈ {L, R} and i, j, k refer to the unit vectors in spin-space andn b s =n b s↑ +n b s↓ . tJ models, where the on-site interaction term vanishes as double site occupations are adiabatically eliminated, have been originally employed to describe magnetic phenomena in condensed matter physics [57][58][59] and later for the interpretation of some aspects of superconductivity [60][61][62]. Physically, the effective Hamiltonian of Eq. (8) 8)] is that gU b /E b 1 allowing for the interaction-driven interband processes to be safely neglected. Within this model states of different energy classes, n B do not couple and as a consequence all the elements of n B are conserved. As we have previously established within the full MB system (that does not possess this symmetry) such interband effects do not alter the eigenspectrum significantly within the interaction regimes A and B.
Below we argue why this model leads to a metastable, phase separated, state |Ψ(0) , in the case of intermediate repulsions, qualitatively explaining the magnetic order exhibited within the interaction regime B.

B. Magnetic Properties of the Effective Model
Let us first discuss the relevant properties of the N -body eigenspectrum of the tJU model. We operate in the t b /(gU b ) 1 limit, where we can neglect the tunneling term ∝ t b . Indeed, for the system examined in section III the criterion gU b t b 1 is well-satisfied 1 within the interaction regime B, 0.5 < g < 3. In view of the decoupling of different energy classes n B within the effective tJU model we will focus on the particular energetic class that the initial state, |Ψ(0) , belongs to, namely, n B . This class is defined as n b B = 2 for b < N/2 and n b B = 0 otherwise 2 . Focusing on the simplest case of t b = 0, for all involved b, the effective Hamiltonian can be expanded in two intra-well Hamiltonian termsĤ eff =Ĥ R +Ĥ L that are decoupled. among them. By projecting these Hamiltonian terms into the energy class n B the former readŝ whereP B is the projection operator into n B . Equation (9) corresponds to a ferromagnetic Heisenberg model, incorporating additional energy shifts depending on the particle occupation ∝n b s , s ∈ {L, R}. For g > 0 the sum of these energy shifts contained withinĤ L andĤ R is minimized in the case that no double occupations of a particular site occur, i.e. Ψ|n b s↑n b s↓ |Ψ = 0 for all b and s. The spin configuration for t b = 0 can be characterized by the quantum numbers S, S L , S R , whereŜ 2 s , refers to the total spin within the s ∈ {L, R} well. It is well-known that ferromagnetic Heisenberg models exhibit ferromagnetic ground states [63] and as a consequence the ground states of Eq. (9) correspond to the largest possible values of S L and S R . Notice, also, that |Ψ(0) is characterized by maximum S L and S R , since the spin-state within each well is fully spin-polarized. As a consequence, we can conclude that |Ψ(0) belongs to a degenerate manifold of dimension N/2 at an energy E = E B = 2 N/2−1 b=0 b . This manifold consists of the states |Φ(t b = 0); S with quantum numbers S L = S R = N 4 but varying total spin S ∈ {0, 1, . . . , N 2 } (see also below). In addition, the eigenstates |Φ(t b = 0); S get energetically well-separated from the other states with n B = n B as the gap between them scales linearly with g, see Eq. (9).
The inclusion of the tunneling term for t b = 0 induces couplings between the above-mentioned degenerate states resulting in the lifting of their degeneracy. Indeed, by treating the tunneling term in Eq. (8) within second order perturbation theory [see Appendix C], we can show that in the t b gU b limit the effective Hamiltonian projected on the manifold of degenerate states spanned by |Φ(t b = 0); S readŝ whereP D is the projection operatorP D = b=0 J b0b . Equation (10) provides great insight into the structures emanating in the intra-well ferromagnetically correlated states within the tJU model in the case of non-vanishing tunneling. Indeed, the inclusion of tunneling for t b = 0 results to an apparent antiferromagnetic Heisenberg exchange interaction for g > 0, known as the Anderson kinetic exchange interaction [42]. Note that the total spinŜ =Ŝ L +Ŝ R commutes withP DĤeffPD implying that the eigenstates of the tJU model |Φ; S reduce within the zeroth order approximation to the ones for t b = 0, i.e. |Φ; S = |Φ(t b = 0); S + O( t b gU b ). Regarding their eigenenergies, the tJU eigenstates, |Φ; S are expected to be energetically ordered in terms of increasing S due to the antiferromagnetic character of the Anderson exchange interaction and be quasi-degenerate possessing energy shifts among them of the order of Ω The above properties of the tJU eigenspectrum imply that the initial state |Ψ(0) being a superposition of the eigenstates |Φ; S dephases during the time-evolution with a slow timescale ∼ min b (π/Ω b d ). In particular |Ψ(0) can be expanded in terms of these eigenstates by utilizing the Clebsch-Gordan coefficients leading to Φ; Moreover, the maximum values of S L = N 4 and S R = N 4 which characterize the eigenstates |Φ; S ≈ |Φ(t b = 0); S imply intra-well ferromagnetic correlations for particles occupying the same well 3 and stem from the ferromagnetic Hund exchange interactions emanating in Eq. (8). Therefore, the emergence of the interaction regime B can be attributed to the dominant contribution of the intra-well ferromagnetic correlations when compared to the above-mentioned effective antiferromagnetism stemming from t b , see Eq. (10). We remark here that the nature of this effective antiferromagnetism has been identified and studied by employing spin-chain models tailored to operate in the vicinity of the Tonks-Girardeau limit, g → ∞ [27][28][29][30][31][32]. Note also that this notion of antiferromagnetism does not conflict with our notion of ferromagnetism as the first is an effective magnetic phenomenon induced by the tunneling, t b , while the second is a result of the exchange interaction term in Eq. (7).

C. Comparison with ML-MCTDHX
Before analyzing further the magnetic properties of the system and their connection to the emergent tunneling dynamics let us first establish that the magnetic properties exhibited in the framework of the tJU model carry forward to the fully correlated case. To this end we shall compare the eigenspectra obtained within the tJU model with the ML-MCTDHX method.
The relevant eigenenergies within the tJU model [Eq. (8)] appear in Fig. 4(a) as colored lines referring to the case N = 4. Here the three eigenstates |Φ; S , with S = 0, 1, 2 possess three distinct eigenenergies at g ≈ 0. The energy difference between the S = 0 and S = 1 states is given by t 0 and the one between the S = 1 and S = 2 corresponds to t 1 . This decrease of the single-particle energy of |Φ; S = 1 and |Φ; S = 0 stems from the occurrence of one and two doublons respectively in the g ≈ 0 case. The formation of these doublons implies the double occupation of the For increasing g the energy of the S = 1 and S = 0 eigenstates is larger due to the involvement of these doublons which contribute a substantial amount of interaction energy. Most importantly, for 0.5 < g < 2.5 (interaction regime B) the energies of the eigenstates |Φ; S converge towards the eigenenergy of |Φ; S = 2 , possessing E S=2 (g) = E B ≈ 13.424, and this leads to the formation of the quasi-degenerate manifold, identified also in Fig. 2(a). It can also be checked that the energy differences between the states |Φ; S are consistent with Eq. (10) possessing a characteristic energy scale of Ω 1 d ≈ 0.033/g. Figure 4(a) further reveals that the eigenstates of the tJU model follow closely the behaviour of the eigenstates of the MB system, represented as dots in Fig. 4(a), within both the interaction regimes A and B. There are a few discrepancies associated with the avoided crossings emerging in the interaction regimes A and C which, as also mentioned in section III, stem from the couplings between states with different n B . Such couplings are indeed neglected within the tJU model. Nevertheless, the agreement within the interaction regime B is almost perfect and it can be further shown that the key ingredients of the magnetic order within the tJU model are also exhibited within the fully correlated case. Indeed, the overlaps between the MB eigenstates |α , |β , |γ and the initial state, |Ψ(0) agree well with those found within the effective tJU description [see Fig. 4 Fig. 4(b) we demonstrate the large overlap of the MB eigenstates |α , |β and |γ with the eigenstates, |Φ(t b = 0); S , of the tJU model for t b gU b within the interaction regime B. Indeed, this overlap is in excess of 95% [see also the inset of Fig. 4(b)], a result which is also consistent with the values obtained within the tJU model for the overlaps | Φ; S|Φ(t b = 0); S | 2 (not shown for brevity). The above mentioned findings explicitly showcase that the magnetic order exhibited in the interaction regime B within the tJU model carries forward to the MB case. However, for stronger interactions and as the interaction regime C is approached, e.g. see Fig. 4(b) at g ≈ 2.5, the overlap of the MB eigenstates and the |Φ(t b = 0); S states decreases. This feature is beyond the tJU model description and occurs due to the interband coupling caused by the presence of the cradle mode. 3 Note that within this particular configuration in terms of n B and n D , the total spin within the s well solely depends on the corresponding spin-spin correlator,PŜ 2 sP = N

D. Relation of the Magnetic Properties to the Tunneling Dynamics
Having identified the magnetic order of the interaction regime B within the full MB approach by comparing to an effective model, we subsequently showcase the relation of these magnetic properties to the tunneling dynamics of the system, see also Fig. 3(b 1 ) and 3(b 2 ). To unravel this interplay we define the states with S L = S R = N 2 and a definite spin projection S z;s = b S b z;s within each of the wells, namely HereŜ ±;s = bŜ b x;s ± iŜ b y;s refer to the spin increasing and lowering operators within the s ∈ {L, R} well. Note that it can be verified that the states |Ψ; N ↑L , N ↑R are related to the states |Φ(t b = 0); S by a unitary transformation. But in contrast to |Φ(t b = 0); S the states |Ψ; N ↑L , N ↑R have a definite number of spin-↑ and spin-↓ particles within each well. The expansion of |Ψ; N ↑L , N ↑R in the basis |Φ(t b = 0); S can be easily obtained by employing the Clebsch-Gordan coefficients [64]. The introduction of this new basis relates the phenomenon of quasidegeneracy of the states |Φ(t b = 0); S exhibited both within the tJU model (|Φ; S ≈ |Φ(t b = 0); S ) and the full MB case (|α ≈ |Φ(t b = 0); S = 0 , |β ≈ |Φ(t b = 0); S = 1 and |γ ≈ |Φ(t b = 0); S = 2 ), see Fig. 4(b), with the emergent tunneling processes. Owing to the unitary transformation between the states |Ψ; N ↑L , N ↑R and the approximate eigenstates |Φ(t b = 0); S , the accumulation of relative phases between the eigenstates during the dynamics [due to their different eigenenergies, see Fig. 4(a) and Eq. (10)], corresponds to a population transfer between the |Ψ; N ↑L , N ↑R states and hence to an apparent tunneling dynamics within the spin components. For the particular case of N = 4 particles this mechanism is illustrated in Fig. 5. Notice that due to the strongly-correlated nature of the involved states [see also Fig. 5] such a mechanism is absent within the Hartree-Fock mean-field theory since |Ψ; N ↑L = 1, N ↓R = 1 cannot be written as a single Slater determinant.
To explicitly demonstrate the occurrence of this tunneling mechanism we present in Fig. 4(c) the overlap of the time-dependent wavefunction, |Ψ(t) , with the states |Ψ; N ↑L , N ↑R for g = 1 (interaction regime B). For 0 < t < 500 we observe a population transfer process from the initial state |Ψ(0) = |Ψ; 2, 0 to the states |Ψ; 1, 1 and |Ψ; 0, 2 . Recall that these two processes have, also, been identified in the time evolution of ρ (1) σ (x; t) [see Fig. 3(b 1 ), 3(b 1 ) and section III]. Most importantly, the intricate relation of the tunneling dynamics to the magnetic properties of the system is now evident via employing the unitary transformation connecting the states |Ψ; N ↑L , N ↑R to the eigenstates |Φ(t b = 0); S . For later times t ≈ 2600 an almost perfect revival of the state |Ψ; 2, 0 is exhibited owing to the commensurability of the frequencies of these particle transfer processes. Indeed, the two-body tunneling process |Ψ; 2, 0 ↔ |Ψ; 0, 2 is found to possess a roughly three times smaller frequency than the single-particle tunneling mode |Ψ; 2, 0 ↔ |Ψ; 1, 1 [see Fig. 4(c)].

E. Spin-Spin Correlations
Having described in detail the interconnection of the magnetic properties of the system and its tunneling dynamics we are able to shed light onto the relation of ferromagnetism and phase separation. The ferromagnetic order of a Fermi gas is characterized by the spin polarization and the spin-spin correlations of the system. The total spin polarization and the total spin of the system, the latter being related to the spin-spin correlator [34], are constant during the dynamics due to the symmetries of the Hamiltonian [Eq. (1)]. As a consequence no global ferromagnetic order can appear during the dynamics due to the conservation laws stemming from the above symmetries. However, as the tJU model reveals the intra-well magnetic properties are important for the adequate description of the system. In this spirit, the quantity M besides being a measure of the phase separation also quantifies the spin polarization within each well [Eq. (4)]. An adequate quantity that captures the intra-well magnetic correlations is also hinted by the effective tJU model. This refers to the total spin within each of the wells, Ψ(t)|Ŝ 2 s |Ψ(t) , with s ∈ {L, S}. In particular, we can employ a more refined quantity by involving some of the magnetic properties of the system identified within the tJU model. As we have previously discussed, the subset of states |Φ(t b = 0); S are characterized by ferromagnetic spin-spin correlations within each well since S L = S R = 1. Specifically, |Φ(t b = 0); S are the only states within the configuration n B = (2, 2, 0, . . . ) that exhibit this property (see also section IVB). It is thus instructive to evaluate the overlap of the MB wavefunction, |Ψ(t) with the states |Φ(t b = 0); S , i.e. C F F = 2 S=0 | Φ(t b = 0); S|Ψ(t) | 2 . C F F is an adequate quantity for studying the spin-spin correlation properties of the system, as it constitutes a lower bound for the values of the intra-well spin-spin correlator Ψ(t)|Ŝ 2 L |Ψ(t) ≥ 2C F F and Ψ(t)|Ŝ 2 R |Ψ(t) ≥ 2C F F . Accordingly, large values of C F F indicate that both wells are simultaneously characterized by intra-well ferromagnetic spin-spin correlations.
The time evolution of C F F is presented in Fig. 4 (d) for varying interaction strength g. For weak interactions, within the interaction regime A, C F F exhibits rapid fluctuations between zero and unity manifesting the periodic decay and revival of the intra-well ferromagnetic spin-spin correlations of the initial state, |Ψ(0) . Recall that the phase separation, and hence the intra-well spin polarization, is unstable in this interaction regime exhibiting decay and revival oscillations, see also Fig. 2(b), 3(a 1 ) and 3(a 2 ). However, in the interaction regime B, we observe that the spin-spin correlations within each well are ferromagnetic since C F F = 1. Indeed, the value of C F F is almost constant and possesses a large value being of the order of C F F ≈ 0.98, see Fig. 4(d). The weak fluctuations of C F F around this average value, further, showcase the stability of the intra-well ferromagnetic order. Note also that the intra-well spin polarization quantified by M is characterized as metastable within the interaction regime B. Entering the interaction regime C (2.5 < g < 4.5) we observe that C F F exhibits multi-frequency oscillations. These oscillations can be explained in terms of the observed resonance of the cradle mode which introduces an interband coupling 4 , see also Fig. 2(a) and 3(d). For even stronger interactions, i.e. within the interaction regime B , the intra-well ferromagnetic order is reestablished and it is characterized by large and almost constant values of C F F during the dynamics.
The above results manifest the close relation between the intra-well ferromagnetic order emanating in a DW trap and the ferromagnetic order emerging in a harmonic trap with weakly broken SU(2) symmetry, as it has been demonstrated in Ref. [34]. This order appears for intermediate interactions where the ferromagnetic Hund exchange interaction, stemming from spin-exchange interaction processes, see e.g. Eq. (7), dominates and leads to largely stable ferromagnetic spin-spin correlations but a fluctuating polarization. The different imposed potential alters the manifestation of this ferromagnetic order during the dynamics. In the case of the DW the ferromagnetic order is exhibited locally within each of the wells and as discussed above implies a metastable phase separated state for the system. In contrast, in the case of a harmonic trap the emergent ferromagnetic order affects the global values of the spin polarization and total spin implying miscibility of the contributing spin components [34]. This difference stems from the Hund exchange interaction between two fermions which is only sizable if the involved single particle states (i.e. the orbitals) possess a significant density overlap, see also Eq. (7). To understand this analogy further in the following section we study the dynamics of the DW by employing an additional potential that breaks the SU(2) symmetry of the system.

V. SU(2) VIOLATING CASE
Up to this point, we have identified the metastability of the phase separated state in a DW due to the presence of the SU(2) invariance of the system. Also we have characterized the emerging metastability of the phase separated state appearing for intermediate interactions and connected it to the magnetic properties of the system. Next we aim to show that the phase separated state is stable within region B even in the case that the SU(2) symmetry is weakly broken. Also, in analogy to Ref. [34] the intra-well ferromagnetic correlations are shown to persist within this region. Moreover, the implications regarding the magnetic order exhibited in a DW are briefly discussed.
To study the case of a system with broken SU(2) symmetry we employ a linear gradient of the magnetic field which shifts the energies of the spin-↑ and spin-↓ fermions in a spatially-dependent manner. The corresponding term which is incorporated in the MB Hamiltonian of Eq. (1), readŝ The value of B 0 determines the energy offset between the two wells for the different spin components. A positive value of B 0 means that it is energetically preferable for the spin-↑ atoms to occupy the right-well and the spin-↓ atoms the left-well. Accordingly, when B 0 < 0 it is favorable for the spin-↑ and spin-↓ atoms to occupy the left and the right-well respectively.  Figure 6 illustrates the time-evolution of the magnetization imbalance M and C F F which quantifies the degree of intra-well ferromagnetic correlations, for varying B 0 at three different values of g corresponding to the interaction regimes A, B and C. We observe that in the weakly-interacting case [belonging to the interaction regime A in Fig.  6(a) and 6(b)] and for B 0 < 0 both M and C F F are stable throughout the time-evolution indicating that the system remains close to its initial state. For B 0 > 0 a multitude of resonances appear at different intervals of B 0 involving prominent tunneling as captured by M [ Fig. 6(a)]. Also, C F F reveals that the state of the system is driven away from the S L = S R = 1 manifold [ Fig. 6(b)] since C F F < 1. These resonances correspond to possible tunneling pathways where the spin-↑ particles occupying initially the left-well of the DW resonantly tunnel to the right-well (or to the opposite direction for the spin-↓ atoms) leading to the decay of the intra-well ferromagnetic order.
For g = 2 (interaction regime B) it can be deduced that besides the very narrow region around the SU(2) symmetric case, i.e. at B 0 = 0, the phase separated initial state is stable for |B 0 | < 0.04 as M (t) ≈ 1 throughout the evolution [see Fig. 6(c)]. Notice also that within these values of |B 0 | the ferromagnetic intra-well order is stable as indicated by C F F (t) ≈ 1 [ Fig. 6(d)]. The stable phase separated state appears due to the quasi-degeneracy of the states |α , |β , |γ in the SU(2) preserving case for the interaction regime B. As stated in the previous sections these states, owing to their intra-well ferromagnetic correlations, lie in an energy region of the MB spectrum where no other eigenstates appear and are quasi-degenerate characterized by a different value of the total spin S [see also Fig. 4(a)]. Recall that these states possess C F F ≈ 1 indicating their intra-well ferromagnetic character. Moreover, their energetic ordering in terms of increasing S manifests the presence of the weak antiferromagnetic Anderson exchange interaction, [see section IV B and Eq. (10)]. By breaking the SU(2) symmetry with the additional spin-dependent potential described by Eq. (12) the states |α , |β and |γ couple with one another resulting in the formation of eigenstates with definite number of spin-↑ and spin-↓ atoms in each of the wells (results not shown here for brevity). Therefore, for decreasing B 0 < 0 the initial state, |Ψ(0) , becomes the lowest-in-energy state with S L = S R = 1, while it corresponds to the highest-in-energy eigenstate of the same manifold of states for B 0 > 0. In both cases the phase separation of this state is stable as imprinted also in the time evolution of M (t) for |B 0 | < 0.04 [see Fig. 6(c)]. In the vicinity of B 0 ≈ 0, M (t) is depleted during the time-evolution while C F F (t) ≈ 1 throughout the dynamics. The appearance of this region is explained by the fact that the couplings between the states |α , |β and |γ associated withĤ g are smaller than their energy differences due to the Anderson kinetic exchange interaction (being of the order of t b gU b ). The latter implies a large but finite life-time of the phase separation of the initial state, in agreement with the SU(2) preserving case B = 0. In addition, further resonances appear when |B 0 | > 0.04 for g = 2 [see Fig. 6(c)]. More specifically, the resonances at B 0 > 0.04 correspond to tunneling resonances in a similar fashion to the case of the interaction regime A [ Fig. 6(a)]. The positive shift of these resonances when compared to the corresponding ones appearing for g = 0.5 is attributed to the increased interaction energy of the states accessed by tunneling. For B 0 < −0.04 another set of resonances occurs in Fig. 6(c) that correspond to interband processes similar to the aforementioned cradle mode. These resonances emerge due to the coupling of different bands induced by the interactions.
Within the interaction regime C the stability properties of the phase separation are similar to the corresponding ones of the interaction regime A, compare in particular Fig. 6(e) and 6(f) to Fig. 6(a) and 6(b) respectively. For large B 0 < 0 the initial state is stable [see Fig. 6(e) and 6(f)], however, for B 0 ≈ 0 the phase separation and intra-well ferromagnetic order as imprinted in M (t) and C F F (t) respectively fluctuate during the dynamics. This fluctuating behaviour can be explained by the inter-band coupling that occurs within this interaction regime suppressing the intra-well ferromagnetic order of the initially phase separated state [see also  Fig. 2(a) for g ≈ 3 is observed.
The above discussed stability properties of the phase separated state, especially within the interaction regime B, provide direct insight into the magnetic properties of the SU(2) violating system. First, the fact that the phase separated state, |Ψ(0) , which is not an eigenstate ofŜ 2 , becomes an eigenstate of the system,Ĥ =Ĥ 0 +Ĥ I +Ĥ g , even for a relatively small breaking of the SU(2) symmetry shows that, as also identified previously, for a DW there is no global ferromagnetic order imprinted in S 2 . This is in contrast to the case of the harmonic confinement as it has been demonstrated in Ref. [34]. Instead, for a DW trap the instability of the S 2 becomes more pronounced for intermediate interactions. This property can be understood by inspecting the effective tJU model [Eq. (8)]. For fermions confined in a DW, ferromagnetic Hund interactions occur only between particles that reside in the same well and as a consequence only the intra-well ferromagnetic correlations are robust within each well. An observation that is also supported by the apparent stability of the phase separated state except for the cases within the interaction regimes B and C where inter-band couplings are involved, see Fig. 6(b), 6(c), 6(d) and 6(f). Most importantly, for intermediate interactions supporting the intra-well ferromagnetic order [see Fig. 6(c)] the phase separated state is stabilized even for a very weak breaking of the SU(2) symmetry. This feature of the DW system can be understood by the fact that the extremely weak Anderson kinetic exchange interaction is the only magnetic mechanism that can possibly prohibit the coupling of states with different S for a system with broken SU(2) symmetry. On the contrary, the intra-well ferromagnetic order is stable independently of whether the SU(2) symmetry is preserved or it is weakly broken as the intra-well ferromagnetic correlations are protected by the much stronger Hund exchange interaction. The above imply that within the interaction regime B dominated by ferromagnetic intra-well correlations an instability occurs which is triggered by the breaking of the SU(2) symmetry. This instability leads to the formation of two polarized ferromagnetic domains of the spin components as the system phase separates almost perfectly among the two wells.

VI. CONCLUSIONS
We have explored the stability of the phase separated state of interacting spin-1/2 fermions confined in DW potentials. Most importantly, we have revealed an interaction regime characterized by a metastable phase separation for moderate interactions. By invoking an effective tight-binding model, we unveil that the metastability of the phase separation is related to the formation of a quasi-degenerate manifold of states described by ferromagnetic intra-well spin-spin correlations but varying total spin. The formation of this quasi-degenerate manifold of states can be intuitively understood by the inclusion of an effective ferromagnetic Hund interaction, stemming from the spin exchange interaction between two interacting particles residing at the same well. This exchange interaction cannot be neglected due to the large spatial overlap of the particles occupying different bands but the same well of the DW. The breaking of the SU(2) symmetry is found to substantially alter the behaviour of the system in this interaction regime where the ferromagnetic correlations dominate. Indeed, the phase separated state becomes stable even when we break the SU(2) symmetry by employing a very weak linear gradient of the magnetic potential.
The description of the magnetic properties of 1D fermions in terms of the ferromagnetic Hund interaction provides a unifying viewpoint on the relation between phase separation and ferromagnetism. Most importantly, it provides a theoretical framework via which the stability of ferromagnetic correlations in the absence of SU(2) symmetry (see also [34]) can be understood. In particular, the ferromagnetic correlations are found to be stable only within the spatial regions where the Hund interaction is strong, i.e. within each of the wells of a DW and not between them. In this picture the ferromagnetic correlations of the system are not directly related with the phase separation in contrast to the conventional Stoner instability viewpoint. Instead, the effective antiferromagnetism induced by the Anderson kinetic exchange interaction is responsible for the absence of phase separation in SU(2) symmetric systems. Indeed, when this effective antiferromagnetism is weak the system is found to be unstable towards phase separation. More precisely, in the case of a DW potential these two phenomena are indeed related. In the interaction regime where the ferromagnetic Hund interaction dominates the Anderson kinetic exchange interaction leading to stable intra-well ferromagnetic correlations, even a weak breaking of the SU(2) symmetry enforces the system to phase separate.
Our work sets several avenues of further study that can be pursued. First, notice the absence of any obvious limitation of the underlying mechanisms that would make them incapable of describing higher dimensional settings. The examination of higher dimensional settings is therefore a promising next step for understanding the ferromagnetic properties emerging in DW systems. Also, the tunability of the phase separation by weakly breaking the SU(2) symmetry gives rise to the prospect of controlling the formation of ferromagnetic domains in the case of DW or lattice systems. Finally, the inclusion of various inherent effects that break the SU(2) symmetry of a Fermi system such as spin-orbit coupling or weak spin-dependent interactions might allow cold atoms to form realistic models that better emulate the ferromagnetic properties encountered in real materials. The discussion in section IV B reveals that within the effective model description an overall similar dynamical behaviour of the system is expected independently of the particle number [see also footnote 2]. To verify this expectation within the fully correlated approach we investigate the dynamical behaviour of a system consisting of N = 6 fermions and identify the underlying phenomenology associated with the different interaction regimes A, B and C in the corresponding spin-α one-body densities, ρ (1) α (x; t), illustrated in Fig. 7. In particular, for weak interactions (g = 0.05) the one-body density of both spin-↑ and spin-↓ fermions exhibits a tunneling dynamics among the wells, see Fig. 7(a 1 ) and Fig. 7(a 2 ). In this case, each of the particles occupying the three energetically lowest bands performs an individual tunneling oscillation with a frequency close to the one associated with the band it occupies, t b , see for instance the fast tunneling of the three-humped structure emerging in ρ (1) α (x; t) in comparison to the overall slower tunneling dynamics. This observed dynamics is in line to the one emerging within the interaction regime A for N = 4 particles [compare Fig. 3(a) with Fig. 7(a)]. For increasing interactions, g = 1.75, no tunneling oscillations are observed and the phase separation appears to be almost completely stable within the time scales we have studied, see Fig. 7(b 1 ) and 7(b 2 ). This behaviour of the one-body density is characteristic for the interactions belonging to the regime B, where as identified in the N = 4 particle case the tunneling dynamics slows down dramatically [see also Fig. 3(b)] as a consequence of the formation of the quasi-degenerate manifold of eigenstates with ferromagnetic intra-well correlations. Finally, the cradle mode being the characteristic feature of the interaction regime C [see also Fig. 3(c) and 3(d)] can also be observed for N = 6. Inspecting the dynamics of the one-body density for g = 3.5 [ Fig. 7(c)], we observe a collective tunneling mode of the density among the wells, as well as, deformations of the one-body density within each of the occupied sites possessing a much larger frequency than the tunneling mode. By employing the temporal fluctuations of the total one-body density, δρ (1) (x; t) [see Fig. 7(d)] these deformations can be related to the emergence of the cradle mode, verifying the existence of the interaction regime C in the N = 6 case. As we have discussed in the main text (see section IV E and V) the relation of the phase separation phenomenon and the ferromagnetism depends on the shape of the external potential imposed on the atoms. Indeed, it is found that despite the fact that the same microscopic mechanisms are at play for a parabolically or a DW trapped spin-1/2 Fermi system, the manifestation of the above-mentioned phenomena differs significantly. The purpose of this section is to study the dependence of the stability properties of the phase separated state, |Ψ(0) , Eq. (3) on the barrier height of the DW potential. To achieve this we study the case of a shallower DW with V 0 = 5 and w = 0.5 and compare with the case of V 0 = 8.
The eigenspectrum for a shallow DW is presented in Fig. 8(a). The qualitative structures emerging in the eigenspectrum for V 0 = 5 are similar to the case of V 0 = 8 [compare Fig. 8(a) with Fig. 2(a)]. However, there are also prominent quantitative differences which, as we shall explain below, lead to a different dynamical behavior. Within the regime A, 0 < g < 1, the role of eigenstates with high energy (see Fig. 8(a) for E > 12.2 and g < 1) is very pronounced as they accumulate a population larger than in their deep DW counterpart, see also Fig. 2(a). In the dynamics of the shallow DW this translates to a much faster loss of M (see e.g. Fig. 8(b) for g < 1) when compared to the case of the deep DW [ Fig. 2(b)] which is accompanied with the loss of intra-well ferromagnetic correlations imprinted in C F F , see Fig. 8(c). Of course, this difference is simply caused by the larger tunneling rates, t b involved in the V 0 = 5 case [ Fig.  1(a)]. The differences between the two setups become more interesting in the intermediate interaction regime, B, for 1 < g < 2.5. In the shallow DW case three eigenstates dominate similarly to V 0 = 8, but their spacing is quantitatively larger in the shallower DW [compare Fig. 8(a) with Fig. 2(a) for g ≈ 2]. This is not surprising since the spacing of these eigenstates (see also section IV B) is proportional to t b gU b which decreases with increasing V 0 . In addition, and in direct contrast to the V 0 = 8 case higher-lying eigenstates [see Fig. 8(a) for 1 < g < 2.5 and E > 12.3] and most importantly lower-lying ones [see Fig. 8(a) for 1 < g < 2.5 and E > 11.8] are involved in the dynamics within this regime. Accordingly the dynamics of M and C F F shows that in the shallow DW case the initial state cannot be characterized as metastable for any interaction in the regime B. Indeed, the magnetization imbalance M (t) [ Fig.  8(b)] is greatly suppressed for t > 100 possessing values M (t) < 0.6 for all interactions in 1 < g < 2.5. Regarding the spin-spin correlations it can be seen that C F F (t) is almost stable during the dynamics except for a very fast decay at initial times t < 4 [see Fig. 8(c)]. During the time-evolution it acquires values of the order of C F F ≈ 0.8, for all interactions within the regime B showcasing predominantly ferromagnetic intra-well correlations. The above implies that while the mechanisms at play in the shallow DW case are similar to the ones emerging in the case of a deeper DW, the apparent phenomenology is altered due to the pronounced involvement of lower-lying states. These lower-lying states are able to alter the dynamics within the regime B because, as it can be seen by inspecting the eigenspectrum for g ≈ 3 the cradle resonances are much wider in the case of a shallower DW thus affecting a broader interaction regime than for V 0 = 8.
In the case of V 0 = 5 the regime C appears in the interaction range 2 < g < 4.5. The phenomenology taking place within C is completely analogous to the case of V 0 = 8. Indeed, the tunneling is prevalent within this regime as imprinted in the fluctuating behavior of the magnetization imbalance M (t) [see Fig. 8(b)]. In addition, the intra-well spin-spin correlations imprinted in C F F (t) can be also seen to fluctuate similarly to the case of V 0 = 8 [compare  Fig. 8(a)] it can be deduced that in this regime a quasi-degenerate manifold of the predominantly occupied eigenstates begins to form similarly to the regime B encountered for V 0 = 8 .
In conclusion, the nature of the microscopic mechanisms that govern the stability properties of phase separation are not altered as the depth of the DW changes. However, because of their direct competition, in particular between the exchange interaction and the combined effects of the tunneling and the cradle mode, the observed dynamics differs significantly as the barrier height, V 0 decreases. Indeed, the mechanisms competing with the exchange interaction become more prevalent for a shallower DW as it is also clearly imprinted in the corresponding eigenspectrum. This renders the intra-well ferromagnetic order unable to completely dominate the dynamics for every interaction strength, resulting in the absence of stable ferromagnetic intra-well correlations and its direct imprint on the dynamics i.e. the metastability of the phase separation.

Appendix C: Anderson Kinetic Exchange Interaction
The purpose of this section is to provide the explicit derivation of the effective antiferromagnetic interaction acting upon the different wells of our DW setup. This antiferromagnetic interaction is similar to the Anderson kinetic exchange interaction emanating among the different sites of a lattice within the Hubbard model [42]. Although such an effective magnetic term can be derived within the Resolvent formalism by invoking less assumptions, here we opt to employ the standard Reyleigh-Schrödinger second order perturbation theory due to its mathematical (and physical) clarity.
The terms appearing in the Hamiltonian of the tJU model [Eq. (8)] can be separated into two Hamiltonian terms 5 that solely act within each of the wells,Ĥ s , with s ∈ {↑, ↓}, and a Hamiltonian part corresponding to the coupling between them,Ĥ LR . By performing this separation the effective Hamiltonian readsĤ eff =Ĥ R +Ĥ L +Ĥ LR . The intra-well Hamiltonian terms, correspond to ferromagnetic Heisenberg models with additional occupation dependent terms ∝n b s . The intra-well describes the tunneling among the wells. Our intention is to perturbatively treatĤ LR and show that it acts as an effective antiferromagnetic interaction between the particles occupying the same band but different wells. According to the discussion in sections II B and IV B we are particularly interested in the configuration with no doublons and a single occupation of each Wannier state up to the b = N 2 − 1 band. The projection ofĤ s to this particular configuration results in the Heisenberg model which possesses the degenerate ground states Note here that the action ofĤ LR on the basis of the ground state manifold |Ψ; N ↑L , N ↑R is rather simple due to its product state character. IndeedĤ LR |Ψ; N ↑L , N ↑R can be expressed via the action of the creation and annihilation operators on the single-well ferromagnetic states | Ns 2 , N ↑s − Ns 2 s . Indeed the annihilation operator creates a vacancy to the single-well ferromagnetic stateŝ whereŨ b0 = N 2 −1 b=0 J b0b (recall that U b0 = J b0b0 ), and as a consequence also constitute eigenstates ofĤ s . By using Eq. (C5) and (C7) we can show that each term appearing inĤ RL couples each state |Ψ; N ↑L , N ↑R to a single state containing one doublon. This coupling scheme is schematically depicted in Fig. 9. Here the state |Φ Therefore, by employing the basis states |Ψ; N ↑L , N ↑R and |Φ b0 s ; N ↑L , N ↑R the couplings between states possessing no and one double occupation induced byĤ RL are intuitive. Indeed, the tunneling termsâ b0 † Rσâ b0 Lσ (â b0 † Lσâ b0 Rσ ) create a double occupancy on the right (left) well of the b 0 -th band and shift N ↑L − N ↑R by two in the case that σ =↑. For instance, the tunneling termâ b0 † R↑â b0 L↑ (blue arrows in Fig. 9) transfers the spin-↑ particle of the state |Ψ; N ↑L , N ↑R from the left to the right well of the b 0 -th band resulting in the formation of a double occupancy on the right well of this band and modifying the occupation of spin-↑ particles to N ↑L − 1 for the left and N ↑R + 1 for the right well.
The approach followed to obtain the dominant perturbative correction to the eigenstates ofĤ L +Ĥ R , |Ψ; N ↑L , N ↑R , in the presence of the coupling termĤ RL is explicated below. First, we define the Hilbert space spanned by the degenerate eigenstates |Ψ; N ↑L , N ↑R as H 0 . Obviously since the states within H 0 are not directly coupled byĤ LR the first order perturbative correction to their energy vanishes. In order to obtain the first non-trivial correction to the energy of those degenerate states we have to treat the coupling termĤ LR within second order perturbation theory. Let us define the perturbative eigenstates up to second order in perturbation theory as |S i ≈ |S 1 the occupation of these states becomes highly suppressed and as a consequence these corrections can be neglected. Indeed, as Fig. 4 (b) reveals such corrections even beyond the effective tJU model contribute to a correction less than 2% to the fully correlated many-body eigenstates within the interaction regime B. Within the above mentioned approximation theĤ RL coupling term can then be substituted with the one of the effective Anderson exhange interaction which corresponds exactly to the form of the effective antiferromagnetic interaction appearing in Eq. (10) of the main text.

Appendix D: The Computational Method: ML-MCTDHF
To solve the MB Schrödinger equation i ∂ t −Ĥ |Ψ(t) = 0 we rely on the multilayer multiconfiguration timedependent Hartree method for atomic mixtures [26] (ML-MCTDHX). More specifically, a reduction of the ML-MCTDHX method for spin-1/2 fermions is employed which is referred to as the spinor-variant of the multiconfiguration timedependent Hartree method for Fermions (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][66][67][68][69][70]) and recently also applied in the field of ultracold atoms [26,[34][35][36][71][72][73]. MCTDHF is a variational method the key idea of which is to employ a time-dependent (TD) and variationally optimized MB basis set, which allows for the optimal truncation of the MB Hilbert space. The ansatz of the MCTDHF method can be summarized as follows. First, the MB wavefunction, |Ψ(t) is expanded on a TD number-state basis where A n (t) are the corresponding TD expansion coefficients. The TD number states | n(t) each possessing different occupation numbers n = (n 1 , . . . , n D ) read As Eq. (D2) reveals the time-dependence of this MB basis stems from the utilization of D different TD creation operators,â † j (t), j = 1, . . . , D. These operators create a fermion in the TD and variationally optimized single particle function (SPF) where the variational parameters φ jα (x; t) refer to the spatial distribution of the spin-α part of the j-th SPF andψ α (x) is the spin-α fermionic field operator. The operatorsâ j (t) satisfy the standard fermionic anti-commutation relations {â i (t),â † j (t)} = δ i,j and thus the MCTDHF ansatz takes explicitly into account the particle symmetry of the system. Note here that we have used the term spinor-variant when referring to our implementation of MCTDHF as each SPF, |φ j (t) , in our case is a general spinor wavefunction [see Eq. (D3)]. By employing the above mentioned ansatz Eq. (D1), (D2) and (D3) the time-evolution of the N -body wavefunction, |Ψ(t) under the effect of the HamiltonianĤ reduces to the determination of the coefficients A n (t) and the components of the SPFs, φ j↑ (x; t), φ j↓ (x; t). The latter in turn follow the variationally obtained MCTDHF equations of motion [26]. In the limiting case of D = N , the method reduces to the time-dependent Hartree-Fock approach neglecting all two-body and higher-order correlations. In the opposite limiting case of D = 2M p , where M p is the dimension of the basis for the SPF coefficients, MCTDHF is equivalent to a full configuration interaction approach (commonly referred to as "exact diagonalization" in the literature). The major advantage of the MCTDHF method when compared to methods employing a stationary single-particle basis is that the employed time-dependent basis is able to adapt to the correlation patterns emerging in the system during the dynamics and thus a smaller set of basis states is required for numerical convergence.
For our implementation we discretize the spatial coordinate by employing a harmonic oscillator discrete variable representation (DVR), which results after a unitary transformation of the commonly employed basis of harmonic oscillator eigenfunctions. 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 testify convergence, we increase the number of SPFs and DVR basis states such that the observables of interest (M , C F F ) do not change within a given level of accuracy which is in our case 10 −4 . More specifically, we have used M p = 80, D = 16 and M p = 80, D = 18 for the N = 4 and the N = 6 case respectively. Note that a full configuration interaction treatment of the above-mentioned systems in the employed primitive bases would require 2.63 × 10 7 number states for N = 4 and 2.12 × 10 10 ones for N = 6.