Creation and dynamics of remote spin-entangled pairs in the expansion of strongly correlated fermions in an optical lattice

We consider the nonequilibrium dynamics of an interacting spin-1/2 fermion gas in a one-dimensional optical lattice after switching off the confining potential. In particular, we study the creation and the time evolution of spatially separated, spin-entangled fermionic pairs. The time-dependent density-matrix renormalization group is used to simulate the time evolution and evaluate the two-site spin correlation functions, from which the concurrence is calculated. We find that the typical distance between entangled fermions depends crucially on the onsite interaction strength, and that a time-dependent modulation of the tunneling amplitude can enhance the production of spin-entanglement. Moreover, we discuss the prospects of experimentally observing these phenomena using spin-dependent single-site detection.


Motivation
The tremendous experimental progress with ultra cold atoms in optical lattices makes them a unique playground for studying nonequilibrium dynamics of many-body systems. The remarkable features of these systems are the precise and dynamical control of the interparticle interaction and external trapping potentials, as well as long coherence times [1]. Different nonequilibrium initial states can be generated by a quantum quench (for a recent review see [2]): a parameter in the Hamiltonian is suddenly changed such that the system, initially in the ground state of the Hamiltonian, is afterwards in an excited state of the new Hamiltonian. Quenches have been experimentally realized, for instance, in the interaction strength [3][4][5] and in the additional trapping potential by either switching it off [6] or displacing its centre [7][8][9]. This led to observations such as the collapse and revival of the coherence in a Bose-Einstein condensate after a quench from the superfluid to the Mott insulating regime [3].
New detection schemes [10][11][12] provide the possibility of observing the many-body state at the single-site and single-atom level and make these systems even more suitable for studying the dynamics of (spatial) correlations in nonequilibrium situations. These methods have already been used to monitor the propagation of quasi-particle pairs [13,14] and a single spin impurity [15] in a bosonic gas. They have also inspired theoretical work on many-body dynamics subject to observations, treating issues such as entanglement growth in a bosonic system [16], destructive single-site measurements [17] and quantum Zeno effect physics for spin [18] and expansion [19] dynamics.
In the present paper, we study a different aspect of these systems: the creation and dynamics of spin-entanglement during the expansion of a strongly interacting, spin-balanced fermionic gas in an optical lattice. We find that the expansion out of an initial cluster of fermions can automatically generate long-range spin-entanglement.
Expansion dynamics of interacting fermions in an optical lattice has been realized recently in a three-dimensional optical lattice [6]. In that experiment, the authors observed a bimodal expansion with a ballistic and diffusive part and were able to explain its anomalous behaviour using a Boltzmann-based approach. In addition, expansion dynamics of this kind has been studied numerically in one dimension, going beyond a kinetic description. These studies revealed the dependence of the expansion velocity [20,21], the momentum distribution function, and the spin and density structure factors [22] on the onsite interactions and the initial filling. Furthermore, the effects of the different quench scenarios [23] and of the gravitational field [24] on the time-evolution of the density distribution have been discussed. The sudden expansion of a spin-imbalanced fermionic gas has been recently considered for observing Fulde-Ferrell-Larkin-Ovchinnikov correlations [25][26][27].
In general, the dynamics is crucially affected by the difference in behaviour between a single fermion and that of a doublon (i.e. a pair of fermions at the same lattice site). For large interaction strengths doublons are very stable against decay into fermions (analogously to the repulsively bound boson pairs [28]) and move slowly. These properties can lead to the condensation of doublons [29] and a decrease of the entropy [30] in the centre of the cloud during the expansion. They play also a role in the decay dynamics of doublon-holon pairs in a Mott insulator [31][32][33]. In the present work, we will show how spin-entanglement is generated by the decay of a doublon into single fermions.
While we focus here on the build-up of correlations in a many-body state (which is mainly driven by the creation of correlated fermions out of doublons), the efficient production of entangled atom pairs is by itself an important topic, especially in the context of atom interferometry [34]. Using such nonclassical atom pair sources would allow matter-wave optics beyond the standard quantum limit. Recent experiments succeeded in generating large ensembles of pair-correlated atoms from a trapped Bose-Einstein condensate, using either spinchanging collisions [35,36] or collisional de-excitation [37]. In the context of such experiments, detection methods able to image single freely propagating atoms have also recently become available [38].
The paper is organized as follows. First, we will describe in detail the expansion protocol, starting from a finitely extended band insulating state, i.e. a cluster of spin-singlet doublons. Moreover, we give details about the numerical time-dependent density matrix renormalization group (tDMRG) simulation. Before studying spin-physics, we first discuss correlations in the density of the expanding cloud (section 2), where the influence of interactions is already significant. We find that large onsite interactions lead to positive correlations between certain lattice sites. The dynamics of spin-entangled pairs, the main focus of our work, is then presented in section 3. First of all, we relate the concurrence to the spin-spin correlation functions. Then we discuss the propagation of spin-entangled pairs in the cloud. Furthermore, we compare the cumulative 'amount' of spin-entanglement in different regions of the cloud and for various interaction strengths. We find that spin-entanglement between remote lattice sites is only found for large interaction strength, while for small onsite-interactions there is entanglement preferentially only between nearby sites. Moreover, the Pauli-blocked core of the cluster favours both partners of a spin-entangled pair to be emitted into the same direction when compared to the decay of a single doublon, where they almost always are emitted into different directions. In addition, we discuss the expansion under the action of a time-dependent (modulated) tunnelling amplitude. We find that the modulation can enhance the production of spin-entanglement. Finally, we discuss a possible experimental setup for observing this spin-entanglement dynamics in the near future using spin-dependent single-site detection (section 4).

