Ground states and dynamics of population-imbalanced Fermi condensates in one dimension

By using the numerically exact density-matrix renormalization group (DMRG) approach, we investigate the ground states of harmonically trapped one-dimensional (1D) fermions with population imbalance and find that the Larkin-Ovchinnikov (LO) state, which is a condensed state of fermion pairs with nonzero center-of-mass momentum, is realized for a wide range of parameters. The phase diagram comprising the two phases of i) an LO state at the trap center and a balanced condensate at the periphery and ii) an LO state at the trap center and a pure majority component at the periphery, is obtained. The reduced two-body density matrix indicates that most of the minority atoms contribute to the LO-type quasi-condensate. With the time-dependent DMRG, we also investigate the real-time dynamics of a system of 1D fermions in response to a spin-flip excitation.


Introduction
One-dimensional atomic gases have been vigorously investigated over the last decade. By using a strongly anisotopic Ioffe-Pritchard-type magnetic trap [1] or pairs of counterpropagating laser beams [2], elongated traps of atoms have been realized. Tens to hundreds of ultracold atoms have been trapped in such an elongated potential that the maximum kinetic energy of the atoms are much smaller than the energy separation between the two lowest transverse levels of the potential [3,4,5,6,7]. Thus the physics of the system becomes essentially one-dimensional [8,9]. For the case of bosons, a strongly interacting Tonks-Girardeau regime was realized [3]; no equilibration of an integrable system -"quantum Newton's cradle" -was demonstrated [5]; and the crossover from thermal to quantum noise was also observed [6].
For the case of fermions, the 1D atomic gas allows one to study an analogue of a system of electrons in solids. The numbers of atoms in trappable hyperfine states can be controlled individually, so that by populating two hyperfine states simultaneously, one can study effects of population imbalance on the properties of one-dimensional systems. While a true longrange order is forbidden in 1D even at T = 0 by the Hohenberg-Mermin-Wagner theorem, (quasi) condensation and superfluidity are possible in a finite-sized, trapped system. The search for Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) [10,11] pairing states has been of major experimental interest in 1D fermionic systems.
The confinement-induced two-particle bound state was observed [4] in an array of quasi one-dimensional tubes with equal populations of two hyperfine states, validating the 1D treatment of these tubes. More recently, the Rice group [7] has realized quasi one-dimensional systems of population-imbalanced cold atoms and measured the density profile for each of the hyperfine species.
In this article we first report our study on the ground state of a system of populationimbalanced fermions in harmonic traps [12]. The system has also been numerically studied by using the density-matrix renormalization group (DMRG) [13,14,15,16,17,18,19], Bogoliubov-de Gennes (BdG) equations [20,21], and Quantum Monte Carlo (QMC) method [22,23,24]. The phase diagram including an FFLO superconductor region for a 1D homogeneous system was obtained by bosonization [25]. Studies on trapped systems based on the exact solution of the Gaudin-Yang model were carried out in Refs. [26,27,28,29,30]. Other related studies include the cases with unequal (effective) masses and the balanced or imbalanced numbers of atoms [31,32,33,34]; an angular FFLO state in a toroidal trap [35]; effects of particle-correlated hopping in an optical lattice [36]; imbalanced Fermi gases on two-leg ladders [37] and two-dimensional optical lattices [38]; composite particles of more than two particles [39]; universality and Efimov states in systems where two or more species of atoms are confined in traps having different dimensions [40,41]; transport through a Yjunction and ring-geometry systems [42].
We note that, for the ground state that is symmetric under time reversal, the condensate cannot be of a Fulde-Ferrell type, because the Fulde-Ferrell state is described by a pair correlation function that is a complex-valued function of the spatial distance. On the other hand, the Larkin-Ovchinnikov state is not forbidden by the time-reversal symmetry because the correlation functions can be real-valued.
Here we numerically study the ground-state wavefunction of the trapped system by discretising the system and the DMRG [43,44,45] method. We calculate the density distribution of individual spin components in the trap, which is a physical quantity measured in in-situ light absorption or phase shift measurements [46,47,48], and show that the densitydifference distribution exhibits oscillations which are the hallmark of the Larkin-Ovchinnikov state [49].
In the second part of the paper, we study the real-time evolution of a system of trapped Fermi atoms after a sudden spin flip. Because of the large interparticle distance and low speed of the atoms, the dynamics is extremely slow as compared to electrons in solid state physics. Moreover, the atoms are trapped in an optical and/or magneto-optical potential in the vacuum, so that the dynamics can be optically observed in real time. As the state of the atoms can be manipulated by microwaves, it is of interest to study the response of the 1D ground state to a change in population.

