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

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–5] and in the additional trapping potential by either switching it off [6] or displacing its centre [7–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–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–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–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 spin-changing 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).

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]

$$\begin{equation} \hat{\mathcal{H}}=-J\sum_{i=1}^{L-1}\sum_{a=\uparrow,\downarrow}\left\{ \hat{c}_{i,a}^{\dagger}\hat{c}_{i+1,a}+\hat{c}_{i+1,a}^{\dagger}\hat{c}_{i,a}\right\} +U\sum_{i=1}^{L}\hat{n}_{i,\uparrow}\hat{n}_{i,\downarrow}. \end{equation} \tag{ 1 }$$

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 is $\hat {n}_{i,a}=\hat {c}_{i,a}^{\dagger }\hat {c}_{i,a}$, with the creation (annihilation) operator $\hat {c}_{i,a}^{\dagger }$($\hat {c}_{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).

Zoom In Zoom Out Reset image size

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 = −2Jcos(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²/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].

For the numerical evaluation of the nonequilibrium time evolution, we employ the tDMRG method [42–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⁻⁵ or 10⁻⁶, 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.

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 $\hat {\rho }_{ij}(t)=\mbox {Tr}_{{\mathrm { sites}}\neq i,j}\{|\Psi (t)\rangle \langle \Psi (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 $\hat {\rho }_{ij}(t)$ onto those states yields

$$\begin{equation} \hat{\rho}_{ij}^{\mathrm{s}}(t)=\frac{1}{\mbox{Tr}\{\hat{\rho}_{ij}(t)\,\hat{n}_{i}^{\mathrm{s}}\hat{n}_{j}^{\mathrm{s}}\}}\hat{n}_{j}^{\mathrm{s}}\hat{n}_{i}^{\mathrm{s}}\,\hat{\rho}_{ij}(t)\,\hat{n}_{i}^{\mathrm{s}}\hat{n}_{j}^{\mathrm{s}}. \end{equation} \tag{ 3 }$$

Here, $\hat {n}_{i}^{\mathrm {s}}=\hat {n}_{i,\uparrow }+\hat {n}_{i,\downarrow }-2\hat {n}_{i,\uparrow }\hat {n}_{i,\downarrow }$ 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, $\mbox {Tr}\{\hat {\rho }_{ij}(t)\hat {n}_{i}^{\mathrm {s}}\hat {n}_{j}^{\mathrm {s}}\}=\langle \hat {n}_{i}^{\mathrm {s}}(t)\hat {n}_{j}^{\mathrm {s}}(t)\rangle$, is the probability of finding at time t a single fermion at each of the lattice sites i and j. The state described by $\hat {\rho }_{ij}^{\mathrm {s}}(t)$ is the state obtained after a successful projective measurement. The reduced density matrix $\hat {\rho }_{ij}^{\mathrm {s}}(t)$ is equivalent to a two-qubit density matrix, which can be expressed in the form $\hat {\rho }=\sum _{\alpha ,\beta =0}^{3}\lambda ^{\alpha \beta }\,\hat {\sigma }_{(1)}^{\alpha }\otimes \hat {\sigma }_{(2)}^{\beta }$, where $\hat {\sigma }^{1,2,3}$ are the Pauli matrices, $\hat {\sigma }^{0}=\mathbf {1}$, and the factors λ^αβ are determined by the correlation functions $\lambda ^{\alpha \beta }=\frac {1}{4}\langle \hat {\sigma }_{(1)}^{\alpha }\hat {\sigma }_{(2)}^{\beta }\rangle$ [52]. Consequently, the reduced density matrix of single fermions at lattice sites i and j can be written as

$$\begin{equation} \hat{\rho}_{ij}^{\mathrm{s}}(t)=\frac{1}{\langle\hat{n}_{i}^{\mathrm{s}}(t)\hat{n}_{j}^{\mathrm{s}}(t)\rangle}\sum_{\alpha,\beta=0}^{3}\langle\hat{S}_{i}^{\alpha}(t)\hat{S}_{j}^{\beta}(t)\rangle\,\hat{\sigma}_{i}^{\alpha}\otimes\hat{\sigma}_{j}^{\beta}, \end{equation} \tag{ 4 }$$

where $\hat {S}_{i}^{1,2,3}=\frac {1}{2}\sum _{a,b=\uparrow ,\downarrow }\hat {c}_{i,a}^{\dagger }(\hat {\sigma }^{1,2,3})_{\ b}^{a}\hat {c}_{i,b}$ is the x-, y-, and z-component of the spin operator and $\hat {S}_{i}^{0}$ is, for compactness, defined as half the single fermion number operator, $\hat {S}_{i}^{0}:=\frac {1}{2}\hat {n}_{i}^{\mathrm {s}}$. Note that $\langle \hat {S}_{i}^{\alpha }(t)\hat {S}_{j}^{\beta }(t)\rangle$ 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 $\hat {\rho }_{ij}^{\mathrm {s}}(t)$. As detailed now, the reduced density matrix $\hat {\rho }_{ij}^{\mathrm {s}}(t)$ depends only on $\langle \hat {S}_{i}^{z}(t)\hat {S}_{j}^{z}(t)\rangle /\langle \hat {n}_{i}^{\mathrm {s}}(t)\hat {n}_{j}^{\mathrm {s}}(t)\rangle$ for the cluster initial state shown in figure 1. The Hubbard Hamiltonian (1) preserves the spin-dependent particle number, i.e. $[\hat {\mathcal {H}},\hat {N}_{\uparrow ,\downarrow }]=0$, $\hat {N}_{\uparrow ,\downarrow }=\sum _{i=1}^{L}\hat {n}_{\uparrow ,\downarrow }$. 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 $\langle \hat {n}_{i}^{\mathrm {s}}(t)\hat {S}_{j}^{x,y}(t)\rangle =0$, $\langle \hat {S}_{i}^{x,y}(t)\hat {S}_{j}^{z}(t)\rangle =0$ and $\langle \hat {S}_{i}^{x}(t)\hat {S}_{j}^{x}(t)\rangle -\langle \hat {S}_{i}^{y}(t)\hat {S}_{j}^{y}(t)\rangle \pm \mathrm {i}[\langle \hat {S}_{i}^{x}(t)\hat {S}_{j}^{y}(t)\rangle +\langle \hat {S}_{i}^{y}(t)\hat {S}_{j}^{x}(t)\rangle ]=0$. 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 $\langle \hat {S}_{i}^{x}(t)\hat {S}_{j}^{x}(t)\rangle =\langle \hat {S}_{i}^{y}(t)\hat {S}_{j}^{y}(t)\rangle$ and $\langle \hat {S}_{i}^{x}(t)\hat {S}_{j}^{y}(t)\rangle =-\langle \hat {S}_{i}^{y}(t)\hat {S}_{j}^{x}(t)\rangle$. Moreover, the Hamiltonian (1) is fully rotationally invariant, $[\hat {\mathcal {H}},\sum _{i=1}^{L}\hat {S}_{i}^{x,y,z}]=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 follows that $\langle \hat {S}_{i}^{z}(t)\hat {S}_{j}^{z}(t)\rangle =\langle \hat {S}_{i}^{x}(t)\hat {S}_{j}^{x}(t)\rangle =\langle \hat {S}_{i}^{y}(t)\hat {S}_{j}^{y}(t)\rangle$, $\langle \hat {n}_{i}^{\mathrm {s}}(t)\hat {S}_{j}^{z}(t)\rangle =\langle \hat {n}_{i}^{\mathrm {s}}(t)\hat {S}_{j}^{x,y}(t)\rangle =0$, and $\langle \hat {S}_{i}^{x}(t)\hat {S}_{j}^{y}(t)\rangle =\langle \hat {S}_{i}^{x,y}(t)\hat {S}_{j}^{z}(t)\rangle =0$. In summary, the reduced density matrix for the expansion from a cluster of doublons reads

$$\begin{eqnarray} \hat{\rho}_{ij}^{\mathrm{s}}(t) & = & \frac{1}{4}\cdot\mathbf{1}+\frac{\langle\hat{S}_{i}^{z}(t)\hat{S}_{j}^{z}(t)\rangle}{\langle\hat{n}_{i}^{\mathrm{s}}(t)\hat{n}_{j}^{\mathrm{s}}(t)\rangle}\left[\hat{\sigma}_{i}^{x}\otimes\hat{\sigma}_{j}^{x}+\hat{\sigma}_{i}^{y}\otimes\hat{\sigma}_{j}^{y}+\hat{\sigma}_{i}^{z}\otimes\hat{\sigma}_{j}^{z}\right]\nonumber \\ & = &\left(\frac{1}{4}+\frac{\langle\hat{S}_{i}^{z}(t)\hat{S}_{j}^{z}(t)\rangle}{\langle\hat{n}_{i}^{\mathrm{s}}(t)\hat{n}_{j}^{\mathrm{s}}(t)\rangle}\right)\left[|T_{ij}^{1}\rangle\langle T_{ij}^{1}|+|T_{ij}^{0}\rangle\langle T_{ij}^{0}|+|T_{ij}^{-1}\rangle\langle T_{ij}^{-1}|\right]\nonumber \\ & & +\left(\frac{1}{4}-3\frac{\langle\hat{S}_{i}^{z}(t)\hat{S}_{j}^{z}(t)\rangle}{\langle\hat{n}_{i}^{\mathrm{s}}(t)\hat{n}_{j}^{\mathrm{s}}(t)\rangle}\right)|S_{ij}\rangle\langle S_{ij}|, \end{eqnarray} \tag{ 5 }$$

where $|S_{ij}\rangle =\frac {1}{\sqrt {2}}(|\uparrow _{i}\downarrow _{j}\rangle -|\downarrow _{i}\uparrow _{j}\rangle )$ is the singlet state and |T^m_ij〉 is the triplet state with the m denoting the spin projection in z-direction.

Instead of $\langle \hat {S}_{i}^{z}(t)\hat {S}_{j}^{z}(t)\rangle$ one could in principle evaluate the spin–spin correlation in any other direction. However, for cold atomic gases in optical lattices the correlation function $\langle \hat {S}_{i}^{z}(t)\hat {S}_{j}^{z}(t)\rangle =\frac {1}{4}\langle [\hat {n}_{i,\uparrow }(t)-\hat {n}_{i,\downarrow }(t)][\hat {n}_{j,\uparrow }(t)-\hat {n}_{j,\downarrow }(t)]\rangle$ 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(\hat {\rho })$ [50] to quantify the entanglement. Given the time-reversed density matrix $\hat {\tilde {\rho }}=\sigma ^{y}\otimes \sigma ^{y}\hat {\rho }^{*}\sigma ^{y}\otimes \sigma ^{y}$, with the complex conjugation $\hat {\rho }^{*}$ taken in the standard basis {|↑↑〉,|↑↓〉,|↓↑〉,|↓↓〉}, the concurrence is defined by $C(\hat {\rho })=\mbox {max}\{0,\sqrt {\lambda _{1}}-\sqrt {\lambda _{2}}-\sqrt {\lambda _{3}}-\sqrt {\lambda _{4}}\}$, where the λ_i are the eigenvalues of $\hat {\rho }\hat {\tilde {\rho }}$ in descending order. In our case, the concurrence of the reduced density matrix (5) is given by

$$\begin{equation} C_{i,j}(t)=\mbox{max}\left\{ 0,-\frac{1}{2}-6\frac{\langle\hat{S}_{i}^{z}(t)\hat{S}_{j}^{z}(t)\rangle}{\langle\hat{n}_{i}^{\mathrm{s}}(t)\hat{n}_{j}^{\mathrm{s}}(t)\rangle}\right\}. \end{equation} \tag{ 6 }$$

Spin–spin correlations $\langle \hat {S}_{i}^{z}(t)\hat {S}_{j}^{z}(t)\rangle /\langle \hat {n}_{i}^{\mathrm {s}}(t)\hat {n}_{j}^{\mathrm {s}}(t)\rangle \geqslant -1/12$ result in a vanishing concurrence. The concurrence approaches 1 as $\langle \hat {S}_{i}^{z}(t)\hat {S}_{j}^{z}(t)\rangle /\langle \hat {n}_{i}^{\mathrm {s}}(t)\hat {n}_{j}^{\mathrm {s}}(t)\rangle \searrow -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):

$$\begin{equation} P_{ij}^{\rm Singlet}(t)=\mbox{Tr}\left[\hat{\rho}_{ij}^{\mathrm{s}}(t)|\mathcal{S}_{ij}\rangle\langle\mathcal{S}_{ij}|\right]=\frac{1}{4}-3\frac{\langle\hat{S}_{i}^{z}(t)\hat{S}_{j}^{z}(t)\rangle}{\langle\hat{n}_{i}^{\mathrm{s}}(t)\hat{n}_{j}^{\mathrm{s}}(t)\rangle}. \end{equation} \tag{ 7 }$$

Each of the spin triplet states is measured with the probability $1/4+\langle \hat {S}_{i}^{z}(t)\hat {S}_{j}^{z}(t)\rangle /\langle \hat {n}_{i}^{\mathrm {s}}(t)\hat {n}_{j}^{\mathrm {s}}(t)\rangle$.

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 $\langle \hat {n}_{i}^{\mathrm {s}}(t)\hat {n}_{j}^{\mathrm {s}}(t)\rangle$, with the single fermion number operator already defined above ($\hat {n}_{i}^{\mathrm {s}}=\hat {n}_{i,\uparrow }+\hat {n}_{i,\downarrow }-2\hat {n}_{i,\uparrow }\hat {n}_{i,\downarrow }$). Figure 3 shows the numerical results for $\langle \hat {n}_{i}^{\mathrm {s}}(t)\hat {n}_{j\neq i}^{\mathrm {s}}(t)\rangle$ for different evolution times Jt and onsite interaction strengths U/J.

Zoom In Zoom Out Reset image size

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 $\langle \hat {n}_{i}^{\mathrm {s}}(t)\hat {n}_{j}^{\mathrm {s}}(t)\rangle$ 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, $\sum _{i=1}^{L}\langle \hat {n}_{i}^{\mathrm {s}}\rangle$, approaches within a few Jt an almost constant value, cf area under the curves for Jt = 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, $\langle \hat {n}_{i}^{\mathrm {s}}(t)\rangle$ becomes relatively flat in the centre of the cloud for larger evolution times Jt.

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 Jt and onsite interaction strengths U/J.

Zoom In Zoom Out Reset image size

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 spin-entanglement 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 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}\propto \sum _{lm}\mathrm {e}^{-\mathrm {i}k(l-m)}\langle \hat {S}_{l}^{z}\hat {S}_{m}^{z}\rangle$. 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.

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, $\langle \hat {n}_{i}^{\mathrm {s}}(t)\hat {n}_{j}^{\mathrm {s}}(t)\rangle$. 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–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

$$\begin{equation} C_{{\rm tot},i}(t)=\sum_{j\neq i}\langle\hat{n}_{i}^{\mathrm{s}}(t)\hat{n}_{j}^{\mathrm{s}}(t)\rangle\, C_{i,j}(t). \end{equation} \tag{ 8 }$$

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 $\langle \hat {n}_{i}^{\mathrm {s}}\rangle$. 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.

Zoom In Zoom Out Reset image size

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

$$\begin{equation} C_{\rm ss}(t) = \sum_{i+1<j<l\:\vee\: i-1>j>r}\langle\hat{n}_{i}^{\mathrm{s}}(t)\hat{n}_{j}^{\mathrm{s}}(t)\rangle\, C_{i,j}(t), \end{equation} \tag{ 9 }$$

$$\begin{equation} C_{\rm ds}(t) = \sum_{i<l\:\wedge\: j>r}\langle\hat{n}_{i}^{\mathrm{s}}(t)\hat{n}_{j}^{\mathrm{s}}(t)\rangle\, C_{i,j}(t), \end{equation} \tag{ 10 }$$

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, $C_{d}(t)=\sum _{i}\langle \hat {n}_{i}^{\mathrm {s}}(t)\hat {n}_{i+d}^{\mathrm {s}}(t)\rangle C_{i,i+d}(t)$, and the summed concurrence of all sites, $C_{{\mathrm { tot}}}(t)=\sum _{i<j}\langle \hat {n}_{i}^{\mathrm {s}}(t)\hat {n}_{j}^{\mathrm {s}}(t)\rangle C_{i,j}(t)$. The time evolution of these summed concurrences is shown in figure 6.

Zoom In Zoom Out Reset image size

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 spin-entanglement 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₈(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 Jt ≳ 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 Jt 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 ⩽ 4Jt and roughly decays as (Jt)⁻¹ for the simulated times. Note that C₁(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 (Jt ≈ 20) for the comparison of the summed concurrences are reached only for U/J ≳ 6. The summed concurrence of all sites, C_tot, 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.

Zoom In Zoom Out Reset image size

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 amount of spin-entanglement as shown in figure 8. For the presented case, the repeated quenches lead to the generation of more and more single fermions, see figure 8(a). This results in an enhanced spin-entanglement between distant lattice sites, while the spin-entanglement between nearest-neighbour sites is suppressed (compare figures 8(b) and (c) with figures 6(c) and (f)).

Zoom In Zoom Out Reset image size