Model
In this paper, we consider an ultracold gas of fermionic atoms loaded into a one-dimensional optical lattice. The atoms can be prepared in two different hyperfine states, which we label ↑ and ↓, and fermions of different 'spin' interact via s-wave scattering. This fermionic mixture represents a realization of the standard fermionic Hubbard Hamiltonian [39,40] The first term describes the tunnelling of fermions between adjacent lattices sites with amplitude J . The second term encodes the effective onsite interaction energy U of fermions with different hyperfine states, which can be controlled using a Feshbach resonance. The particle number operator isn i,a =ĉ † i,aĉ i,a , with the creation (annihilation) operatorĉ † i,a (ĉ i,a ) satisfying the fermionic commutation relations.
In the following, we focus on the expansion from a band-insulating state, i.e. the sites in the centre of the lattice are doubly occupied. Experimentally, this is achieved by using an additional trapping potential to confine the fermions to a central region of the optical lattice [6]. At time t = 0 the trapping potential is switched off and the fermions expand into the empty lattice sites, as depicted in figure 1(a).
The time evolution of the average fermion density and the average doublon density has already been studied for this expansion protocol, by numerical simulations for 1D lattices [20]. For large onsite interactions strengths U/J 4 two wave fronts with different velocities are visible in the time-dependent density profile, while there is a single wave front for small onsite interaction strengths. It turns out that the expansion can be basically understood as a mixture of propagating single fermions and doublons, see also [21]. For small onsite interaction strengths the initial state quickly decays into single fermions, which move ballistically through the lattice. Due to the cosine dispersion relation of a single particle in the lattice, k = −2J cos(k) with wavenumber k ∈ (−π, π], its maximal velocity is given by |v max | = 2J . This fact leads to light cones in the density distribution, as indicated in figure 1(a). For large onsite interaction U/J 4, only a small fraction of doublons decay into fermions, which move away ballistically. As the effective tunnelling amplitude of a doublon is 2J 2 /U , they initially remain in the central region and then move through the lattice on a much larger time scale.
However, the propagation of the single fermions is highly correlated, as they are always created as spin-singlet pair when a doublon decays. As sketched in figures 1(b) and (c), the two particles of such a spin-entangled pair may be either detected on the same side of the initial cluster or on different sides.
We briefly point out that the dynamics of the spin-entanglement as well as that of the density-density correlation is invariant under the change of the sign in the interaction strength U. This is a direct consequence of a transformation property of the (spin-or density-) correlators and of the initial state under time reversal and π -boost (translation of all momenta by π). This dynamical U → −U symmetry in the Hubbard model is discussed in more detail in [6,41].

Numerical simulation
For the numerical evaluation of the nonequilibrium time evolution, we employ the tDMRG method [42][43][44][45][46]. The initial state is a cluster consisting of ten doublons, which are located in the centre of a lattice of size L = 100 with open boundary conditions. Our tDMRG simulation uses a Krylov subspace method with time step J δt = 0.125 and the discarded weight is set to either 10 −5 or 10 −6 , depending on the interaction strength. For all evolution times shown in the figures, the density remains negligible at the boundaries of the lattice. We also verified that the features presented here does not depend on the precise position of the cluster in the lattice nor change when the truncation error is modified within the given range.
For noninteracting particles, we have used the exact expressions for the time-dependent correlation functions and concurrence, as given in appendix A.