The ground state
First, we study the ground state of the trapped, population imbalanced fermionic system in 1D within the single-channel model. The effect of a molecular state in the BCS-BEC crossover in 1D systems [50] has also been studied in population imbalanced systems [18,51], but here we focus on the case in which the contribution of the molecular state can be neglected.
The interaction between ultracold atoms is described by the s-wave scattering which can be approximated by the Lee-Huang-Yang pseudopotential [52] with the bare 3D coupling constant g 3D s . The relation between g 3D s and the s-wave scattering length a 3D s is given by g 3D s = 2π a 3D s /µ, where µ is the reduced mass.
The relation between a 3D s and the scattering length in a 1D trap a 1D s has been given in [9], for the case in which the transverse confinement can be approximated by a 2D harmonic potential with (angular) frequency ω ⊥ . Provided that the size of the ground state of the transverse confinement a ⊥ ≡ ( /µω ⊥ ) 1/2 is much larger than an effective range r 0 of the bare potential, a 1D s is given by with C = 1.4603 . . .. The 1D scattering amplitude for the relative wavenumber k z is expressed as which, in the low-energy limit (|k z a ⊥ | ≪ 1), reduces to the scattering amplitude derived from the one-dimensional δ potential with g 1D s = − 2 /(µa 1D s ) [9]. When we have two fermionic species with the equal mass m 0 , the reduced mass is µ ≡ m 0 /2. The two hyperfine states can be regarded as a pseudo spin-1/2 system with spin up (| ↑ ) and down (| ↓ ) states. We assume that the up-spin state is the majority state. We start with an effective 1D Hamiltonian for the continuum, where m is the atomic mass, V(z) ≡ κz 2 /2 is the harmonic potential in the z direction and g is the coupling constant. We discretize this Hamiltonian to perform calculations in the DMRG framework.
We take a characteristic size of the potential, l, as the unit of length, and denote a characteristic trap depth by A: V(±l) = A. We discretize the system by introducing a lattice of L sites with lattice constant d ≡ 2l/L. The transfer amplitude between the neighboring sites is chosen to be J ≡ 2 /(2md 2 ). The band dispersion, where k = p/ is the wave number of the atom with momentum p, reproduces the energy dispersion of the atom in the free space, E(k) = p 2 /(2m) = 2 k 2 /(2m), in the vanishing limit of the filling factor (L → ∞). The interaction g 1D δ(z a,↑ − z a ′ ,↓ ) is approximated by an on-site interaction, U L−1 i=0n i,↑ni,↓ , where i runs over all lattice sites 0, 1, . . . , (L − 1), with coupling constant U ≡ g 1D /d. In the following we take = m = 1. We consider the case of negative U, which corresponds to a positive a 1D s and attractive interaction between atoms with opposite spins.

Relations between the experimental and discretized model parameters
Here, we describe how the various parameters introduced above can be related to the experimental parameters. For the 40 K atoms with ω r = 2π × 69 kHz and ω z = 2π × 256 Hz [4], a ⊥ = 86 nm, and at a 3D Feshbach resonance (a 3D s → ∞), we have a 1D s = 62.8 nm.
If there are N atoms per spin, the Fermi energy for the noninteracting case is E F = N ω z . For the same number of atoms per spin, the dimensionless strength of interaction ξ at the Fermi level is Similarly, for n atoms per spin per unit length l, the Fermi wave vector is k F = πn/l, and we obtain

