Engineering entanglement Hamiltonians with strongly interacting cold atoms in optical traps

We present a proposal for the realization of entanglement Hamiltonians in one-dimensional critical spin systems with strongly interacting cold atoms. Our approach is based on the notion that the entanglement spectrum of such systems can be realized with a physical Hamiltonian containing a set of position-dependent couplings. We focus on reproducing the universal ratios of the entanglement spectrum for systems in two different geometries: a harmonic trap, which corresponds to a partition embedded in an infinite system, and a linear potential, which reproduces the properties of a half-partition with open boundary conditions. Our results demonstrate the possibility of measuring the entanglement spectra of the Heisenberg and XX models in a realistic cold-atom experimental setting by simply using gravity and standard trapping techniques.

Introduction. The study of entanglement in quantum many-body systems [1,2] has become one of the major efforts in the physics community, not only because it is a central feature of quantum theories, but also due to its potential for describing quantum phases of matter and topological order [3][4][5][6][7][8][9]. Directly measuring entanglement in experiments, on the other hand, has proven to be a challenging task. This is due, in particular, to the difficulty of obtaining the full-state tomography [10] of many-body systems. Nevertheless, outstanding progress has been made in recent years with respect to the extraction of the entanglement properties of quantum systems, both in terms of theoretical proposals [11][12][13][14][15][16][17] and experiments [18][19][20][21]. In common, these works have the feature of employing indirect measurement protocols, either through the probing of correlations or by interference of identical copies of the system.
In this context, it is highly desirable to have at hand alternative proposals [22] for the measurement of entanglement that are at the same time direct -in terms of quantities which are ordinarily accessible in experiments -and scalable. A significant step in this direction has been taken recently [23][24][25][26], with works showing that the entanglement spectrum [27,28] of lattice Hamiltonians can be reproduced by obtaining the physical spectrum of a Hamiltonian with the same general properties, but with a set of spatially varying coupling parameters. This notion is based on the Bisognano-Wichmann (BW) theorem, which originally describes entanglement Hamiltonians for continuous systems in the context of quantum field theory [29,30]. Later developments have shown [31] that particular features of the original model, such as boundary conditions, have a direct effect on the functional form of the couplings in the physical Hamiltonian. Such a relation hints at the prospect of reproducing the * barfknecht@lens.unifi.it entanglement properties of discrete systems with Hamiltonians with specifically designed non-homogeneous couplings.
An immediate possibility that arises in this context is the simulation of entanglement Hamiltonians with trapped systems of cold atoms [32], where optical confinement, atomic interactions and internal states can be manipulated with remarkable precision. Such setups are therefore ideal candidates for quantum simulations of condensed matter systems. [33][34][35][36]. On the theoretical side, it has been demonstrated that if the interactions between atoms are strong enough, the system can be mapped from the continuum to a spin chain, where nearest-neighbor couplings are determined by the local trapping geometry [37][38][39][40]. Experimentally, the strongly correlated regime is accessible with cold atoms both for fermionic [41][42][43][44][45] and bosonic [46][47][48][49][50] gases. In fact, the validity of the theoretical approach described above has been verified in recent experiments [51].
Inspired by this perspective, in the present work we show how to engineer the entanglement Hamiltonians of spin systems, such as the Heisenberg (XXX) and XX models, by considering a strongly interacting two-component system of cold atoms in effectively onedimensional optical traps. Our first application is the case of harmonically trapped atoms, which gives rise to a spin chain where the couplings approximately follow a parabolic distribution. As we will show, this remarkably ordinary assumption is enough to reproduce the universal ratios of the entanglement spectrum of a partition embedded in an infinite system with periodic boundary conditions. The second application assumes the presence of a linear potential, which in turn results in a set of linearly increasing couplings for the spin chain. In this regime, the physical Hamiltonian can reproduce the ratios of the entanglement spectrum of a half-partition of a spin system with open boundary conditions. The simple linear form of the potential additionally poses the intriguing possibility of using gravity as a tool for probing entanglement. To account for other experimental-related aspects, we also include Density Matrix Renormalization Group (DMRG) simulations of the continuum away from the strongly interacting limit, as well as the expected effect of finite temperatures on measurable quantities such as the dynamical structure factor.
System description We calculate the entanglement properties of one-dimensional spin chains such as the Heisenberg model, , where σ denotes the Pauli vector and J is a homogeneous coupling. After finding the ground state solution |Ψ for one of these models, the reduced density matrix for a subsystem A can be calculated by tracing over the remaining subsystem B with ρ A = Tr B |Ψ Ψ|= e −H A /Z A , where H A is called the entanglement Hamiltonian and Z A is a normalization constant. The entanglement spectrum for partition A is then obtained as the set of eigenvalues { n } = − ln{ A n }, where A n denotes the eigenvalues of the reduced density matrix.
The BW theorem for discrete systems [24] states that the entanglement Hamiltonian H A can be equivalently calculated as where J i is a set of position-dependent couplings. Remarkably, predictions from conformal field theories predict [31,52] that these couplings should be given by for entanglement Hamiltonians corresponding to partitions embedded in an infinite system with periodic boundary conditions (which we label T 1 ) and halfpartitions in systems with open boundary conditions (T 2 ), respectively (see Fig. 1). These results indicate that the entanglement Hamiltonians of spin systems can be simulated by constructing a physical Hamiltonian with a set of properly engineered couplings. A suitable platform for such an endeavor is an effectively one-dimensional cold atomic gas with a Hamiltonian described by where the first term on the right-hand side includes the presence of a trapping potential V (x) and the second term accounts for contact interactions between atoms in different internal states (we assume a two-component gas where the internal states are labeled as ↑, ↓). The strength of the interactions is set by the parameter g (in units ofh 2 /ml where l is the characteristic length of the system), and we initially consider the fermionic case where interactions between atoms in the same internal state are forbidden by the Pauli principle (we will relax FIG. 1. The entanglement properties of one-dimensional spin chains can be simulated with strongly interacting systems of cold atoms in optical tubes. These systems can be mapped to spin chains where the couplings are determined by the local trapping geometry. a) A partition embedded in an infinite spin system with periodic boundary conditions (T1), can be simulated with b) a strongly interacting (g 1) twocomponent system of cold atoms in a harmonic trap. We model this system as a spin chain where the exchange coefficients αi follow the distribution of an inverted parabola. c) A half-partition of a finite spin chain with open boundary conditions is reproduced by d) cold atoms in a linear potential, which can be realized by gravity. This system is mapped into a spin chain where the coefficients increase linearly with the position.
this restriction later when discussing the simulations of the XX model).
In the limit of strong repulsion (g 1), the wave function of a system described by Hamiltonian is the wave function in the limit of infinite repulsion and each term of the sum includes a permutation P k (with amplitude a k ) of the coordinates. By employing the Hellmann-Feynman theorem dE dg = Ψ| dH dg |Ψ along with boundary conditions for the (3) can be mapped into the following spin chain [38,53,54] (see Supplemental Materials [55] for details) where E 0 denotes the energy of the system at the fermionization limit where g = ∞. This mapping can be interpreted as a perturbation with respect to the limit of infinite repulsion (notice that we can rewrite this Hamiltonian in terms of the permutation operator [38]). By diagonalizing this Hamiltonian we find the amplitudes a k for the wave function Ψ and the energy spectrum in the limit of strong interactions. Another fundamental aspect of the approach described above is that, in Eq. (4), the positiondependent exchange coefficients α i are obtained from the properties of the spatial wave function Φ 0 using the fol-lowing relation which is independent of spin. The many-body wave function described by Φ 0 is constructed as the Slater determinant of the N lowest single-particle states in the trapping potential V (x). Aside from the symmetry considerations regarding permutations of particles, this wave function is the same as the one that describes a Bose gas in the limit of infinite repulsion [56]. In Eq. (4), the energy E 0 is thus simply given by the sum of the energies of each of these single-particle states. It is therefore clear that once we determine the geometry of the trap V (x), we can obtain the exchange coefficients that describe it in Eq. (4). While solving the integrals given by Eq. (5) can be cumbersome for large N , efficient methods which exploit the determinant properties of Φ 0 are available [57,58]. In Fig. 1 we show a sketch of the protocol adopted as a proposal for the simulation of entanglement with trapped systems of cold atoms. We have thus established our proposal for simulating the entanglement spectrum of a spin chain with a twocomponent atomic system, where the energy levels can be obtained by standard spectroscopy techniques [59]. The presence of strong interactions and external trapping potentials not only allows for benchmarking the results against the limit of infinite repulsion, but also automatically generates the set of coupling needed for engineering entanglement Hamiltonians. The relation between entanglement and inhomogeneities in the underlying geometry of the system has been also explored, particularly for non-interacting fermions, in [60][61][62][63].
Results. We now calculate the properties of the atomic systems described above, including the particle distributions and numerical values of the exchange coefficients for a given trapping potential. In Fig. 2, we show the single-particle densities for the wave function Φ 0 for spinless fermions in a) the harmonic trap V (x) = x 2 /2 and in b) the case of a finite system confined by hard walls at x = 0 and x = l and exposed to a linear potential V (x) = V 0 (l − x). Here we observe how the overlaps between neighboring particles are affected by the underlying geometry. This feature is reflected in the numerical values of the exchange coefficients α i , which are shown in Fig. 2 c) and d). The distribution of these values, for the harmonic trap, is symmetric across the origin (therefore it is enough to calculate at most the first N/2 coefficients) and has the shape of an inverted parabola [64,65]. For the tilted potential, we find that the coefficients increase linearly and are simply given by α i = V 0 i. In both Figs. 2 c) and d), we additionally show a comparison with the results obtained by fitting the functions in Eq. (2), for two different sizes of partitions. We find excellent agreement between the set of couplings given by Eqs. (2) and those generated by the two choices of trapping potential.
We are now able to calculate the exact entanglement spectrum for these different types of partitions and com- pare to results obtained with a physical Hamiltonian with spatially varying couplings. For the Heisenberg model, the physical Hamiltonian is naturally given by Eq. (4). The XX model can analogously by reproduced by initially considering a bosonic system in the continuum, where interaction between atoms in the same internal states are allowed. This results in an additional term in Eq. (4) given by − 1 . For a case of imbalanced interactions (κ = 2) this bosonic system maps into the XX model [53].
Our focus is on obtaining universal ratios for the entanglement spectrum, defined as where n denotes the energy level, 0 is the ground state energy and r is a reference energy level (unless stated otherwise, we fix r = 3 -the second excited state). In Fig. 3 a)-d) we show the results for this quantity obtained by exactly diagonalizing different systems. The physical spectrum of the harmonically trapped strongly interacting Hamiltonian is compared to a partition embedded in a system with periodic boundary conditions, while the case of linear potential is compared to a half-partition of the same size in a system with open boundary conditions. Additionally, we include the results for an "ideal" physical Hamiltonian with a set of couplings provided by Eq.
(2). The agreement between the results obtained with three different approaches is particularly remarkable for small values of n, corresponding to the low-energy part of the entanglement spectrum. The discrepancy with the exact results obtained from the reduced density matrix (black solid curves) at higher energies stems mainly from finite-size effects. In Fig. 3 e)-f), we include also results obtained with DMRG, where we see that the agreement between the entanglement spectrum of homogeneous partitions of the Heisenberg model and the physical spectrum of trapped fermions also holds for larger systems.
Experimental details. Effectively one-dimensional systems of cold atoms are currently realized in the lab by loading the atoms into tight optical waveguides, where the confinement along the transverse direction is much larger than the longitudinal one. These waveguides can be provided by optical lattices, where the transverse confinement could easily reach values on the order of ω ⊥ = 2π50 kHz, withhω ⊥ much larger than temperature and atomic chemical potential.
Different spin states can be simulated by exploring the internal degrees of freedom of particular atomic species. Usually, these are hyperfine states in alkali elements (such as fermionic 6 Li and bosonic 87 Rb), but nuclear-spin states in alkaline-earth elements like 173 Yb or 87 Sr can also be used. These last cases also provide the opportunity of exploring additional internal states with SU(N ) symmetry) [44,66].
With state-of-the-art atom trapping techniques there are ample possibilities of controlling the potential parameters to engineer the required exchange coupling α i . Typical axial harmonic trapping frequencies for the protocol shown in Fig. 1 b) lie in the range ω T = 2π(10 − 10 3 ) Hz. The protocol shown in Fig. 1 d) can be realized either by means of magnetic/optical gradients or by exploiting the effect of gravity, where the linear potential could be tuned by changing the tilting angle of the tubes. An additional focused laser beam (or other spatial-lightmodulator-based techniques) is required to create hard walls at the bottom of the tubes for keeping the atoms confined.
We also address two other important experimental aspects which can be relevant in the detection of the physical spectrum of trapped systems: finiteness of interactions and temperature. To investigate the first effect, we realize simulations of the continuum with the Hubbard model in an underlying harmonic trap [55]. In Fig. 4  a) we show a comparison of the universal ratios obtained with this approach to the expected results for a Heisenberg model with a set of coefficients given by Eq. (2). We find that even at a moderate interaction regime (well within the experimental capabilities) the measured spectrum agrees with the entanglement spectrum of the target model (especially in the low-energy sector at small n). We also note that this calculation provides an independent check of the validity of the spin-mapping approach described in the main part of this letter. In experimental setups, interactions between atoms can be manipulated by means of Feshbach [67] or confinement induced resonances [68].
We quantify the effect of finite temperature on the energy spectrum by calculating the temperature-dependent dynamical structure factor [69][70][71] where Z(T ) = Tr(e −H/k B T ) is the canonical partition function, k B is the Boltzmann constant and |i and |j denote the eigenstates; S z q = i (2) sin(qi)S z i / (N + 1) is the Fourier transform of the operator S z i = σ z i /2 and q = nπ/(N + 1) are discrete momentum values. In Fig. 4  b) we compare the excitation spectrum of the dynamical structure factor (at values of k B T corresponding to different fractions of the Fermi energy F ) to the position of the energy gaps corresponding to the universal ratios. Particularly, at lower temperatures we find pronounced peaks located precisely at the energy values predicted by Eq. (6), that are clearly visible for experimentallyachievable temperatures 0.05 F (note the vertical logarithmic scale). As expected, for larger temperatures such FIG. 4. a) Comparison of the universal ratios of the Heisenberg model with the ideal coefficients in Eq. (2) (gray dashed curves) to DMRG simulations of the continuum for fermions in a harmonic trap. Blue, yellow and red circles correspond to interaction strength g = 5, 10 and 15 (in units ofh 2 /ml), respectively. The universal ratios κn are calculated having as a reference state r = 2. b) Temperature-dependent dynamical structure factor (summed over q) for harmonically trapped fermions (Eq. (4)) with g = 25. Blue, yellow and red curves correspond kBT = 0.002 F , 0.05 F and 0.2 F , respectively, where F is the system's Fermi energy. The vertical gray dashed lines denote the position of the energy gaps corresponding to the universal ratios obtained for the Heisenberg model with couplings given by Eq. (2). The reference states for these calculations is r = 4, and the frequency ω0 is analogously defined as ω0 = ( 4 − 0)/h. In both panels, we assume N = 7 in a sector of fixed magnetization +1/2. results are washed out by the contribution of several additional frequencies. We point out that not all excitation peaks present in Fig. 4 b) are contemplated with a given definition of the universal ratios given by Eq. (6).
Concluding remarks. We have presented a theoretical proposal for the realization of entanglement Hamiltonians with strongly interacting cold atoms in effectively one-dimensional optical traps. A key feature of these models is the possibility of mapping the Hamiltonian into a spin chain with a set of couplings that depend on the underlying geometry. Particularly, by simply assuming a harmonic confinement, we have shown that the universal ratios of the entanglement for a partition embedded in an infinite system can be reproduced. Analogously, a linear potential can produce the results expected for a halfpartition in a finite system. The energy spectrum of these systems can thus be obtained by standard spectroscopy of cold atoms in elongated harmonic traps or box-like traps under the effect of gravity. We have benchmarked the robustness of our predictions against important experimental effects such as finite interaction strength and finite temperature, evidencing the experimental feasibility of our proposal. Such a protocol can be extended to the study of various spin Hamiltonians, like systems with higher internal symmetries [72]; in such cases, the particular details (e.g. interactions and internal states) of a given atomic model will determine the structure of the spin chain [73,74]. Nevertheless, the relation between the couplings between nearest neighbors and the geometry of the external trap remains valid in the limit of strong interactions.

