{'fraction_non_alphanumeric': 0.044727940038470894, 'fraction_numerical': 0.015322027902451966, 'mean_word_length': 4.387089090042878, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Many-body molecule formation at a domain wall in a one-dimensional strongly interacting ultracold Fermi gas

We analyze how the presence of the bound state on top of strong intercomponent contact repulsion affects the dynamics of a two-component ultracold Fermi gas confined in a one-dimensional harmonic trap. By performing full many-body numerical calculations, we retrieve dynamics of an initially phase separated state that has been utilized to excite the spin-dipole mode in experimental settings. We observe an appearance of pairing correlations at the domain wall, heralding the onset of a molecular faction at the interlayer between the components. We find that such a mechanism can be responsible for the stabilization of the phase separation.

Introduction-For decades, investigations of multicomponent mixtures-from alloys and polymers to biological systems and glasses-have provided a deep insight into details of intercomponent interplay. Specifically, an effective repulsion between constituents may result in their spatial separation [1,2] and induce an itinerant ferromagnetism in metals such as iron or nickel [3,4], where electrons spontaneously form extended, spin-polarized domains. This phenomenon can be understood within a simple mean-field framework proposed by Stoner, where a short-ranged screened Coulomb repulsion overcomes the Fermi pressure that favors a paramagnetic state [5]. While the simplified Stoner approach allows for qualitative description of many-electron systems, it does not capture effects related to beyond short-range interactions that may promote different, competing mechanisms suppressing ferromagnetism [6,7].
Since late 2000s experiments have tried to settle whether pairing processes prevent the ferromagnetic domains from appearing. In the initial attempts, some signatures, such as increase of a kinetic energy, suggested the onset of ferromagnetism, however these efforts proved inconclusive [29,31]. Only after the system was initialized in an artificial domain structure-in which components reside in their respective halves of a harmonic trap-the phase separation undoubtedly persisted for some finite time in a strong interaction regime [30,33]. Recently, time-resolved investigation of the competition between pairing and ferromagnetic instabilities provided deeper insight into many-body physics governing the dynamics in a quenched system [34].
On the theory side, approaches neglecting pairing processes showed that a ferromagnetic transition should take place in a three-dimensional geometry [10,12,14,15,[36][37][38][39][40][41]. However, inclusion of pairing has been done only at the mean-field level or by the introduction of phenomenological terms. In the context of Stoner ferromagnetism, an analysis of many-body eigenstates has been performed in few-body one-dimensional (1D) systems with contact interactions [42][43][44][45]. Nevertheless, in such a case the pairing process is supported only by purely attractive interactions in the absence of the repulsive core characterizing realistic atomic potentials.
To capture short-range details of the interatomic potential in 1D, we utilize the so-called three-delta potential [46][47][48][49][50][51][52][53]: where x denotes the relative position between fermions belonging to different spin components. The parameters c 0 > 0 and c < 0 describe contact repulsion and finite-range attraction at a distance , respectively. It can be understood as a first correction to the contact repulsion due to finite-range interactions and was shown to reproduce van der Waals forces under quasi-1D confinement [53]. Additionally, the intracomponent interactions are assumed to be negligible, as the system is brought close to the Feshbach resonance of the opposite spins and the temperature is very low. The paradigmatic scenario considered here involves a strongly repulsive core with a weakly attractive well (0 < −c c 0 ). Our aim is to thoroughly analyze dynamics of two initially separated Fermi clouds-a setup inspired by former experiments [33,54,55], see Fig. 1 (b). We argue that the potential W (x) can give insight to many-body processes involved in the competition between pairing and ferromagnetic instabilities.
Model -Mapping continuous, extended interactions onto discrete set of delta potentials was first analyzed in 1980s [46][47][48]. It was soon realized that the threedelta potential is a minimal extension of a realistic shortrange potential between atoms [49,50] and later utilized to study extended versions of the Gross-Pitaevskii and Lieb-Liniger models [51,52]. Recently, it was shown that the three-delta potential provides valuable insight into qualitative behavior of the many-body system that goes beyond standard contact interactions [53].
It should be stressed that in quasi-1D, the scattering theory implies a relationship c 0 + 2c = g 1D , where g 1D is a 1D mean-field coupling constant. However, in our analysis these parameters are kept free, since we consider full many-body model incorporating interplay between contact repulsion and beyond-contact attraction.
Note that such a model should also be attainable in experiments involving optical lattices, where on-site and nearest-neighbor interactions can be optically tuned.
For much smaller than other system length scales, it was shown that a simple contact interaction description can be retrieved with the effective coupling constant c eff = g 1D + corrections depending on , c 0 , c (for details see [52]). In our work, the length scale associated with is finite which affects the system beyond the contact approximation.
Our considerations are restricted to a balanced system consisting of M identical fermions with mass m in each of spin-↑ and spin-↓ components. At the beginning we analyze a two-body case (M = 1) determining a ground state wave function Ψ(x ↑ , x ↓ ) both for a free space and a harmonic potential confinement given by ω frequency, U (x) = 1 2 mω 2 x 2 , for details see Supplemental Materials (SM). While the point-like repulsion characterized by c 0 is assumed to be extremely strong, we modify the weak attraction strength c . Consequently, Ψ has a distinct cusp along x ↑ = x ↓ . However, an increase of |c | entails an enhancement of |Ψ| 2 around |x ↑ − x ↓ | ≈ , which is a premise of a bound state formation. Indeed, by energy considerations we found that for each c 0 there exists some critical value of c below which Ψ represents a bound state (see SM). As we show later, many-body dynamics of initially separated components with M > 1 reveals a similar probability density accumulation for finding two fermions with opposite spins at a distance , heralding the onset of the bound state.
For large M the many-body ground state can be studied within the lowest-order constrained variational (LOCV) approximation [36,[56][57][58][59][60][61], where the manybody wave function of the homogeneous, balanced twocomponent Fermi system is assumed to be represented by a product of two Slater determinants and a Jastrow factor, managing two-body correlations explicitly [62] (see SM). Within the LOCV approximation, the Jastrow fac-tor is varied to minimize the lowest order expanded total energy. Such an approach has proved to be very accurate in the description of three-dimensional quantum mixtures, rivalling quantum Monte Carlo's accuracy.
Taking the limit c 0 → ∞ we found that the interaction energy in the ground state within the LOCV approximation becomes negative for c ≤ c crit = − 2 m . Note that the critical value c crit does not depend on M , corresponding to c eff = 0 from [52] and agrees with a M = 1 exact solution (see SM). Therefore, when crossing c crit , we expect to observe a change of dynamical behavior of the initially separated ferromagnetic state, as signatures of bound structures may appear.
We now proceed to study dynamical properties of the initially phase separated system, where the fermionic components are confined to their respective halves of the harmonic trap by means of an extremely high potential barrier in the center. Both components do not have any meaningful overlap and are effectively noninteracting. The barrier is then instantaneously released and the system evolves freely with two components starting to push into each other. As the intercomponent overlap starts to accumulate, the corresponding correlations begin to appear around the trap center.
Two-body dynamics-First, we focus on M = 1 case which can be analyzed analytically (see SM). We stick to the case in which /d ≈ 0.238, where d = /mω is the harmonic trap length. For convenience, the following dimensionless interaction parameters, γ α = 2 1/4 cα d ω with α = 0, , crit, are utilized. Throughout our analysis, the repulsive core constant is set to γ 0 = 100 (γ 0 |γ |), corresponding to large g 1D , which guarantees a phase separation of fermionic clouds at the mean-field level.
Despite being very weak in comparison to the zerorange repulsion, the beyond-contact attraction has dramatic consequences on the system dynamics. When it is sufficiently weak, the stabilization of the initial phase separation is guaranteed due to the central core repulsion, cf. Fig. 1(d). Surprisingly, the stabilization is also present when the finite-range attraction is sufficiently strong, see Fig. 1(f). Such a behavior can be tracked back to a system interacting via only contact potential undergoing an interaction quench from an infinitely repulsive Tonks-Girardeau (TG) gas [63,64] to an excited super-Tonks-Girardeau (sTG) gaslike state of infinitely attractive atoms [65,66]. The initial state we consider can be decomposed into few lowest TG eigenstates due to the spatial separation between atoms, while γ < γ crit provides an effective strong attraction, despite large g 1D . From the exact solution one can see, that in the limit γ → −∞ the dynamics of the system approaches the dynamics under Hamiltonian with the same but with γ → ∞ (see SM). Further analysis of this mechanism will be provided in the subsequent work [67]. On the other hand, intermediate regime, γ ≈ γ crit , does not support stabilization of the spatial separation, as atoms mix with each other, cf. Fig. 1(e). Let us now proceed to a similar analysis in the few-body system. Many-body dynamics-In order to tackle the manybody problem, we study the Fermi-Hubbard model representing a discretized version of the continuous model discussed above. The contact interaction becomes the on-site one, while peripheral deltas correspond to the nearest-neighbor attraction. More specifically, we restrict ourselves to the lattice of 50 sites, where the lattice spacing ∆ equals the distance between deltas, ∆ = . To prepare the system in the initially separated state we employ a standard density matrix renormalization group (DMRG) method in the presence of the harmonic trap and a Gaussian barrier. It yields a matrix product state (MPS) in which the clouds of atoms belonging to different components are well separated and they are mirror images of each other with respect to the trap center, cf. Fig. 1(b). After removing the barrier at time t = 0 the system dynamics is investigated with the help of a recently developed algorithm combining one-site timedependent variational principle (TDVP) procedure [68] and a global basis expansion [69], cf. Fig. 1(c). For both DMRG and time evolution we use ITensor C++ library [70], where the codes for modified TDVP we employed are available thanks to the authors of [69] (for numerical details see SM).
In contrast to the M = 1 case, the domain structure remains stable for all the interaction strengths (see Fig. 1). It can be explained by analyzing the intercom-ponent density-density correlation function G(x, y; t) = Ψ(t)|n ↑ (x)n ↓ (y) |Ψ(t) /M 2 , withn σ (x) =ψ † σ (x)ψ σ (x), whereψ σ (x) is a canonical Fermi field operator for spin σ at x and |Ψ(t) is the time-evolved MPS at time t. The correlation function is normalized such that G(x, y; t)dxdy = 1. In Fig. 2 we compare G (x ↑ , x ↓ ) obtained for M = 1, 3 and different finite-range attraction strength γ , where ... denotes a temporal average over 5 trap periods. Far from γ crit both systems evolve similarly, suggesting that stabilization mechanisms described in the previous section are present also for M > 1.
On the other hand, in the intermediate regime γ ≈ γ crit , an escalation of G at |x ↑ − x ↓ | ≈ is revealed for M = 3. We interpret this enhancement as the onset of the bound state contribution (molecule faction) appearing in the course of time evolution. It can be viewed as an extension of similar structures present in the relative wave function of the bound state in the two-body problem (see SM). It is also supported by the fact that this growth takes place at the critical interaction strength predicted by the LOCV approximation and M = 1 exact solution. Note that in the M = 1 case, the dynamics reveals no molecular formation as there are no additional atoms to absorb excess kinetic energy.
This effect can be further investigated by analysis of the pair distribution function g(r; describes spatial correlations between spin-↑ and spin-↓ fermions at a distance r. The value of g(r; t) corresponds to the probability density of finding two fermions of opposite spins at a distance r. A temporal variability of g(r; t) obtained for M = 1, 3 and different γ is presented in Fig. 2(k)-(o). For M = 1, independently of the γ value, g(r; t) oscillates, closely following the dynamics of single-particle densities, cf. Fig. 1(e)-(f) and upper panels of Fig. 2(k)-(o). In contrast, for γ ≈ γ crit , the M = 3 case reveals a steady and gradual growth of g(r; t) at r = , being the distance at which signatures of anticipated bound pairs are expected. Since the dominant contribution to g( ; t) comes from the trap center, cf. panels (h) and (i) of Fig. 2, we interpret this result as a footprint of molecular faction in the interlayer between components. Far from γ ≈ γ crit , the resulting g( ; t) is inappreciable, suggesting that stabilization mechanisms noticed and described in the two-body case are also present in larger systems.
To study the molecular faction in more detail, one can compute N ↑-↓ = M g( ; t) representing the expectation value of the number of ↑-↓ pairs of size ≈ . Since the components in the initial state are spatially separated, the anticipated ↑-↓ bound structures may appear in the course of time evolution. Therefore, in Fig. 3(a) we investigate how the temporal average N ↑-↓ depends on M and the beyond-contact attraction γ . It is striking that while N ↑-↓ turns out to be almost M -independent for FIG. 2. Temporally averaged over 5 trap periods densitydensity correlation functions G (x ↑ , x ↓ ) calculated for M = 1, panels (a)-(e), and M = 3, panels (f)-(j), and for different finite-range attraction strengths γ . In accordance with the results shown in Fig. 1, mixing between the two Fermi clouds takes place only for M = 1 in the presence of intermediate attraction, see (b)-(d). In contrast to the M = 1 case, for M = 3 with moderate attraction one can observe an accumulation of G (x ↑ , x ↓ ) in the trap center, see (g)-(i). In right panels (k)-(o) we show a temporal behavior of g(r; t) obtained for M = 1, 3 and different γ . While for M = 1 the pair distribution g(r; t) oscillates together with single-particle densities plotted in Fig. 1(d)-(f), for γ ≈ γcrit the M = 3 case reveals a steady and gradual increase of correlations between opposite-spin fermions separated by r ≈ . Such correlations are strictly related to the presence of the finite-range attraction and are responsible for a domain stabilization. γ > γ crit , at γ ≈ γ crit one observes an abrupt splitting between the results obtained for M = 1 and M > 1. For both M = 2 and M = 3, N ↑-↓ is very similar and reveals a dramatic growth up to ≈ 0.2. When γ becomes more negative, it rapidly decays approaching the values obtained for M = 1. The escalation of N ↑-↓ for M > 1 around γ crit can be better understood when looking at the spatial distribution K (x). That is, as shown for  Fig. 3(b), in the vicinity of γ crit there is a significant accumulation of ↑-↓ pairs of size in the center of the trap. This is an additional signature of an appearance of a molecular faction that resides between the two components forming a domain wall and thus providing a phase separation mechanism. It is reminiscent of the microemulsion of two fermionic components and dimers investigated in Refs. [34,71], where in the threedimensional trap bound pairs were also present at the interlayer separating components.

M = 3 in
Short-range atomic potentials-Despite the fact that the model Hamiltonian we consider could be analyzed experimentally through geometries involving optical lattices, the paradigmatic case of the Stoner instability in a two-component Fermi gas not constrained to the additional lattice is of particular interest. In quasi-1D settings the perpendicular confinement may be very strong and confining length scales become comparable to the range of beyond-contact interactions making effects of finiterange coupling vastly pronounced [53].
In such settings, finite-range corrections to the usual contact interaction can be modelled by the pseudopotential W p (r) = g 1D (1 + g p 2 )δ(r), where p = −id/dr is the 1D momentum operator and g 1D and g can be determined in a realistic scenario involving the Feshbach resonance and strong perpendicular confinement (see Ref. [53] and SM for details). When p is discretized, the three-delta potential (1) is recovered, however there is some ambiguity-there are multiple choices of c 0 , c and that correspond to the same parameters g 1D and g describing identical scattering properties (see SM for a detailed calculation). To map the realistic potential onto the three-delta one unambiguously, the energy of the true bound state has to be also reproduced.
However, in our work we analyze short-range correla-tions emerging at a domain wall of the phase-separated state where the strength of the beyond-contact attraction is treated as a free parameter. Intuitively, length scale should reflect a range of the real interatomic potential, but in our analysis it is much larger. Thus, the scaling with decreasing needs to be addressed. Many-body simulations we performed are not attainable for significantly larger number of sites that would allow for studies of the system behavior when is smaller. Nevertheless, we have found no evidence of qualitative differences when was increased twice and other parameters were adequately rescaled. This may suggest that our findings hold also if smaller-and more realistic in unconstrained atomic systems-values of were to be used.
Conclusions and outlook -We have proposed to utilize the three-delta interaction potential to study competition between many-body effects, such as pairing and ferromagnetism, in an ultracold two-component Fermi gas.
We have analyzed the dynamical stability of the experimentally inspired artificial domain structure and have found three stabilization mechanisms. First one is due to usual mean-field repulsion of the components, while the other two involve many-body processes. The process dominating in the presence of strong attraction can be associated with the stability of the super-Tonks-Girardeau gas, while in the regime of moderate attraction, the stabilization is due to the molecule formation at the interlayer between the components. As a future line of work, the three-delta potential can become a powerful asset in the many-body studies of systems where repulsive and attractive interactions compete to form novel quantum phases. In this section we investigate a one-dimensional problem of two particles interacting via a three-delta potential and determine a criterion for a bound state formation in the absence of external confinement. The Hamiltonian of our system reads and we assume that c 0 > 0 and c < 0. We decouple center of mass (CM) and relative motion by introducing new We focus on H rel and seek for localized eigenstates (bound states) describing a pair bound due to attractive part of W (r). In general the potential may have at most two bound states -one symmetric (+) and one antisymmetric (−). The problem is solved by imposing the following standard continuity and δ-potential-related conditions where k = 2m|E|/ , with E being the binding energy. The abovementioned conditions lead to the fol- where the dimensionless parameters u 0 = mc 0 / 2 , u = mc / 2 and t = k were introduced. The energies E can be determined thanks to relations . As the lowest energy state is always symmetric, our potential hosts a bound state only if (3) has a positive solution. Assuming a very strong contact repulsion c 0 → ∞, we take the limit t → 0 + in (3) and find the critical value of attractive coupling u crit = −1 below which there exists a bound state in the system. While the regime of bound state existence is shown in Fig. 1(a), exemplary shapes of the corresponding bound state wavefunctions can be found in Fig. 1(b).

II. 1 + 1 PROBLEM IN A HARMONIC TRAP
Let us now study the problem of two initially separated particles interacting via three-delta potential. We show that exact solutions provides us with a lot of insights about the time evolution and the structure of the energy spectrum. Additionally, it will serve as a benchmark for numerical methods employed when investigating larger systems.

A. Solutions
Here, applying methods similar to [1,2], we solve the problem that may be viewed as an extension of the problem analyzed in [1]. Namely, we consider the Hamiltonian As in the previous section we introduce R and r variables obtaining H = H CM + H rel , where H CM describes the standard Hamiltonian of harmonic oscillator and thus we focus only on H rel = − 2 2m ∂ 2 r + mω 2 2 r 2 + W (r). In the absence of potential W (r) the corresponding Schrödinger equation describes a standard harmonic oscillator in 1D, i.e., Eψ(r) = − 2 2m ∂ 2 r ψ(r)+ mω 2 2 r 2 ψ(r). The most general solutions of the problem can be grouped into symmetric and antisymmetric ones and cast into the following forms where d = /mω is the oscillator length, and the energy-related coefficients κ = 3/2 − E/ ω and σ = 1/2−E/ ω [2,3]. Additionally, for convenience we introduced functionsF ξ,ζ (r) = exp −r 2 /2d 2 1 F 1 ξ; ζ; r 2 /d 2 andŨ ξ,ζ (r) = exp −r 2 /2d 2 U ξ, ζ, r 2 /d 2 , where U and 1 F 1 are confluent hypergeometric functions of the first and second type, respectively. It is important for further analysis to list the most important properties of 1 F 1 and U in order to construct physical eigenstates (taking delta functions into the account) similarly to the case without confinement [2,3]. If σ = −n with n = 0, 1, 2, . . ., then 1 F 1 − n 2 ; 1 2 ; x 2 ∼ H 2n (x) and for κ = −n (n = 0, 1, 2, . . .) we get 1 F 1 − n 2 ; 3 2 ; x 2 ∼ H 2n+1 (x). In this way we retrieve free oscillator solutions out of general forms (5) and (6). However, when interactions shift values of energies so that σ and κ are no longer nonpositive integer numbers, the functionsF σ which are continuous at r = 0 but reveal discontinuities in derivatives at this point, which we use to fulfill the abovementioned Dirac delta conditions. These properties motivate us to propose solutions in the form Now, we proceed similarly as in the case with no confinement, i.e., we impose continuity and δ-potentialrelated conditions at r = 0, ± . Let us start with the symmetric solutions, where continuity of the wave function at r = ± implies and δ-induced condition at r = 0 yields (see [3] for details) Applying the same conditions but at r = ± together with relations (9) and (10) one finds where Ξ σ = √ π ωdσ/c 0 Γ 1 + σ 2 + 1/Γ 1 2 + σ 2 was introduced for convenience. By plugging these equalities into the continuity condition, we determine the equation for σ (directly related to energies) When dealing with antisymmetric solutions, the corresponding wave functions are not affected by δ-potential at r = 0 and thus it is enough to focus on r = ± cases, where from the continuity and δ-related conditions we find Finally, κ gives the energies of antisymmetric states and can be determined from Our eigenproblem is now solved. Energies are readily obtained as numerical solutions of (13) and (16). Wave functions are determined from Eqs. (11), (12) and (14) together with the normalization condition. Ground states for different parameters are presented in Fig.2. In summary, let us note that full eigenstates of the two body problem can be cast into the following form where φ n , n = 0, 1, 2, . . . is a harmonic oscillator eigenstate whereas ψ µ represents an eigenstate of H rel . Here µ denotes either κ or σ depending on the symmetry of a given state in the relative coordinate r.

B. Super Tonks-Girardeau-like behavior
As discussed in the main text, the system dynamics in the presence of strong finite-range attraction turns out to be very similar to the one observed in the case of a purely repulsive contact interaction. We argue that this counterintuitive behavior of the system can be related to the phenomenon of the super Tonks-Girardeau gas (sTG) [4][5][6][7].
The sTG is a highly excited, metastable state of matter where the particles do not form bound states despite a strong interparticle attraction. Such a state has a large overlap with the ground state of infinitely repulsive Tonks-Girardeau gas (TG). This phenomenon was understood over ten years ago in the Lieb-Liniger model [5]. It turns out that a sudden quench of interactions from strongly repulsive to strongly attractive does not drive the system out of the equilibrium. A trace of such a behavior can be observed already for two confined particles interacting via contact potential [1,6]. As found in Ref. [1], the spectrum of the relative Hamiltonian in the strongly attractive regime equals the spectrum of the strongly repulsive one with an additional very deep bound state.
In this section we show that the situation is similar in the case of the three-delta interaction. Let us consider the limit of c 0 / ωd → ∞ in Eqs. (13) and (16), where the spectrum becomes double-degenerate and all the states can be grouped into pairs consisting of one symmetric and antisymmetric characterized by equal energies. The resulting eigenenergies of the relative motion as a function of c for fixed are presented in Fig. 3. We find that when c / ωd → −∞ the spectrum equals the spectrum of the system in the limit c / ωd → ∞ with additional two infinitely deep bound states (symmetric and antisymmetric). In general, the dynamics of the system is given by where E n denotes the n-th harmonic oscillator eigenenergy, E µ is the µ-th eigenenergy of H rel and C n,µ = Ψ n,µ Ψ 0 , where Ψ 0 and Ψ n,µ represent the initial state and the eigenstate of the form (17), respectively. Let us focus on the case c / ωd → −∞. As the initial state, Ψ 0 , corresponds to two initially separated atoms, the wave function is zero when x = y and is still very small for |x − y| = . On the other hand, the bound states, Ψ b , are tightly localized at |x − y| = and almost zero elsewhere. Consequently, Ψ 0 Ψ b → 0 and both symmetric and antisymmetric bound states effectively do not contribute to the dynamics (18). The remaining states that do have a meaningful overlaps, are at the same time eigenstates of the system with c / ωd → ∞. As a result, the dynamics of the system with strong finite-range attraction closely follows the dynamics of the system with strongly repulsive potential c / ωd → ∞. We will elaborate on the mechanisms presented here in subsequent work [8].

III. LOCV APPROXIMATION
In this section we go beyond the two-body problem and consider a uniform, one-dimensional Fermi gas of N ↑ spin-↑ and N ↓ spin-↓ atoms characterized by mass m. The balanced gas, N ↑ = N ↓ = M = N/2, is described by the following Hamiltonian where x σ i denote spatial coordinates of spin-σ atoms and W (r) is a short-range two-body interaction potential between particles belonging to different spin components, with no spin-flipping allowed. Here, the total wave function is approximated by means of the Jastrow-Slater Ansatz [9]: where D ↑ and D ↓ are Slater determinants: of plane waves ϕ σ i , and J is a Jastrow factor with f being a Jastrow function that accounts for a twobody relative wave function of opposite spin atoms. Such an Ansatz is not valid for general Fermi systems, e.g. neutron matter, however it proved to be a good approximation for cold atomic systems in which short-range interactions are dominant [10]. We assume that the Jastrow factor slightly modifies the whole wave function at short relative distances only and that the parameter Let us proceed to evaluate the normalization of the wave function. By the orthonormality of plane waves, we find that up to the first order in F ij which under the transformation r = x − y, R = x + y can be cast into n dr|f (r)| 2 = N , where n σ = n/2 denotes the σ-particle density. Now, we need to introduce physical assumptions on the behavior of the Jastrow function. We demand that for long separations between the atoms, J asymptotically tends to 1 not to affect the long range behavior of the wave function. Even more than that, we introduce the healing length d that is of order of average interparticle separation 1/n, beyond which the Jastrow function does not modify Ψ. It expresses the intuition that on average only the nearest neighbors of opposite spins are correlated and the correlations with atoms at longer distances are negligible. This constraint and the requirement for smoothness at the healing distance, are two out of three constraining conditions. The last one comes from the normalization of the Jastrow factor. We assume that f deviates from 1 only inside the sphere with radius d, K(d), and on average there is only a single ↑-↓ pair inside K(d). Therefore, starting with n dr|f (r)| 2 = N , assuming L to be the system size one rewrites which results in It gives the last constraint These three constraints are needed in variational calculations in order to reproduce the experimental data [10,11]. Unconstrained calculations were utilized to study nuclear matter, however each time Jastrow function tended to be unphysically long-ranged. On the other hand, constrained studies provided extremely good fit to the experimental data [12,13]. Note that introduction of these constraints will also imply some kind of the effective two-body interaction. If the Jastrow function were to satisfy the usual two-body Schrödinger equation for the relative wave function, it would not tend to 1 for long distances. Therefore, some kind of alternation of interaction is needed in order to satisfy the constraints, which will be addressed later.
We will now proceed to evaluate the energy within the Slater-Jastrow Ansatz. The noninteracting contribution comes from the kinetic energy of the orbitals and in large atom number limit is equal to E 0 /V = nE F /3, where the Fermi energy E F = 2 π 2 n 2 V /8m. The interaction part comes as a two-body cluster expansion of the total energy that takes into account both the kinetic energy from the Jastrow function and the influence of a bare two-body potential. Namely, Now, Eq. (29) needs to be extremized with respect to variations of f . However, in order for f to satisfy the constraints, two-body interaction needs to be altered. The easiest way to do so is by adding additional constant external potential λ in such a way that W → W − λ [12]. It is not the only way of renormalization, however it was shown to provide reliable results and to be equivalent to so called Moszkowski-Scott separation [14] and Brueckner theory [15]. This procedure can be understood as an effect of average pressure of further neighbors and sometimes it is interpreted as a contribution of two-body correlations to average field felt by a given atom.
Therefore, to find a minimum, one needs to solve It yields two-body Schrödinger-like equation that needs to be solved with constraints The interaction energy is then The final result, E int /V = nλ/2, indicates that on average we have n/2 pairs, each one contributing energy λ.
The LOCV approximation proved to provide reliable experimental predictions in many settings, from dense nuclear matter to ultracold quantum mixtures. However, it was mostly utilized in three-dimensional geometries, in which the level of its accuracy matches state-of-theart quantum Monte Carlo calculations [10,13]. Analysis of low dimensional systems was scarce, but the preliminary results in two dimensional Fermi gas suggest that LOCV fares worse in such geometries. The probable reason is a larger role of quantum fluctuations in lower dimensions and resulting need to go beyond the lowest order in the cluster expansion to match the more sensitive approaches. Nevertheless, the qualitative match can be expected as LOCV should reproduce leading order behavior of systems considered.

A. Three-delta potential in LOCV approximation
We now consider a specific potential consisting of three delta contributions where we choose the contact interaction strength c 0 to be large and positive to mimic the repulsive core of realistic interparticle potential. On the other hand, beyondcontact part, c , is chosen to be small and negative to support a weakly bound state and provide the simplest extension of a contact part towards full van der Waals forces. In a three-dimensional two-component Fermi mixture, the contact repulsion given by the Fermi-Huang pseudopotential ∝ δ(r) ∂ ∂r r reproduces both the weakly bound pairs (attractive branch of the many-body system) and the repulsive, anticorrelated, phase-separated Fermi clouds (repulsive branch). This is not the case in a one-dimensional space as repulsive contact interactions do not support a bound state.
To simplify our considerations, we will work in the regime of infinite repulsion, c 0 → ∞, yielding f (0) = 0 We compute the interaction energy within the LOCV approximation solving Eqs. (31)- (34). If the bound state is considered, λ = − 2 κ 2 /2m is negative. Then, we assume the Jastrow function to take form and needs to be self-consistently solved for A, B, D, d, κ with constraints (33)- (34). It simplifies to the following two equations where we define κ 1 = κ , κ 2 = κd, λ 2 = mc / 2 . The interaction energy density then equals It is evaluated numerically and presented in Fig. 4(left). Note that the bound state is not supported for every value of c . Most notably, numerical evaluation allows to write down the simple expression for the existence of the bound state which does not depend on the gas density. On the other hand, similar procedure for the first excited state (the repulsive branch) yields where λ = 2 k 2 2m is positive and we define k 1 = k , k 2 = kd, λ 2 = mc / 2 . The result is presented in Fig. 4(right).
The lowest energy branch of positive solutions to the constrained equations tends to 0 at mc / 2 ≈ −1, i.e. the interaction strength at which the bound state starts to be supported. In Fig. 5 the energy dependence on parameter c for the many-body ground state at a given density is presented. The sharp transition is clearly visible, suggesting nontrivial behavior of the system at the critical value of c .

IV. SHORT-RANGE ATOMIC POTENTIALS
A. Scattering in quasi-one-dimensional geometry Here we sum up the analysis of Ref. [16] where the authors show how finite-range corrections to the standard pointlike delta interaction can be modelled through the effective pseudopotential where g 1D denotes a one-dimensional mean field interaction strength, g is a parameter depending on the details of the physical situation, and following Ref. [16], p represents a one-dimensional momentum operator.
The system under consideration involves ultracold atoms interacting via the attractive part of the van der Waals potential W vdW (r) = −C 6 /r 6 with a characteristic length (mean scattering length)ā = 2π Γ(1/4) 2 (2µC 6 / 2 ) 1/4 , where µ is the reduced mass of the atomic pair [17] and Γ is the Euler gamma function. It can be rephrased in terms of the van der Waals radius R vdW asā = 0.955978R vdW .
In the presence of tight harmonic confinement 1 2 µω 2 ρ 2 ρ = y 2 + z 2 and ω denotes trapping frequency, the usual Fermi-Huang pseudopotential can be generalized to account for an energy dependence where 2 k 2 /2µ is the kinetic energy of the relative motion and δ 3D (k) is the the phase shift due to interactions. Consequently, the corresponding energy-dependent scattering length reads a 3D (k) = − tan δ 3D (k)/k. One notes that the full three-dimensional scattering analysis can be reduced to the 1D problem in the s-wave channel with the pseudopotential [16,[18][19][20][21] where Here, δ 1D (p) is a 1D phase shift depending on 1D momentum p, where k 2 = p 2 + 2/d 2 ⊥ with d ⊥ = /µω. Expanding (47) in small p up to the second order one finds with whereC = ζ(3/2)/8 ≈ 0.3265 and C = −ζ(1/2) ≈ 1.46035. The effective range r 3D can be determined from the expansion In vicinity of the Feshbach resonance, the effective range of a single-channel van der-Waals potential can be determined analytically [22][23][24], where R * = 2 /(2µa bg ∆δµ), a bg is the background scattering length away from the resonance, ∆ denotes the resonance width and δµ is the magnetic moment difference between channels. This expression approximately holds for open channels, but can fail in the presence of the closed ones [24]. Then, it is more accurate to use formula where the magnetic field dependence of a 3D reads with B 0 being the resonance point, and the fitting parameters v, r 0 and a ex depend on a particular resonance and can be evaluated through the numerical solutions of a full multi-channel Schrödinger equation.
B. Mapping onto three-delta potential We will now proceed to discretize effective interaction potential (44). The starting point is to consider two-body problem with a total wave function where R = x + y describes the decoupled center-of-mass motion of two atoms, r = x − y describes the relative motion, and x, y represent positions of atoms. For the delta interaction W a (x−y) = gδ(x−y−a), the interaction energy reads (56) The position representation of potential (44) can be written as After discretizing the second derivative, i.e., the energy associated with the last term of the sum in (57) read After realizing that one finds the interaction energy related to the finite-range correction By Eq. (56), it corresponds to the following potential where ∆x can be arbitrarily chosen. It is convenient to denote ∆ = |g |/α with some arbitrary real number α > 0. In result we reproduce the three-delta potential with and inversely g 1D = c 0 + 2c , g = −4 2 c /(c 0 + 2c ).
Note that g and g 1D can be uniquely determined by c 0 , c and , while α provides an additional information how to perform inverse transformation. It is due to the fact that g and g 1D describes only the scattering properties of the physical system. Effectively, α fine tunes the length scale introduced through the equality = |g |/2α and can be evaluated by considering other properties of the two-body problem, such as the bound state energy.

V. DISCRETIZATION OF THE CONTINUOUS MODEL
The Hamiltonian of the considered system readŝ whereψ σ (x) denotes the canonical Fermi field annihilation operator satisfying standard fermionic anticommutation relations. We discretize the space into sites with spacing ∆ and introduce new fermionic operators acting directly on lattice sites, i.e.,ψ σ (x i ) → √ ∆â σ,i where x i represents a position of i-th lattice site.
In result we arrive at the Fermi-Hubbard model with beyond on-site interactionŝ where T = 2 /2m∆ 2 , U = mω 2 /2,n σ,i =â † σ,iâ σ,i is the spin-σ particle number operator at site i, and s represents the distance in lattice sites, i.e., = s∆.

VI. DETAILS OF NUMERICAL SIMULATIONS
To study the many-body dynamics of the considered problem we employ matrix product states (MPS) [25,26] ansatz with open boundary conditions where the system of size L = 10 has edges located at x OBC = ±L/2. The N s sites, representing the lattice in a discretized space, are located at positions x j=1,...,Ns = −L/2 + ∆(j − 1/2) with ∆ = L/N s . The harmonic confinement is realized by employing the trap of frequency ω = √ 2 centered at x = 0, which guarantees that particles in the systems under considerations (M ≤ 3) never feel the hardwall boundaries. The initially separated state is prepared with the help of a standard density matrix renormalization group (DMRG) method [25][26][27][28][29] in the presence of a Gaussian barrier V b (x) = Ae −x 2 /2σ 2 , starting with a state where fermions of opposite spins occupy sites in different halves of the trap far from x = 0. For this purpose we used A = 10(M + 1) and σ = 0.2.
Generally, time evolution within the MPS approach is error-prone and computationally demanding. It is due to the fact that together with a ballistic increase of an entanglement entropy in the course of dynamics, it is required to deal with increasing bond dimension. This issue was partially overcome by a recent development of the time-dependent variational principle (TDVP) algorithm [30][31][32], which allows for a time propagation with a fixed bond dimension (one-site TDVP). It is much less error-prone in comparison with other earlier techniques like the time evolving block decimation (TEBD) routine [33]. Nevertheless, to reproduce a system dynamics properly, the bond dimension has to be first enlarged enough, which can be done for example as in the "hybrid" TDVP strategy, where the desired bond dimension is obtained through the initial application of the twosite TDVP scheme [33][34][35]. Here, we employ a recent approach in which the basis for MPS Ψ is enlarged by means of a subspace expansion by global Krylov vectors { Ψ ,Ĥ Ψ , . . . ,Ĥ k−1 Ψ }, whereĤ is the Hamiltonian governing the system dynamics and k denotes the so-called Krylov order. The method is meticulously described in Ref. [36], where the authors argue that this strategy is more accurate and more efficient than the twosite TDVP approach. Moreover, this method turns out to be more reliable in the case of beyond-nearest-neighbor interactions. We use the algorithm implementation based on the ITensor library provided by the authors of [36] under the ITensor/TDVP repository. The calculations we performed with the time step ∆t = 0.02 and k = 3, where the truncation error of each application ofĤ to Ψ and the truncation error controlling diagonalization of the sum of the reduced density matrices were chosen to be equal 10 −8 . Additionally, we preserve exact unitarity by imposing no truncation during one-site TDVP sweeps. This choice of parameters guarantees a balance between cost and accuracy, when testing in the most demanding M = 3 case. We let the bond dimension χ of the time-evolved MPS grow up to pre-determined value χ max = 750 and switch off the global basis expansion routine when χ max is reached.
First, we benchmark obtained results in the M = 1 case, where the exact results in the continuous space are known. Namely, we consider the system with different numbers of sites N s = 50, 100, 150, but with the same = 0.2 in the three-delta W (x) potential. Since ∆ = L/N s , in order to compare the numerical results with the analytical predictions, it is required to take into account different beyond-contact attraction term in the discretized space, i.e. different values of the integer s in the Hamiltonian (65). That is, while for N s = 50 we employ nearest-neighbor attraction (s = 1), for N s = 100 and N s = 150 we incorporate next-tonearest-neighbour (s = 2) and second-next-to-nearestneighbour (s = 3) interactions, respectively. The comparison between dynamics of the center-of-mass position of ↑-component, X CM,↑ , determined analytically in the continuous model and numerically in the discretized system for γ = −3, −5, −7, −10 is illustrated in Fig. 6. Note that the results obtained for N s = 100 and N s = 150 are almost identical to the analytical ones, showing the convergence and reliability of the TDVP evolution we employed. On the other hand, in the case of N s = 50 the considered center-of-mass position evolves slightly different in comparison with the analytical results. Nevertheless, the discreteized system behavior is still very similar to the continuous one, which allows us to suppose that conclusions we formulated basing on the analysis performed in the discretized model with N s = 50 are also valid in the continuous systems.
Secondly, we show the convergence of our main results determined for N s = 50 by comparing g( ; t) obtained with different maximal bond dimension χ max = 250, 500, 750 for M = 2 and M = 3 and different beyondcontact attraction strengths γ = −3, −5, −7, −10 (see Fig. 7). In addition, in Fig. 8 we present the comparison between dynamics of the von Neumann entanglement entropies, S j , calculated across the bond between j and j + 1 sites for M = 3, different χ max , and two central bonds j = 20, 25.