Method: Density-matrix renormalization group
After the discretization described above, the Hamiltonian is given bŷ where V(i) = κ(i − C) 2 with κ ≡ 4A/L 2 and C ≡ (L − 1)/2. Note that A/J ∝ L −2 (thus κ/J ∝ L −4 ) and U/J ∝ L −1 . We use DMRG to calculate the ground state of the Hamiltonian (8) within the particle-number sector of N ↑ and N ↓ atoms in spin-up and spin-down states, respectively. The population imbalance parameter is defined as P ≡ (N ↑ − N ↓ )/N total , where N total ≡ N ↑ + N ↓ . The Fermi momentum at the trap center is calculated from the averaged density n σ ≡ n i,σ /d as k Fσ ≡ πn σ , where the overbar denotes the average over 0.1L − 0.2L neighboring sites. Typically m = 300 states per block are retained and 4-5 finite system iterations are required for DMRG convergence after the standard infinite system initialization or the recursive sweep method initialization [53]. We target the ground state with N ↑ and N ↓ atoms in spin-up and spin-down states as well as several states with a few atoms added or removed, in order to improve the convergence of the pair correlation. While we do not assume the reflection symmetry of the wavefunction in the finite system iteration process, the reflection symmetry emerges spontaneously for the ground state wavefunction for any of the set of parameters studied. The sum of the eigenvalues of the reduced density matrix for the discarded eigenstates in each DMRG step is below 2 × 10 −5 for N = 64 and L = 320; it is below 5 × 10 −7 for N = 40 and L = 200 if we target the ground state alone.

Results: the Larkin-Ovchinnikov quasi-condensate 2.3.1. Pair correlation function We first examine the on-site pair correlation function defined as
where O on−site (z j , z j ) = n j,↑ n j,↓ gives the double occupancy ratio. The left column of figure 1 shows O on−site (z j , z j ) in color-coded two-dimensional plots for several different values of the imbalance parameter P. For P = 0, the pair correlation is maximal for z i = z j with the double occupancy per site, decreases with increasing |z i − z j |, and precipitously drops to zero where the atomic density, shown in the right column, vanishes.
As we increase P, the correlation function starts to oscillate. The zeros of the pair correlation O on−site (z, z j ) for fixed z j appear almost periodically, except for the region |z − z j | ≪ l, where O on−site is positive and large. The sign change of the pair correlation reflects that of the pair amplitude. We note that the whole region with a significant density of minority atoms is strongly correlated inside the 1D trap.

Density distribution: Oscillation of density difference
Next, we examine the density distributions of the spin-up and spin-down atoms. For the non-interacting system, the extensions of the majority and minority distributions are determined by those of the N ↑ -th and N ↓ -th eigenstates of the trap, respectively. As |U| is increased, the attractive interaction causes the atomic distributions to shrink towards the trap center, as shown in figure 2. We also find that the density difference n i,↑ − n i,↓ /d is initially spread all over the region with a nonzero atomic distribution, but that it shrinks more rapidly than the atomic distribution for sufficiently large |U|. The oscillation is independent of the lattice constant d for large L, as can be seen in figure 3. The peaks in the frequency spectrum D(q) of the density difference, calculated as and shown in the right column of figure 3, do not show a significant shift upon increasing the system size. From the spin-dependent Fermi momentum at the trap center, k Fσ , we obtain the difference of the Fermi momenta, q ≡ k F↑ − k F↓ . In the LO state, an up-spin atom close to the up-spin Fermi point (in 1D) ±k F↑ and a down-spin atom close to the down-spin Fermi point ∓k F↓ form a pair; therefore the total momentum of the pair must be close to ±q. In figure 4 we plot the cosine transformations of the pair correlation O on−site (z, z L/2 ), where z L/2 = 1/L is the location of one of the two sites closest to the trap center z = 0, and the density difference 2S z (z) ≡ n ↑ (z) − n ↓ (z) against the wave vector. To clearly identify oscillations close to the trap center, the density difference has been multiplied by a slowly varying numerical factor, ζ(z) ≡ exp (−6z 2 ), before Fourier transformation. The Fourier components of the pair correlation and density difference at wave vector q are defined as As shown in figure 4, the oscillation wave vector of the pair correlation function is always close to q. This is consistent with studies of population imbalance in the optical lattice systems, namely, the pair momentum distribution function calculated in [13], Fourier transform of the pairing correlation [14], on-site pair amplitude [15], spatial noise correlations [16], and radio-frequency spectroscopy calculation [20], as well as the pair momentum   distribution calculated in [22,23] for the continuum. The oscillation wave vector of the density difference is close to 2q, as expected in the LO state [49]. This result is consistent with the spin-spin correlation function studied in [24] and distribution of spin polarization studied by solving BdG equations in [21].