SUPPLEMENTAL MATERIAL
Mapping a strongly interacting atomic system to a spin chain The general procedure described in this section has been developed and extensively detailed in different works, such as [37][38][39]. Further discussions and applications to different atomic models can be found, for instance, in [40,73]. For clarity, we reproduce below the essential steps required to obtain the spin chain models considered in the main text. We start by considering a one-dimensional two-component fermionic system with contact interaction, described by where N = N ↑ + N ↓ . In the limit of infinite repulsion (g → ∞), we can write the complete many-body wave function as where P k is the permutation operator and the sum is carried over a total number of L(N ↑ , N ↓ ) = N ↑ +N ↓ N ↑ permutations. In this context Φ 0 is the wave function in the limit of infinite repulsion where {x ↑i , x ↓j } denotes a given ordering of the particles.
For strong finite interactions, the Hellmann-Feynman theorem can be employed to write Additionally, we write the expression for the derivative condition at the contact point between particles as which follows the guidelines of the coordinate Bethe ansatz approach [76]. Plugging Eq. (11) into Eq. (10) and integrating with respect to g, we find where we discard terms of order O(1/g 2 ) and higher. In this expression, E 0 denotes the energy of a system of spinless fermions (which is the same energy expected in the regime of infinite repulsion). By inserting (9) in this equation, we find where a ik denotes the wave function coefficient where neighboring ↑ and ↓ fermions are found in position i and i + 1, while a ik is the coefficient for the wave function where the particle positions are i + 1 and i. In this expression, we have which is not spin-dependent. Thus, Φ 0 can now be treated as the wave function for a system of spinless fermions. We now take under consideration a spin Hamiltonian given by where Π i,i+1 ) is the permutation operator. A general wave function for this Hamiltonian can be written as in analogy with Eq. (9). Using this wave function, we calculate χ|H|χ as where a ik and a ik are as described for Eq. (13). This means that expressions 13 and 17 are equivalent provided that J i = α i /g. We can now rewrite Eq. 15 as which is the Hamiltonain given in Eq. (4).