Density-density correlation
Before turning to the spin-entanglement (in the subsequent section), we will first study the density-density correlation function D i j (t) = n i (t)n j (t) − n i (t) n j (t) , withn i =n i,↑ + n i,↓ . For the expansion, sketched in figure 1, the density at different lattice sites is expected to be correlated for several reasons, in particular since the fermions are always created in pairs out of doublons.
Note that the single-site detection of density correlations (more precisely the parity correlations) has been very recently used to study the quasi-particle propagation in a commensurate ultracold bosonic gas after an interaction quench [13,14]. It has also been employed to study the role of onsite interaction in a bosonic two-body quantum walk, experimentally realized in a nonlinear optical waveguide lattice [47].
As a point of reference, we first consider the density-density correlations for the initial state consisting of a single doublon. The fluctuation in the density at a lattice site i, D ii (t), is the variance of the occupation numbern i . For the expansion from a doublon, it is maximal for those lattice sites located at the edge of the single fermion light cone (doublon light cone) in the noninteracting (strongly interacting) case. In the following we focus on the off-diagonal density-density correlations D i j =i (t), which are shown in figures 2(a) and (b). Most importantly, we find that the onsite interaction U leads to a positive correlation of those fermions propagating with almost maximal velocity |2J | in opposite directions, see figure 2(b). On the other hand, the density-density correlation assumes large negative values between those lattice sites in the centre that are most likely occupied by the doublon (and not by single fermions) (central square region in figure 2(b)). In contrast, the off-diagonal density-density correlations are always nonpositive for clusters of noninteracting fermions, as detailed in appendix A, see also figures 2(a) and (c)-(e).
The bunching effect in the interacting case can be understood by writing the motion of the two fermions in relative and centre of mass (c.o.m.) coordinates, r and R, respectively (see appendix B for details of the discussion). Note that we assume an infinite lattice for this argument, which is compatible with the simulation as the boundary conditions do not play a role for the evolution times considered here. The c.o.m. motion is a plane wave with total wavenumber K = (k 1 + k 2 ) mod 2π , where k 1 and k 2 are the asymptotic wave numbers of the single fermions. The relative motion is described by a K -dependent Hamiltonian. For U = 0, the eigenstates of this Hamiltonian are one bound state and scattering states. The probability for the doublon to decay into one specific scattering state with wavefunction ψ K ,k (r ) (where k = (k 1 − k 2 )/2 is the relative wavenumber) is given by the modulus squared of that wavefunction's probability amplitude at r = 0: |ψ K ,k (0)| 2 . This decay probability is, up to an overall normalization constant (see equation (B.5)): For small onsite interactions |U/4J | 1, the decay probability is almost the same for the different scattering states, except for K ≈ π or k ≈ 0 where it drops to zero. In the strongly interacting case, it exhibits a pronounced maximum for K = 0 and k = ±π/2, that is k 1 = ±π/2 and k 2 = ∓π/2. In other words, for large onsite interaction the doublon decays primarily into two fermions moving in opposite directions with velocity close to |2J |.
Next we study the density correlations for the expansion from a cluster of several doublons into an empty lattice. The results are displayed in figures 2(c)-(h). Just as for the single doublon, the off-diagonal density correlation has regions of positive values in the presence of onsite interaction. However, in contrast to the case of a single doublon, the existence of other doublons now leads to positive correlations also at lattice sites located on the same side of the initial cluster, see figures 2(f)-(h). Positive correlations between sites on different sides of the cluster are only observed for evolution times larger than the time a hole takes to propagate through the cluster, cf figure 2(h).
The results suggest that the initially Pauli blocked core of the cluster leads to an enhanced decay of the outermost doublons into single fermions moving away in the same direction. However, the alternative case, where the edge doublon decays into a fermion and a hole moving into the opposite direction close to the maximal velocity |2J | also leads to positive correlations. Further simulations show that positive density correlations between lattice sites on the same side of the initial cluster positions are found for clusters of about four and more doublons and U/J 2. Note that as the single fermions are created by the decay of an edge doublon, the situation is different from a free single fermion or fermionic wave packet approaching a cluster of doublons. In the latter situations, the fermion is almost perfectly transmitted through the cluster (band insulating state) in the limit of large onsite interaction [48,49].
In the next section, we discuss the spin-entanglement. This will reveal that the positive density correlations indeed stem from singlet pairs, which become delocalized by the decay of a doublon.

Spin-entanglement between different lattice sites within the expanding cloud
In this section, we discuss the entanglement between fermions located at different lattice sites by means of the concurrence [50]. First, the relation between the concurrence and the spin-spin correlations is established. Then we use tDMRG simulations to find the regions where fermions are likely and unlikely to be entangled during the expansion. We examine the role of the onsite interaction and the core of the cluster on the entanglement dynamics.