Phase diagram: Parameter region with phase separation
We find the following two phases: i) the spin-down (minority) atoms occupy a smaller region compared with the spinup (majority) atoms and the periphery of the system is occupied by the spin-up normal component alone, and ii) the densities of spin-up and spin-down atoms are different at the trap center but become almost equal over a finite width close to the periphery. These two situations can be thought of as two "phases" of the system. The transition between the two phases occurs when the trap is filled with the LO condensate. Such a phase transition was predicted in studies in which the exact solution of the Gaudin-Yang model without a trapping potential was plugged into the trapped system via the local-density approximation (LDA) [26,27]. For a given scattering length a or coupling constant g, case i) is realized for greater spin polarization P, where the condensate can no longer hold all the majority atoms and the system gains energy by pushing some of the majority component to the edge of the trap. This is because the condensation energy of the remaining LO condensate increases by a decrease in the Fermi wave-number difference and the number of nodes of the pairing potential.
Because the coupling constant |g| increases with increasing |U|, which corresponds to shorter |a| in 1D, the critical point of phase separation shifts toward greater values of P, as shown in figure 5. The phase separation curve, however, does not exceed P ≃ 0.21.

Effect of asymmetric potential: Accuracy of local density approximation
To study the robustness of the Larkin-Ovchinnikov state with respect to the trap geometry and the slope of the trap, we examine the ground state of the system for a spatially asymmetric potential, where 0 < C < L − 1, κ + = A/C 2 and κ − = A/(L − 1 − C) 2 . The bottom of the trap is located at i = C ≡ L √ r/( √ r + 1). Let r ≡ κ + /κ − be the ratio of the spring constants of the harmonic trap; this ratio equals the ratio between the potential slopes in the right and left sides of the bottom of the trap. Figure 6 shows the density distribution against the shifted trap coordinate z − z 0 , where z 0 = 2C/L − 1 = ( √ r − 1)/( √ r + 1) and V(z 0 ) = 0. Figure 7 shows the same density distribution plotted against the site energy level V(z).
Despite the increased asymmetry of the trap, the atomic distribution changes very little, and the pair correlation exhibits similar oscillations as those shown in figure 8. This is yet   another evidence that the LDA is a rather good approximation in a system of 1D populationimbalanced fermions. Recently the phase separation between the FFLO condensate and fully paired or fully polarized wings of the trap has been studied in [19], where a completely asymmetric setup, which corresponds to r → ∞, was used and compared with analytic results based on the combination of the Bethe ansatz solution and LDA.