Details on the simulations
In Fig. 2 a) and b) of the main text, we calculate the spatial distributions for a system of atoms in a harmonic trap and in a linear potential. This quantity is calculated with the following expression: where the integration is restricted to the sector Γ = x 1 < x 2 < ... < x N . For larger systems, it is convenient to explore the determinant properties of the wave function Φ 0 [77], and write where the elements of the matrix B(x) are written as b mn (x) = x −∞ dy ϕ m (y)ϕ n (y), and ϕ(x) denotes the singleparticle states in a corresponding trapping potential. For the simple case of the harmonic trap, these can be obtained exactly; for the linear potential, we obtain the single particle solutions by numerical diagonalization, using as a basis 50 eigenstates of the box potential. The characteristic length l in the harmonic trap is related to the trapping frequency ω by l = h/mω (we assume ω = 1 in our calculations). The exchange coefficients shown in Fig. 2 c) and d) can be calculated for small systems with Eq. (5). For larger system, we use the open-source program CONAN [57]. Most results found in Fig. 3 are obtained by considering the full energy spectrum calculated by exactly diagonalizing the corresponding Hamiltonians. In c) and d), the black curves are calculated with DMRG simulations of the Heisenberg model, with a total of 6 DMRG sweeps and a maximum bond dimension of 600.
In Fig. 4, we obtain the results at intermediate interactions in a system of N = 7 fermions with DMRG by approximating the continuum with the Hubbard model H = −t j,σ (c † j+1,σ c j,σ + H.c.) + U j n j,↑ n j,↓ + j,σ V j n j,σ , where the last term denotes the underlying harmonic trap potential. The simulation of the continuum is performed with a total of N S = 200 sites. By fixing a length λ we define the lattice spacing as a = λ/N s . The hopping parameter is then related to the kinetic term in the continuum as t = 1/(2a 2 ) (assuming m = 1), while the interaction parameters are related by U = g/a [78]. In these simulations we perform a total of 40 DMRG sweeps, with a maximum bond dimension of 10 6 .

Additional results
In the main text we showed, in Fig. 3 e)-f), the comparison between the universal ratios obtained with ED for trapped atoms and DMRG for partitions embedded in larger systems, focusing on the Heisenberg model. Here we present in Fig. 5 also the case of the XX model, where the size of the system is the same as considered in the Fig.  3 e)-f). As in the main text, we observe an increasing discrepancy for higher energy levels, which can be partially explained by a lack of convergence of the DMRG runs for the choice of parameters described above. This results can be therefore improved by considering additional runs with larger bond dimensions.
In the results contained in the main text for the temperature-dependent dynamical structure factor, we approximate the delta function δ [ω − (E i − E j )] contained in Eq. (7) with a Lorentzian given by f (ω) = 1 π η 2 η 2 + ω 2 (22) where η = 0.002. In Fig. 4 b) we summed over all values of momentum q. In Fig. 6 we also show the separate results for this quantity at particular momentum values.