Reduced density matrix and concurrence of two fermions
In this paragraph we derive the reduced density matrix and the concurrence for fermions located at different lattice sites within the expanding cloud (see [51] for a related discussion in the context of solid state physics). Given the full many-body wavefunction | (t) , the reduced density matrix for the lattice sites i and j isρ i j (t) = Tr sites =i, j {| (t) (t)|}, where all other lattice sites have been traced out. Whenever two fermions are situated at the same lattice site they form a spin-singlet pair as their spatial degrees of freedom are symmetric under particle exchange. Thus, entanglement in the spin degree of freedom between fermions at two different lattice sites can only occur if each lattice site is occupied by a single fermion. Projectingρ i j (t) onto those states yieldŝ Here,n s i =n i,↑ +n i,↓ − 2n i,↑ni,↓ is the single fermion number operator at site i, which projects onto a subspace with exactly one fermion on that site. The normalization factor in the denominator, Tr{ρ i j (t)n s in s j } = n s i (t)n s j (t) , is the probability of finding at time t a single fermion at each of the lattice sites i and j. The state described byρ s i j (t) is the state obtained after a successful projective measurement. The reduced density matrixρ s i j (t) is equivalent to a two-qubit density matrix, which can be expressed in the formρ = 3 α,β=0 λ αβσ α (1) ⊗σ β (2) , wherê σ 1,2,3 are the Pauli matrices,σ 0 = 1, and the factors λ αβ are determined by the correlation [52]. Consequently, the reduced density matrix of single fermions at lattice sites i and j can be written aŝ is the x-, y-, and z-component of the spin operator and S 0 i is, for compactness, defined as half the single fermion number operator,Ŝ 0 is calculated using the full (unprojected) wavefunction since states with vacancies or doublons at site i or j do not contribute to the expectation value.
Symmetries of the initial state and the Hamiltonian can simplify the form of the reduced density matrix, such that only a few correlation functions are needed to determineρ s i j (t). As detailed now, the reduced density matrixρ s i j (t) depends only on Ŝ z i (t)Ŝ z j (t) / n s i (t)n s j (t) for the cluster initial state shown in figure 1. The Hubbard Hamiltonian (1) preserves the spin-dependent particle number, i.e. [Ĥ,N ↑,↓ ] = 0,N ↑,↓ = L i=1n ↑,↓ . Given an initial state with fixed number of spin-up and spin-down fermions, which is the usually case in cold atom experiments, the time-dependent expectation values of operators that do change the spin-dependent particle number vanish. This yields n s i (t)Ŝ x,y The latter condition comes from creating two spin-down (spin-up) fermions while destroying two spin-up (spin-down) fermions at lattice sites i and j and can also be written as Ŝ x,y,z i ] = 0. If the initial state is rotationally invariant, for instance a cluster of doublons, then the many-body state remains SU(2) spin symmetric during the time evolution. It In summary, the reduced density matrix for the expansion from a cluster of doublons readŝ where is the singlet state and |T m i j is the triplet state with the m denoting the spin projection in z-direction.
Instead of Ŝ z i (t)Ŝ z j (t) one could in principle evaluate the spin-spin correlation in any other direction. However, for cold atomic gases in optical lattices the correlation function seems to be experimentally most realistic to access as it could be obtained from snapshots of the spin-dependent single-site detection of the particle number.
The spin-entanglement between single fermions can be derived from the reduced density matrix. In this work, we use the concurrence C(ρ) [50] to quantify the entanglement. Given the time-reversed density matrixρ = σ y ⊗ σ yρ * σ y ⊗ σ y , with the complex conjugationρ * taken in the standard basis {|↑↑ , |↑↓ , |↓↑ , |↓↓ }, the concurrence is defined by where the λ i are the eigenvalues ofρρ in descending order. In our case, the concurrence of the reduced density matrix (5) is given by Spin-spin correlations Ŝ z i (t)Ŝ z j (t) / n s i (t)n s j (t) −1/12 result in a vanishing concurrence. The concurrence approaches 1 as Ŝ z i (t)Ŝ z j (t) / n s i (t)n s j (t) −1/4, i.e. when the two fermions detected at lattice sites i and j always have opposite spin.
The probability that single fermions at lattice sites i and j form a spin-singlet pair can be directly read off the reduced density matrix (5): Each of the spin triplet states is measured with the probability 1/4 + Ŝ z i (t)Ŝ z j (t) / n s i (t)n s j (t) .