Two-body pair correlation matrix
To quantify how much atoms contribute to the (quasi-) condensation in the populationimbalanced 1D system, we examine the behaviour of the eigenvalues of the two-body reduced density matrix, from which we can identify the presence or absence of the off-diagonal longrange order [54]. We consider the two-body reduced density matrix where α, β, γ, δ are one of the M possible states that N fermions in the system can each take. As proved in [54], the upper bound of the eigenvalue λ 2 of ρ full αβ;γδ is N(1 − (N − 2)/M), and the system possesses an off-diagonal long-range order if and only if λ 2 is of the order of N.
Since pairing occurs between up-spin and down-spin atoms, it suffices to calculate the part of ρ full for which α and β as well as γ and δ have different spins. Therefore, we calculate The eigenstates of this matrix with large eigenvalues correspond to states with high pair concentration. It can be shown that they become independent of the lattice discretization in the limit of L → ∞. Note that while the dimension of the matrix is L 2 , the trace of the matrix, which is also the sum of the eigenvalues, equals N ↑ N ↓ ≪ L 2 , and the maximum possible eigenvalue, which is close to N ↓ , is realized when every spin-down atom is paired with a spin-up atom and when these pairs are all condensed into a single state. We have diagonalized the two-body density matrix for the ground state that is obtained in the DMRG simulation for each parameter set. Because the number of independent matrix elements scale as O(L 4 ) with increasing the number of sites L, we used L = 80 and limited N up to 40. The number of states per DMRG block is also limited to m = 100, because for L = 80, the energy and largest eigenvalues of the two-body density matrix for the ground state converges for smaller m compared to the L = 200 case. The i = i ′ , j = j ′ components of this matrix are the on-site pair correlation functions and correspond to the pair amplitude of the condensation of on-site "molecules." In Ref. [13], the matrix of the on-site pair correlation, was diagonalized and the resulting eigenfunctions showed periodic oscillations. A similar feature is confirmed in our analysis of the full two-body pair correlation matrix. In figure 9 we plot the largest eigenvalue of ρ against the imbalance parameter P for different numbers of atoms N. The dependences on P for different values of N are strikingly similar. This indicates that, in the trapped system, the lowest energy state of atom pairs is occupied only by a few pairs, before N reaches 20, and they do not accommodate additional pairs of atoms, probably due to the effect of the trap potential which localizes the low energy states.
In the left column of figure 10 we plot the n-th largest eigenvalue of ρ against n. While the distribution is by definition a decreasing function of n, we find a kink in the distribution at different values of n for different P. Without the interaction, the eigenvalues of the reduced density matrix are always unity or zero, because each single-atom state is either filled or empty. The kink in the distribution suggests that the eigenvalues larger than the kink value are enhanced by the effect of pair quasi-condensation. In the right column of figure 10 we plot the sum of such enhanced eigenvalues against P. The sum is close to the number of minority atoms in the system, which indicates that almost all the minority atoms contribute to the quasi-condensation and the oscillating pair correlation.

Dynamics
Now we study the dynamics of a 1D, harmonically trapped Fermi gas after a population imbalance is introduced by flipping the spin of a single atom. The flipping can be realized by applying a microwave radiation, whose energy matches the hyperfine energy level separation.

Formulation: Real-time DMRG simulation after a spin-flip excitation
We simulate a sudden spin flip by applying a linear combination of local spin raising operators under spatial discretization. We obtain the wavefunction of the system for the ground state in the N ↑ = N ↓ = N-spin sector |Ψ 0 as discussed in the previous section, and then applyŜ + to obtain the normalized (with a constant N) wavefunction as the initial state of the spin-flipped system, The spin flip does not in general map the ground state in the (N ↑ , N ↓ ) = (N, N) sector to an eigenstate in the (N ↑ , N ↓ ) = (N + 1, N − 1) sector. Therefore, the system will evolve in time, showing some real-time dynamics. As a realistic case, we consider a situation in which the coefficient c l is spatially localized, which corresponds to the local spin excitation. This can be implemented experimentally by using a magnetic field gradient.
In the case of an attractively interacting system, the spin flip leads to pair-breaking and increases the energy of the system relative to the ground state in the (N+1, N−1). As discussed above, this ground state consists of a Larkin-Ovchinnikov condensate at the trap center and an equal-population condensate near the trap fringes, due to the low (P = 1/N < P c ) spin polarization.
For the exact calculation (after the spatial discretization), the full Hilbert space has to be considered, and the size of the Hilbert space rapidly increases. For example, for N = 4 in the L = 16-site system, the dimension of the Hilbert space is already 2446080. Several methods to apply DMRG to simulate real-time dynamics of quantum systems have been implemented [55,56,57,58,59,60,61,62]. Comparison of some of these methods have been conducted in [63]. From the ground state in the restricted Hilbert space obtained in the DMRG simulation, we calculate the spin-flipped state, and take this state as the state at t = 0. Then we apply the time-stepping targeting method [61] to simulate the real-time evolution of this initial wavefunction |Ψ(t = 0) .
In the time-stepping targeting method [61], the Schrödinger equation is treated as the equation of motion for the many-body wavefunction |Ψ(t) , and the Runge-Kutta method is used to obtain |Ψ(t + ∆t) , where ∆t is the time step, within the Hilbert space selected in the previous finite-size system DMRG step. At each step, |Ψ(t) , |Ψ(t + ∆t) , and several other wavefunctions |Ψ(t + a j ∆t) are targeted in the choice of basis. After one iteration of the finite-system algorithm DMRG, t is increased by ∆t and the process is repeated. Following [61], here we target four time steps |Ψ(t) , |Ψ(t + ∆t/3) , |Ψ(t + 2∆t/3) , and |Ψ(t + ∆t) , whose weights in the density matrix are 1/3, 1/6, 1/6 and 1/3, respectively. In this method, as in the methods based on the Suzuki-Trotter decomposition [57,58], the reduced Hilbert space in the subsystems after several DMRG iterations no longer contains much of the information on |Ψ(t = 0) . Use of the Runge-Kutta method in simulating realtime dynamics has also been discussed in [57].
Here we note that the dynamics of the superfluid to Mott insulator transition [64,65], spin-charge separation of a 1D Fermi gas [66,67], dynamics of spinless fermions after an interaction quench [68,69], expansion of a (bosonic or fermionic) Mott insulator [70,71,72], relaxation of atoms after a superlattice is modified [73,74], quench dynamics of the Bose-Hubbard model [75], and damping of dipole oscillations of a 1D Bose gas after a sudden displacement of the trap [76] have been studied by the time-dependent DMRG or a related method, time-evolving block decimation (TEBD) [77]. For a mixed state, the real-time evolution of the density matrix of the system can be simulated in DMRG. Simulations of the dynamics for finite-temperature spin systems have also been carried out [78,79]. Here, we show the results of our simulation for the dynamics of pure states at T = 0.