Time evolution of the spin-entanglement within the expanding cloud
Spin-entanglement between different lattice sites, which we denote by i and j, requires that both sites are singly occupied as discussed in section 3.1. The probability for this is n s i (t)n s j (t) , with the single fermion number operator already defined above (n s i =n i,↑ +n i,↓ − 2n i,↑ni,↓ ). Figure 3 shows the numerical results for n s i (t)n s j =i (t) for different evolution times J t and onsite interaction strengths U/J . Shortly after the quench, single fermions are created at the edges of the cluster, see figures 3(a), (d) and (g). For the noninteracting system, the fermions escape the cluster successively (i.e. starting at the edges, and finally from the centre of the cluster). They move ballistically and the correlation function displays a fourfold symmetric structure, figures 3(b) and (c). In the interacting case, in contrast, the decay of a doublon into fermions is heavily suppressed. This results in a correlation function n s i (t)n s j (t) that has its main contributions within two square regions given by the light cones of fermions emitted by the outermost doublons of the cluster, see figures 3(f) and (i). Within the light cones, single fermions are more likely to be found at lattice sites corresponding to a motion with almost maximal velocity |v max | = 2J into opposite direction. The density correlation between single fermions attains its largest values for nearest-neighbour lattice sites within the cloud, cf figures 3(e), (f), (h) and (i). This may be due to virtual transitions between two configurations: either a doublon with a neighbouring vacancy, or a state of two adjacent single fermions. Note that the mean total number of single fermions, L i=1 n s i , approaches within a few J t an almost constant value, cf area under the curves for J t = 3.5,7.5 and U/J = 3, 6 in the last column of figure 3. That means only the edges of the cluster evaporate for large interactions, releasing a finite number of single fermions. Moreover, n s i (t) becomes relatively flat in the centre of the cloud for larger evolution times J t.
Let us now consider the spatial distribution of spin-entangled fermions during the expansion. For this purpose, the concurrence between two fermions at different lattice sites is numerically evaluated using equation (6). Figure 4 shows the concurrence for any pair of lattice sites i and j, at different times J t and onsite interaction strengths U/J .
For the expansion of noninteracting fermions, figures 4(a)-(c), the concurrence is finite only for nearby lattice sites. It is almost 1 within the outermost wings of the expanding cloud. This can be physically understood the following way: due to the Pauli principle, fermions with the same spin become spatially antibunched during the expansion. Thus, fermions at neighbouring sites are more likely to have opposite spin, cf figures 2(c)-(e). In addition, the outermost region of the cloud lies only within the light cones of the doublons close to the edge of the initial cluster. Two fermions detected in this region are almost certainly emitted by the same edge doublon and in consequence have a high probability to be spin-entangled.
When increasing the onsite interaction U , spin-entangled pairs are formed on remote lattice sites, too, cf figures 4(d)-(i). Indeed, spin-entanglement is found between lattice sites within and outside the initial cluster position, figures 4(e) and (h), as well as on different sides of the initial cluster position, figures 4(f) and (i).
The figures are a fingerprint of the creation of a counter propagating hole and single fermion by the decay of an edge doublon. When the hole moves through the cluster, the spinentanglement with the single fermion outside the cluster is swapped sequentially from one fermion to the next in the cluster. In this way fermions become entangled that have never been on the same lattice site and have never directly interacted with each other. At the end of this process, a spin-singlet pair is created with a fermion at each side of the cluster. The concurrence for two fermions on different sides of the initial cluster position increases with the interaction strength, compare figures 4(f) and (i). This can be understood by the suppression of the decay probability of doublons with increasing interaction strengths. For larger interaction strength, a hole is more likely to cross the cluster without being disturbed by another hole.
For nearest-neighbour sites, which have a relatively large probability to be simultaneously singly occupied (see figures 3(f) and (i)), we observe a concurrence close to 1. This implies Concurrence C i, j (t) of two (single) fermions located at lattices sites i and j after different expansion times J t, see equation (6). Initially, ten doublons are located at the sites indicated by the dashed square. We emphasize that dark red indicates strictly zero concurrence (also see the cuts displayed to the right). The black colour code is used whenever the probability of finding fermions at sites i and j is too small, n s i (t)n s j (t) < 10 −5 . In those cases, the concurrence is not computed as it would become susceptible to numerical inaccuracies. (a), (d), (g) For short evolution times, here J t = 1, single fermions are mainly created at the edge of the cluster since the central fermions are initially Pauli blocked. Fermions close to the same edge of the cloud are likely to be entangled as they are likely to originate from the same doublon. Fermions at opposite edges are not entangled since they cannot have been emitted from the same doublon. (b), (c) At larger evolution times and without onsite interactions, the concurrence is 0 for most sites i and j. However, it is close to 1 when two fermions are located at the outermost part of the cloud. (e), (f), (h) and (i). By increasing the onsite interaction U/J , entanglement of fermions across the cluster becomes possible. In this case, the concurrence is highest when the fermions result from the decay of a doublon at the edge and escape with opposite velocity, indicated by the dotted line (see also cut through panel (i)). Moreover, the concurrence of fermions at neighbouring sites becomes almost 1. Within the central region, approximately given by the overlap of the light cones shown in figure 1, the concurrence remains relatively small. the creation of vacant lattice sites within the cloud during the expansion. The singlet pairs on nearest-neighbour sites come from virtual transitions between a doublon with neighbouring vacancy and a state of two single fermions. We verified this behaviour in addition by numerically creating an ensemble of snapshots for the distribution of fermions as described in [19].
We emphasize that for inhomogeneous systems, such as the one discussed in this paper, the spatially resolved measurement of two-point correlations provides more information about the dynamics than structure factors, such as S k ∝ lm e −ik(l−m) Ŝ z lŜ z m . While the latter could be used to determine the average spin-spin correlation of fermions at fixed distance, it contains no information where these spin-correlated pairs are located in the cloud.