Results: Diffusion of the excess majority atom and destroyed correlation
In the following, a spin-flipping excitation is parameterized by two parameters, the center of excitation z c and width w. The excitation factor for each site l, s l , is determined by where C is the normalization constant. First, we examine the case in which the spin is flipped close to the trap center (z c = 0). We consider the initial state obtained by flipping one of the down spins from the ground state in the (N ↑ , N ↓ ) = (8,8) sector, for U = −8 and A = 1. We take L = 32 and w = 1/16 and retain m = 300 states per DMRG block. Here, J = 2 L 2 /(8m 0 l 2 ) = 39.4 nK × k B , and we take /J = 194 µs as the unit of time. In figure 11 we plot the evolution of the energy, defined by where |Ψ(t) is the calculated wave function at time t after the spin flip at t = 0. We observe that the energy goes up rapidly for ∆t = 0.01 after t ≃ 1.3. Because the Hamiltonian is always the same, this increase in energy is due to the accumulated error in the Runge-Kutta numerical integration. For ∆t = 0.002, the increase of the energy occurs later and it is less significant. For ∆t = 0.001, the energy is almost constant until t ≃ 4. We present the result for ∆t = 0.001 up to t = 4.
For U = −8, we plot the density distribution for different time steps in figure 12. Because of the spin flip, initially (at t = 0) the majority density is largest close to the trap center, and the minority density has a deep dip there, making the density distribution localized there. As the time elapses, the majority density decreases and the minority density increases at the center, and the density difference exhibits two peaks at both sides of the central, decreasing peak. Up to t = 4, these two peaks move away from the trap center. There remain three peaks of the density difference. This is to be compared with the density distribution for the (N ↑ , N ↓ ) = (9, 7) sector of the Hamiltonian (8) (see figure 13). The majority and minority density profiles are not as smooth as those after the spin flip. Such oscillations in the density profiles have also been observed for larger L; they originate from the formation of tightly bound pairs, rather than the lattice discretization [15]. In the excited state after the spin flip, such oscillations in the ground state are smeared out because of the weaker binding between spins and the excess kinetic energy.
In figure 13, we see two peaks of the density difference that correspond to the pair amplitude nodes of the LO condensate. While the profile at t = 3 after a spin-flip excitation in figure 12 might seem similar, the peak of the density difference continues to move outward and the profile at t = 4 has three peaks (including the one at z ≃ −0.45). While the local spin excitation develops into an oscillating density difference profile, the oscillation does not guarantee an LO condensate.
Next we examine the effect of the location of the spin flip. In figure 14, we plot the density profile after the spin flip at two different z c off the center, z c = −0.3125 and −0.625, respectively. The resulting density difference profile at t = 0 does not have a peak at z c ,  because in calculating the excited state (18), the amplitude of the spin-flipped wave function S + l |Ψ 0 , Ψ 0 |Ŝ − lŜ + l |Ψ 0 , depends on the site number l. As the state evolves in time, the peak splits into two and they move to the left and right, as in the case of the spin flip at z = 0. However, the movement of the density difference peak to the direction of z > 0 is considerably faster than that to the direction of z < 0, as the kinetic energy of the atoms is greater for z ∼ 0.
In figure 14, as the spin flip location is removed from the trap center, we observe that the change of the minority population in the side of the spin flip (z < 0) becomes less pronounced, and that in the other side of the trap, the densities of both components stay almost the same. This is because the spin-imbalanced gas in the z < 0 side pushes the spin-balanced region (z > 0) away, thus compensating for the reduced kinetic energy by increasing the potential energy in the trap, but infiltrating this region only slowly.
Finally we look at the effect of the interaction strength. For a smaller |U|, which corresponds to a smaller |g 1D | and a larger |a 1D s |, an increase in energy due to the spin flip is smaller, but at the same time, the effect of Fermi pairing might become weaker.
In figure 15 we show the evolution of the system for U/J = −2 and A/J = 1, after the spin flip at the trap center (z c = 0). Observe that, while the atoms are extended over a wider region due to the weaker attraction, the speed of the movement of the two density-difference peaks is almost the same as in the case of U/J = −8. The insensitivity of the speed to the interaction strength, which is also observed for other spin flip parameters, suggests that the momenta of the excess up-spin atoms are determined primarily from the shape of the center peak whose width is independent of |U|, and not much affected by U, which determines the pairing strength. Indeed, if |U| is kept constant and w is increased, the width of the initial distribution of the population imbalance increases, and the speed of the diffusion of the density difference rapidly decreases; this suggests that the speed may be determined from the uncertainty principle in the initial peak where the excess up-spin atoms are introduced by the spin flip.

Conclusion
We have studied the ground state and dynamics of harmonically confined 1D Fermi gases with population imbalance. We have discretized the 1D system and applied the densitymatrix renormalization group method, which is numerically exact at least for the ground state in our setup, and also allows the simulation of the dynamics of large quantum systems at an unprecedented precision. For the ground state, we have demonstrated that for finite imbalance and attractive interaction, a Larkin-Ovchinnikov type condensate forms at the trap center, for a wide range of the polarization parameter up to almost unity. Qualitatively this result is consistent with both numerical studies [13,14,15,20,21,22,23], and studies that utilize the exact solution for a uniform system together with local density approximation [26,27,28,29,30].
Beyond a critical value of P, the LO condensate close to the trap center phase-separates from the spin-polarized Fermi gas at the periphery, regardless of the strength of interaction. By studying the dependence of the eigenvalues of two-body density matrix on P and N, we have shown that the oscillating pair correlation originates from quasi-condensation of atom pairs, to which almost all the minority atoms contribute. The FFLO-like feature is found to be rather insensitive to the asymmetry, which again supports the validity of the local density approximation.
We have studied the dynamics after a spin-flip excitation by time-dependent DMRG simulations. Of particular interest is whether the density difference oscillates as in the LO condensate in the ground state. We have found that, while the snapshots of the density difference show two or more peaks that originate from a single peak after the localized spin flip, the peaks move in time and are not related with an FFLO-like state.