Summed concurrences
Above, we examined the spin-entanglement between two lattice sites for fixed time points. In this subsection we aim to quantify the spin-entanglement for entire regions in the lattice. In particular, we address the questions Are there sites which share more spin-entanglement with the rest than other lattice sites? How does the spin-entanglement in different regions built up as function of time? Which locations are most entangled in the weakly and strongly interacting case? How does the size of the cluster affect these results?
In the following we discuss the amount of pairwise spin-entanglement of a lattice site or a region in the lattice in terms of the summed concurrence. We define it as the sum over the concurrences, C i, j (t), which are weighted by the probability of detecting a single fermion at both lattice sites, n s i (t)n s j (t) . The weights are introduced to accommodate for the possibility of vacant or doubly occupied lattice sites, in contrast to the summed concurrence used in spin systems [53]. For a system consisting of spin-singlet pairs whose wavefunctions do not overlap, i.e. C i, j (t) is either zero or one for all sites i and j, the summed concurrence equals the average number of delocalized spin-singlet pairs. Note that the summed concurrence is by no means a measure for the total entanglement of the system. It neither includes multipartite entanglement nor entanglement in the occupation numbers. For a detailed discussion on entanglement in many-body systems we refer to [54], see also [16,[55][56][57] for recent proposals on detecting entanglement in cold atom systems.
Let us consider the spin-entanglement of a site with all the other lattice sites. The summed concurrence for site i is defined by The time evolution of C tot,i (t) is shown in figure 5 for different interaction strengths U/J . For noninteracting fermions ( figure 5(a)), lattice sites close to the edge of the cloud display the strongest entanglement, while sites in the rest of the cloud are hardly entangled. For increasing interaction strengths a central region with almost uniform C tot,i (t) builds up during the evolution, see figures 5(b) and (c). For large U/J , we find that C tot,i (t) approaches the expectation value n s i . This turns out to be related to the fact that in this case the probability of having two doublons decay is negligible, and the contribution comes almost entirely from the decay of a single doublon.
In figure 4(i) it is apparent that onsite interactions can lead to spin entanglement across the expanding cluster. In the following we compare the summed concurrences of lattice sites on the same side and on different sides of the initial cluster location, C ss and C ds , respectively. They are defined by where l and r denote the leftmost and rightmost occupied lattice sites of the initial state. Note that nearest-neighbour lattice sites are excluded from C ss (t). That means we do not take into account contributions from virtual transitions of doublons (decaying virtually into two adjacent fermions) that move away from the cluster initial position. In addition, we consider the summed concurrence of sites at fixed distance d, , and the summed concurrence of all sites, C tot (t) = i< j n s i (t)n s j (t) C i, j (t). The time evolution of these summed concurrences is shown in figure 6.
For noninteracting fermions and small evolution times, the total summed concurrence is dominated by contributions from nearest-neighbour sites. For larger times more and more spinentanglement is transferred to fermions found on the same side of the initial cluster position (C ss ), see figure 6(a). The spin-entanglement remains relevant only for small distances, reflected in C 8 (t) = 0 for all simulated times, cf figure 6(d).
By contrast, in the interacting case, spin-entanglement is generated via fermions propagating away on different sides of the cluster. This is seen as a finite value of C ds (t) for times J t 5, which is the time a hole needs to propagate through the cluster, cf figures 6(b) and (c). The total summed concurrence and concurrence between nearest-neighbour sites quickly settle into a damped oscillation around a constant value. The time evolution of the summed concurrence of sites at distance d 2 is displayed in figures 6(e) and (f). C d (t) remains zero for small times, peaks at times J t slightly exceeding d/4, and decreases for larger times. This shows that spin-entangled pairs mainly propagate at almost maximal (relative) velocity 4J . For large U/J , C d (t) is approximately uniform for distances d 4J t and roughly decays as (J t) −1 for the simulated times. Note that C 1 (t) plays a special role. Its main contribution does not stem from 'free' fermions. Rather, it is generated by a doublon virtually dissolving into adjacent fermions. Let us finally discuss the impact of the cluster size on the spin-entanglement dynamics. For very weak interactions, a larger number of doublons means that more delocalized singlet pairs are created shortly after switching off the confining potential. For large interaction strengths up to U/J = 40, we simulate the expansion and compare the summed concurrences for different cluster sizes, including the case of a single doublon. We summarize the results in figure 7. Note that reasonably large evolution times (J t ≈ 20) for the comparison of the summed concurrences are reached only for U/J 6. The summed concurrence of all sites, C tot , The summed concurrences of lattice sites at same side (C ss ) and at different sides (C ds ) of the cloud (see the main text for the definitions) disagrees for a single doublon (dots) and a cluster of doublons (squares). Note that clusters of sizes 4, 7 and 10 give similar values. For strong interactions, a cluster prefers the emission of (delocalized) singlet pairs into the same direction compared to a single doublon. The summed concurrence between nearest-neighbour sites (C 1 ) approaches C tot /2 in both cases. Solid lines show analytical results for C 1 as well as the contribution of scattering states to C ss and C ds for a single doublon (see appendix C).
agrees for all considered cluster sizes and matches the analytical result for a single doublon, see figure 7(a). Apparently, the initially Pauli-blocked core has no effect on the number of created single fermions for the considered times. For C ss and C ds , however, we find a clearly different behaviour for a single doublon and a cluster of four and more doublons ( figure 7(b)): a cluster is more likely to emit (delocalized) spin-entangled pairs into the same direction.

Expansion with modulated tunnelling amplitude
In the previous section we have seen for the interacting case that total summed concurrence, C tot , approaches a fixed value shortly after switching off the potential, via the escape of a few fermions from the cluster edges. Here, we discuss a way of 'continuously' generating single fermions and enhancing C tot compared to the free time evolution. We consider an expansion during which the tunnelling amplitude is repeatedly varied in time, while the interaction strength is constant. Such modulation may be experimentally realized by either varying the laser intensity (the tunnelling amplitude decreases much faster with increased laser intensity than the onsite interaction strength, see e.g. [1]) or by shaking the lattice sinusoidally [58]. We find for certain values of the tunnelling amplitudes and time intervals between the quenches an increased Jt Jt

Remarks on observing the spin-entanglement in experiments
In the main part of this article we have analysed the dynamics of the spin-entanglement for the expansion from a cluster of doublons. As discussed in section 3.1 the concurrence between two lattice sites i and j can be determined by the single fermion expectation value n s i (t)n s j (t) and the spin-spin correlation Ŝ z . Experimentally, both expectation values could be obtained by averaging over many snapshots of the spin-dependent single-site fermionic particle number analogous to the already implemented single-site detection of bosonic particles [11,12]. Since doubly occupied lattice sites do not contribute to these expectation values, it would suffice to be able to detect singly occupied sites. Thus, the loss of atom pairs due to inelastic light-induced collision during the imaging process, showing up in the bosonic case, would not be an issue.
Measuring correlations of the spin-z component will be sufficient for obtaining the concurrence under the following assumptions: (i) the initial state is of the type we have described (total singlet, i.e. total spin 0); (ii) the dynamics proceeds according to the model Hamiltonian, i.e. a SU (2) symmetric Hamiltonian with no decoherence or entanglement with external degrees of freedom. A scenario where the initial state violates condition (i) would be a cluster with impurity sites that contain only single fermions. If there is only one impurity containing a single fermion, then the problem can still be diagnosed because the final snapshot will reveal unequal numbers of spin-up and spin-down fermions. In general, however, if these two conditions are in doubt, the experiment would have to measure the full spin-spin density matrix of two lattice sites, by repeated runs and measurements of spin projections in different directions. This could be done by implementing a rotation in spin space before measurements, in analogy to the state tomography realized with trapped ions [59]. While being more challenging than measuring the spin-dependent fermion number, such a kind of coherent spin control seems to be possible in future experiments. Spin flips at single lattice sites have already been shown experimentally for ultracold bosons [60]. This setup could be in principle extended to coherent spin control, replacing the Rabi frequency sweep by driving a Rabi oscillation [61].

Summary and outlook
In this paper, we have analysed the creation and time-evolution of spin-entanglement during the sudden expansion from a cluster of doublons into an empty lattice. Interestingly, remote spin-entangled pairs are created for large onsite interaction. In addition, an extended cluster favours the emission of the two fermions of a pair into the same direction when compared against the decay of a single doublon. Finally, we found that a time-dependent modulation of the tunnelling amplitude can be used to increase the 'production' of spin-entangled pairs. Our results provide a starting point for studying the spin-entanglement dynamics for more complex initial states, in different dimensionalities, e.g. the crossover from one to two dimensions, or for spin-imbalanced fermionic gases.
In the scenario considered here, the spin-entanglement can be extracted from the twosite spin-z correlation functions. Thus, it will become experimentally accessible once spindependent single-site detection has been implemented.
Here, J denotes the Bessel function of the first kind and we assumed an infinite lattice. For the previously discussed initial state of localized doublons, the Green function reads where the summation is taken over all initially occupied lattice sites O. For both spin species, the density is given by When the fermions are initially localized at the lattice sites, equal time, normal ordered products of operators can be expressed as a Slater determinant of the equal time one-particle Green functions, G i j (t), [62]. The correlation function Ŝ z i (t)Ŝ z j (t) can be evaluated by writing it in terms of density-density correlations and making use of Note that the probability of finding two fermions with the same spin at lattice sites i and j is by 2|G i j (t)| 2 smaller than the probability for fermions with antiparallel spin. This difference leads to a nonvanishing spin-spin correlation even in the absence of onsite interactions. For noninteracting fermions, the spin-spin correlation is related to the density-density correlation D i j (t) = n i (t)n j (t) − n i (t) n j (t) via Thus, the density-density correlation D i j (t) is smaller or equal to zero for different lattice sites i and j, cf figures 2(a) and (c)-(e). Analogously, the density-density correlation of single fermions is calculated using the relations for n i,a (t)n j,b (t) given above. We find The concurrence C i, j (t) (defined by equation (6)) can be evaluated using the expressions (A.5) and (A.9). It approaches one in the limit Ŝ z 2 can only assume values between 0 and 1/2).

Appendix B. Decay of a doublon into scattering states
This appendix presents the derivation of the decay probability of a fermionic doublon into different scattering states |ψ K ,k . In doing so, the wavefunctions of the scattering states are calculated mainly adopting the procedure for two-particle states in the Bose Hubbard model presented in [63].
For one spin-up and one spin-down fermion, and an infinite lattice the Hamiltonian (1) can be expressed in the form (B.1) Analogously, we write the two-fermion states in terms of the basis For the decay of a doublon, the fermions are always in the spin-singlet subspace. Thus, we can write the two-particle wavefunction as (↑ i , ↓ j ) = ϕ s (↑, ↓) · ψ(i, j), with the spin-singlet wavefunction ϕ s (a, b) = δ a,↑ δ b,↓ − δ a,↓ δ b,↑ and a symmetric spatial wavefunction, ψ(i, j) = ψ( j, i). Plugging | into the Schrödinger equation for Hamiltonian (B.1) yields This relation is simplified by introducing c.o.m. and relative coordinates, R = (i + j)/2 and r = i − j, respectively. The wavefunction factorizes into a plane wave motion of the c.o.m. with total wavenumber K ∈ [−π, π ) and a K -dependent relative motion, i.e. ψ(i, j) = e iK R ψ K (r ). The relative motion satisfies the recurrence relation with K -dependent tunnelling amplitude J K = 2J cos(K /2) and energy E K . For vanishing interaction strength U equation (B.3) is solved by plane waves ψ K ,k (r ) = e ±ikr with corresponding eigenenergies E K ,k = −2J K cos(k). In the interacting case, we make a scattering ansatz by writing the wavefunction as a superposition of incoming and reflected plane waves, ψ K ,k (r 0) = e ikr + c e −ikr . Here, k is the relative wavenumber and ψ K ,k (r < 0) is determined by the symmetry condition ψ(r ) = ψ(−r ). The boundary condition at r = 0 in equation (B.3) fixes the coefficient c and we obtain ψ K ,k (r ) = ψ K ,k (0) cos(kr ) + U 2J K sin(k) sin(k|r |) .

(B.4)
Finally, we compare the decay probability of a doublon for different scattering states |ψ K ,k . In doing so, we express |ψ K ,k (0)| 2 in terms of the average density in the relative coordinate,n, which is obtained by averaging |ψ K ,k (r )| 2 over one period 2π/k. Note thatn does only depend on the systems size and is independent of k, K , and U . We find that the decay probability into the scattering state |ψ K ,k equals |ψ K ,k (0)| 2 =n 1 + U 2 / 16J 2 cos 2 (K /2) sin 2 (k) −1 .
(B.5) lattice sites. It can be expressed as 1 − P D (t), with the doublon survival probability P D (t), i.e. the probability of finding the doublon intact at time t. In the long time limit, only fermions in a bound state remain localized close to each other. For a singlet pair with finite onsite interaction U = 0 and c.o.m. wavenumber K , it exists one bound state ψ b K (r ) (localized solution of equation (B.3), details are given below), where r is the relative coordinate. Thus, we obtain P D (∞) = 1 2π π −π dK |ψ b K (0)| 4 in the limit of an infinite lattice, where ψ b K (0) is the overlap between the doublon and the bound state. Analogously, we find C 1 (∞) = 1 2π π −π dK |ψ b K (0)| 2 |ψ b K (1)| 2 , where we have used that C 1 (t) = i n s i (t)n s i+1 (t) equals the probability of finding the fermions at time t at nearest-neighbour lattice sites. The bound state can be calculated using the exponential ansatz ψ b K (r ) = 1/ √ N K α |r | K for the wavefunction in equation (B.3). This gives α K = U/2J K − sign(U/2J K ) 1 + (U/2J K ) 2 and N K = (2J K /U ) 1 + (U/2J K ) 2 , with J K = 2 cos(K /2). Inserting this solution into the doublon survival probability yields P D (∞) = 1 2π π −π dK N −2 K = [1 + 16J 2 /U 2 ] −1/2 . In the limit of infinite times, the total summed concurrence and summed concurrence of nearest-neighbour lattice sites are, hence, given by In the long-time limit, the summed concurrences of lattice sites at same side (C ss (t)) and at different sides (C ds (t)) of initial doublon position can be related to the scattering states calculated in appendix B. We denote by C (scat) ds and C (scat) ss the contributions to the summed concurrences stemming from the scattering states. C (scat) ds is the sum of the occupation probabilities of scattering states corresponding to fermions moving in opposite direction [k 1 ∈ (0, π), k 2 ∈ (−π, 0) or k 1 ∈ (−π, 0), k 2 ∈ (0, π)], and C (scat) ss is the sum of the occupation probabilities of those states with fermions moving into the same direction [k 1 , k 2 ∈ (0, π) or k 1 , k 2 ∈ (−π, 0)]. Here, k 1 and k 2 denote the asymptotic wavenumbers of the two fermions. The occupation probabilities are the absolute square of the overlap between the doublon and the scattering state, |ψ K ,k (0)| 2 . The explicit form of |ψ K ,k (0)| 2 is given by equation (B.5). Taking the continuum limit we find (C.6) The summed concurrences C ss (t) and C ds (t) used for the numerical data (equations (9) and (10)) contain additional contributions from bound states. In the definition of C ss (t) we excluded nearest-neighbour lattice sites in order to remove most of these contributions (note that nearestneighbour sites do not appear in C ds (t)). From the expression for the bound state given above follows that all other terms are of the order O([J/U ] 4 ). Thus, C ss (t) and C ds (t) agree with C (scat) ss and C (scat) ds for large J t and U/J , cf figure 7(b). In conclusion, we find following relations between the summed concurrences for long evolution times and large interaction strengths U/J 1: C 1 = C tot /2, C ss = [1/4 − 2/π 2 ]C tot and C ds = [1/4 + 2/π 2 ]C tot .