Quantum simulation of long range $XY$ quantum spin glass with strong area-law violation using trapped ions

TL;DR

This paper demonstrates that a one-dimensional long-range disordered quantum spin glass exhibits strong area-law violation in its ground state, and proposes a trapped ion quantum simulation to study this complex quantum system.

Contribution

It introduces a quantum simulation method for a long-range disordered spin glass exhibiting strong area-law violation, enabling exploration of large systems beyond numerical limitations.

Findings

Ground state shows more than logarithmic entanglement entropy violation.

Strong area-law violation occurs even in short-range regime.

Feasible trapped ion implementation with realistic parameters.

Abstract

Ground states of local Hamiltonians are known to obey the entanglement entropy area law. While area law violation of a mild kind (logarithmic) is commonly encountered, strong area-law violation (more than logarithmic) is rare. In this paper, we study the long range quantum spin glass in one dimension whose couplings are disordered and fall off with distance as a power-law. We show that this system exhibits more than logarithmic area law violation in its ground state. Strikingly this feature is found to be true even in the short range regime in sharp contrast to the spinless long range disordered fermionic model. This necessitates the study of large systems for the quantum $XY$ spin glass model which is challenging since these numerical methods depend on the validity of the area law. This situation lends itself naturally for the exploration of a quantum simulation approach. We present a proof-of-principle implementation of this non-trivially interacting spin model using trapped ions and provide a detailed study of experimentally realistic parameters.

Equations

$$S^A=-Tr[{\rho }_A\ln ({\rho }_A)],$$

$$H_{\mathrm{SG}}$$

$$\mathrm{where}{J}_{ij}$$

$$S^A=-\sum _{k=1}^{2^L}{\zeta }_k\ln {\zeta }_k.$$

$$H_F=\sum _{i<j}^N{J}_{ij}(c_i^{†}{c}_j+h.c.).$$

$$S^A=-\sum _k[{\lambda }_k\ln {\lambda }_k+(1-{\lambda }_k)\ln (1-{\lambda }_k)],$$

$$H_{\mathrm{spin}}$$

$$\simeq \sum _{ij}{J}_{ij}\left({\sigma }_i^{+}{\sigma }_j^{-}+h.c.\right).$$

$$J_{ij}={\Omega }_i{\Omega }_j\sum _m{\eta }_{i,m}{\eta }_{j,m}\left[\frac{{\nu }_m}{{\Delta }_s^2-{\nu }_m^2}\right]\simeq \frac{J}{|z_i-z_j{|}^{\alpha }},$$

$$|{\Delta }_s|\gg |\delta |\gg |J_{ij}|\text{for all}i,j.$$

$$H_{SG}=\sum _{i<j=1}^N{J}_{ij}({\sigma }_i^{+}{\sigma }_j^{-}+h.c.),$$

$$\rho =\sum _{l,l^{\prime }=1}^{\text{textsuperscript}N{C}_{N/2}}{\psi }_l{\psi }_{l^{\prime }}^{\ast }|l⟩⟨l^{\prime }|.$$

$${\rho }_A=\sum _{l_A,l_A^{\prime }}{\rho }_A^{l_A,l_A^{\prime }}|l_A⟩⟨l_A^{\prime }|,$$

$$S^A=-\sum _{k=1}^{2^L}{\zeta }_k\ln {\zeta }_k.$$

$$H_F=\sum _{i<j=1}^N{J}_{ij}(c_i^{†}{c}_j+h.c.)$$

$$H_F=\sum _{k=1}^N{ϵ}_k{b}_k^{†}{b}_k,$$

$$|{\Psi }_0⟩=\prod _{k=1}^{N_p}{b}_k^{†}|0⟩$$

$$C_{ij}=Tr[{\rho }_A{c}_i^{†}{c}_j]=Tr[\rho c_i^{†}{c}_j],$$

$${\rho }_A=\frac{e^{-H_A}}{Z},$$

$$H_A=\sum _{k=1}^L{h}_k{a}_k^{†}{a}_k,$$

$${\rho }_A=\frac{exp[-\sum _{k=1}^L{h}_k{a}_k^{†}{a}_k]}{\prod _{k=1}^L[1+exp(-h_k)].}$$

$$C_{ij}=\sum _{k=1}^L{\varphi }_i^{\ast }(k){\varphi }_j(k)\frac{1}{e^{h_k}+1}.$$

$${\lambda }_k=\frac{1}{e^{h_k}+1},$$

$$S^A=-\sum _{k=1}^L[{\lambda }_k\ln {\lambda }_k+(1-{\lambda }_k)\ln (1-{\lambda }_k)].$$

$$\hat{H}=\sum _{i,m}{\Omega }_i{\eta }_{im}{\sigma }^x\cos ({\Delta }_{s}t)\left[a_m{e}^{-i{\nu }_{m}t}+a_m^{†}{e}^{i{\nu }_{m}t}\right].$$

$$\hat{H}=\sum _{i,m}{\Omega }_i{\eta }_{im}{\sigma }_i^r{G}_m(t)$$

$${\sigma }_i^r$$

$$G_m(t)$$

$$U(t)=e^{-iH(t)}=e^{H^{(1)}(t)+H^{(2)}(t)+\dots },$$

$$H^{(2)}(t)={\int }_0^{t}id{t}_1\left(H_I^{(2)}(t_1)+H_{II}^{(2)}(t_1)\right),$$

$$H_I^{(2)}(t_1)$$

$$\times {\int }_0^{t_1}[{\sigma }_i^r(t_1),{\sigma }_j^r(t_2)]G_m(t_1)G_{m^{\prime }}(t_2)d{t}_2$$

Full text

I Introduction

Entanglement in quantum systems has become one of the most studied subjects in physics Horodecki et al. (2009); Amico et al. (2008); Calabrese et al. (2009); Gong et al. (2017). The scaling of entanglement entropy with the spatial size of the subsystem has proved to be a useful diagnostic to characterize different phases across quantum phase transitions Eisert et al. (2010); Laflorencie (2016). Many-particle ground state wavefunctions of local Hamiltonians adhere to the “area-law” of entanglementHastings (2007); Eisert et al. (2010); Laflorencie (2016); Srednicki (1993), according to which, in a subsystem of volume $L^{d}$ in a $d$-dimensional system, entanglement entropy scales as $S_{A}\propto L^{d-1}$, since only the degrees of freedom near the boundary of the two subsystems contribute to the entanglement entropy (Fig. 1). The well-known logarithmic violation of the area law ($S_{A}\propto L^{d-1}\ln L$) has been reported in a variety of translationally invariant models such as gapless systems (e.g. quantum critical point) Vidal et al. (2003); Latorre et al. (2004); Korepin (2004); Jin and Korepin (2004), free fermionic (hopping) models Wolf (2006); Gioev and Klich (2006); Swingle (2010); Gori et al. (2015), bosonic systems Lai et al. (2013); Ding and Yang (2009) and interacting spin chains Peschel (2004); Eisert and Osborne (2006); Its et al. (2005); Popkov and Salerno (2005); Ding et al. (2008); Latorre and Riera (2009); Koffel et al. (2012); Buyskikh et al. (2016); Frérot et al. (2017); Vodola et al. (2015). Quantum dynamics of correlations in the long range spin chains with $\ln L$-kind violations are known to show many exotic features Liu et al. (2019); Safavi-Naini et al. (2019); Smith et al. (2016); Jurcevic et al. (2014),

Strong area-law violations indicate very strong correlations and highly-entangled ground sates could be useful for quantum computation technologies. However, ground states showing stronger-than-logarithmic violations of the area law ($S_{A}\propto L^{\beta}$), particularly in systems that preserve translational symmetry are rare Gori et al. (2015); Deng et al. (2017). There are also exotic special constructions like the Motzkin spin chain Bravyi et al. (2012); Movassagh (2017) where the Hamiltonian is theoretically engineered so as to make a pre-designed area-law-violating state as its ground state. In some models Pouranvari and Yang (2014); Vitagliano et al. (2010); Roy and Sharma (2018), the presence of disorder in the couplings has been shown to give rise to strong area-law violations (more than $\ln L$-kind). Thorough characterisation of a ground state that strongly violates entanglement entropy requires the study of large systems which are computationally intractable. Furthermore, the study of such systems with numerical schemes such as density matrix renormalization group (DMRG) Schollwöck (2005, 2011), time evolving block decimation (TEBD)  Vidal (2004) and other matrix product states based techniques is not feasible as they crucially depend on the area law. This leads to the motivation for simulating such systems in ultra-cold setups.

Inspired by Roy and Sharma (2018), we want to propose a physically realizable long range spin model with disorder in the spin-spin couplings that shows strong area-law violation of the entanglement. Quantum spin models are regularly realized in ultra-cold systems using Rydberg atoms Löw et al. (2009); Schauß et al. (2015); Weimer et al. (2008); Gil et al. (2014); Labuhn et al. (2016); Mukherjee et al. (2016), polar molecules Gorshkov et al. (2011); Manmana et al. (2013); Yan et al. (2013) and trapped ions Sørensen and Mølmer (1999); Deng et al. (2005); Lee et al. (2005); Foss-Feig et al. (2013); Arrazola et al. (2016); Kim et al. (2009); Britton et al. (2012); Islam et al. (2013); Jurcevic et al. (2014). The advantage of using ions is that they are cooled and trapped very efficiently. This allows us to manipulate the collective vibrational modes along with the internal states with great precision which is key to fine tuning the spin-spin couplings. Although most of the earlier works with trapped ions were used to simulate spin models with long range interactions without disorder, there have been some involving disorder in the potential Smith et al. (2016); Trautmann and Hauke (2018). In this work, we consider a system of trapped ions to model our long range spin model with disorder in the spin couplings.

In this article, we study the one dimensional long-range quantum $XY$ spin glass which contains disorder in its couplings. Using exact diagonalization, we calculate the von Neumann entanglement entropy (VNEE) for its ground state and demonstrate strong area-law violation for 14 spins. Typically, one exactly solves the 1D spin chain with nearest-neighbor spin interactions by mapping it into the free fermionic chain via the Jordan-Wigner transformation. This free fermionic mapping is not useful anymore for spin chains with (disordered) long range interactions like in our quantum $XY$ spin glass Hamiltonian because this only manages to convert one difficult problem into another difficult problem. We compare the VNEE of our model with the spinless fermionic model with disorder which is known to reliably exhibit strong area-law violation for large systems Pouranvari and Yang (2014); Roy and Sharma (2018). We find that the VNEE of quantum $XY$ spin glass is more sensitive to the sign of the coupling than the spinless fermionic model. Also interestingly, we find strong area law violation even for short-range interactions for the quantum $XY$ spin glass system in sharp contrast to the spinless fermionic model with disorder. In the final section, we investigate the possibility of realizing the quantum $XY$ spin glass in real ion trap experiments.

II Entanglement entropy for quantum $XY$ spin glass model

The VNEE between subsystem A and the rest of the system B is defined as

[TABLE]

where $\rho_{A}$ is the reduced density matrix of the subsystem A. VNEE can also be described as the first order Renyi entanglement entropy Eisert et al. (2010); Amico et al. (2008); Horodecki et al. (2009). In this context, second-order Renyi entropy has been measured recently in experiments based on ultracold atoms Islam et al. (2013); Kaufman et al. (2016). The long-range quantum $XY$ spin glass in 1D De Cesare et al. (1992); Cesare et al. (1990); Albrecht et al. (1982) for $N$ spins is given by

[TABLE]

Here $\sigma_{i}^{+}$ and $\sigma_{i}^{-}$ are the spin raising and lowering operators respectively, acting on the $i\textsuperscript{th}$ spin. The coupling between $i\textsuperscript{th}$ and $j\textsuperscript{th}$ spins is given by $J_{ij}$ where $J$ is chosen from a uniform distribution of random numbers $[-V,V]$ and $V$ is an arbitrary positive constant. This implies that between any two spins, there can be either an anti-ferromagnetic or ferromagnetic coupling, chosen randomly as the sign of $J_{ij}$ changes randomly. $r_{ij}=|i-j|$ is the distance between sites $i$ and $j$ in a chain and $\alpha$ is the long-range parameter. We assume open boundary conditions throughout.

In order to calculate the VNEE for Eq.(2), we first diagonalize the Hamiltonian exactly in the basis of states that correspond to an equal number of spins in the up state and in the down state. We identify the many-body ground state $\ket{\Psi_{g}}$ as the eigenstate corresponding to the lowest eigenvalue of the Hamiltonian and express the full density matrix as $\rho=\ket{\Psi_{g}}\bra{\Psi_{g}}$. We partition our system of $N$ spins into subsystem A of $L$ spins and subsystem B of $N-L$ spins. The eigenvalues of the reduced density matrix $\rho_{A}=Tr_{B}[\rho]$ are obtained and denoted as $\zeta_{k}$ (appendices for details). Using Eq. 1, the VNEE is then given by

[TABLE]

Our exact numerical calculations for a system size of $N=14$ reveals power-law violation of area-law, which is depicted in Fig. 2(a) as a function of the parameter $\alpha$. The power-law exponents can be extracted by curve fitting, which is shown in Fig. 2(b) for increasing system sizes $N=10,12$ and $14$ respectively, which indicates that the power-law violation is becoming more robust with increasing system sizes. Also we have checked that the power law violation of VNEE is robust against the increasing number of disorder realizations.

In order to study the role of disorder, we considered a few variant models. We found that a disordered ferromagnet (negative random couplings), does not show stronger than $\ln L$-kind violation no matter how long-range the couplings are, as shown in Fig. 3(a). On the other hand the disordered anti-ferromagnet (positive random couplings) does yield power-law violation of the area law if the long-range parameter $\alpha$ is sufficiently small, just like the spin glass as is clear from Fig. 3(b).

III Entanglement entropy for spinless fermonic model

On the other hand, we are able to calculate the VNEE for a very large system for the spinless fermionic model which is given as

[TABLE]

Here $c_{i}^{\dagger}$ ($c_{i}$) is the fermionic creation (annihilation) operator at site $i$. $J_{ij}$ is defined similarly as in Eq.(3). One can calculate VNEE for large systems for a noninteracting spinless fermionic model because the entanglement entropy is obtained from the eigenvalues of the correlation matrix $C_{ij}=\langle c_{i}^{\dagger}c_{j}\rangle$, where $i,j\in A$. This is numerically much easier to calculate as compared to the calculation of the exact reduced density matrix (for further details see Peschel (2003, 2012); Roy and Sharma (2018) and appendics). The reason is that the ground state has a Slater determinant structure to which the Wick’s theorem may be applied. Thus the two-point correlation is sufficient to obtain all the higher order correlations.

The VNEE for the noninteracting fermionic system is

[TABLE]

where $\lambda_{k}$ are the eigenvalues of $C_{ij}$. This fermionic system shows a strong (power law type) area law violation of entanglement entropy as can be seen in Fig. 4(a) Roy and Sharma (2018). The dependences of the power law exponents $\beta$ on the long-range parameter $\alpha$ are shown in Fig. 4(b). The exponent $\beta$ remains non-zero in the regime $0<\alpha\leq 1$ and then it slowly vanishes as the long-range parameter $\alpha$ increases.

Interestingly, when compared to the exponent obtained for the quantum spin glass system, we find that it does not vanish even for large $\alpha$. This means that there is an indication of strong area law violation even in short range interactions but to confirm this one would need to go to larger systems of the quantum spin glass system. Yet another point of distinction between the ground state of the spinless fermion model and the quantum spin glass is its dependence on sign of the coupling. All the results of the fermionic model remain unchanged even if the random coupling $J$ is chosen randomly from the negative distribution $[-1,0]$ or the positive distribution $[0,1]$ whereas the spin model gives different results for small values of $\alpha$ in cases of the negative and positive random couplings as previously discussed. In order to get power law behavior across all the values $\alpha$ shown in Fig. 2(b), we need both ferromagnetic as well as anti-ferromagnetic couplings. These observations suggest that the underlying frustration in these models has a crucial role in the strong area law violation.

IV Quantum simulation of long range disordered spin models

Next we implement the long range disordered Hamiltonian in Eq.(2) using trapped ions. There exist plenty of proposals involving trapped ions realizing various spin models Sørensen and Mølmer (1999); Deng et al. (2005); Lee et al. (2005); Foss-Feig et al. (2013); Arrazola et al. (2016); Kim et al. (2009); Britton et al. (2012); Islam et al. (2013) including disordered potentials Smith et al. (2016); Trautmann and Hauke (2018). However, in this work we focus on realising a spin model with disorder in its spin couplings. As shown in Fig. 5, our setup consists of a one-dimensional chain of 171Yb+ ions aligned along the $z$-direction in a standard linear Paul trap James (1998). The effective spin-1/2 system constitutes of two hyperfine states, ${}^{2}S_{1/2}\ket{F=0,m_{F}=0}$ and ${}^{2}S_{1/2}\ket{F=1,m_{F}=0}$, with an energy splitting of $\omega_{\rm hf}=12.64$ GHz (we assume $\hbar=1$). Using a pair of counter-propagating Raman lasers with wave vectors $\mathbf{k}_{i=1,2}$ and laser frequencies $\omega_{i=1,2}$ along the transverse $x$-axis, the two internal states are coupled via an intermediate excited state $|e\rangle$ which in our case is the ${}^{2}P_{1/2}\ket{F=0,m_{F}=0}$ state. The pair of lasers are far detuned (by $\Delta^{\prime}$) from this excited state (see Fig. 5) where $\Delta^{\prime}$ is in the order of few hundred GHz. We define the relative wave vector as $\mathbf{k}_{L}=\mathbf{k}_{1}-\mathbf{k}_{2}$ and the laser beat frequency as $\omega_{L}=\omega_{2}-\omega_{1}$. The Hamiltonian for the effective optical coupling of the internal states (including motional states) in the interaction picture after applying the rotating wave approximation (RWA), is $\hat{H}=\frac{1}{2}\sum^{N}_{i=1}\Omega_{i}(\sigma^{+}_{i}e^{i(\mathbf{k}_{L}\cdot\hat{x}_{i}-\Delta t+\phi)}+{\rm h.c.})$, where $\Omega_{i}$ is the Rabi frequency for the $i$th ion; $\Delta=\omega_{L}-\omega_{\rm hf}$ is the detuning and $\phi$ is the relative phase.

The ion motional dynamics is a result of the confining linear trap and the ion-ion repulsion. Ions oscillate about their equilibrium positions giving rise to normal modes of oscillation which are classified as axial and transverse modes. We couple to the transverse modes of motion along the x-axis only, matching the direction of propagation of our Raman lasers. The Hamiltonian for motional dynamics is $\hat{H}^{\rm m}=\sum_{m}\nu_{m}a^{\dagger}_{m}a_{m}$ where $\nu_{m}$’s are the transverse mode frequencies and $\hat{a}^{\dagger}_{m}$ and $a_{m}$ are the motional mode operators. We work in the Lamb-Dicke regime, $|\eta_{im}|\sqrt{\langle\hat{a}^{\dagger}_{m}\hat{a}_{m}\rangle}\ll 1$ where $\eta_{im}=k_{L}b_{i,m}\sqrt{\hbar/2M\nu_{m}}$ is the Lamb-Dicke parameter and $b_{i,m}$ is the transverse normal mode transformation matrix for the $i$th ion in the $m$th mode Kim et al. (2010).

Similar to previous proposals for generating spin models Sørensen and Mølmer (1999); Deng et al. (2005); Lee et al. (2005); Foss-Feig et al. (2013); Arrazola et al. (2016); Kim et al. (2009); Britton et al. (2012); Islam et al. (2013), our beat frequency $\omega_{L}$ is simultaneously detuned to red and blue sidebands of the highest normal mode frequency which in our case is $\nu_{x}$. We define $\Delta_{b}=\Delta+\nu_{x}$ and $\Delta_{r}=\Delta-\nu_{x}$. For symmetrical detuning, $\Delta_{s}=\Delta_{b}=-\Delta_{r}$, we obtain the well known Ising model Sørensen and Mølmer (1999); Islam et al. (2013). However, by adding a small asymmetry $\delta$ to the effective detuning $\Delta_{s}$ (see appendices), we get the following effective spin Hamiltonian

[TABLE]

Here we have adiabatically eliminated the motional modes by assuming $|\Delta^{\pm}_{m}|\gg|\eta_{im}\Omega_{i}|$ where $\Delta^{\pm}_{m}=(\Delta_{s}\pm\nu_{m})$. The spin coupling is given as

[TABLE]

where $z_{i}$ is the equilibrium position of the $i$th ion. The explicit expression for $B_{i}$ is given in the appendix. Eq.(8) is a good approximation to Eq.(7) provided

[TABLE]

For our numerics, we consider a linear Paul trap with an axial trap frequency $\nu_{z}/(2\pi)=1.0$ MHz and transverse trapping frequencies, $\nu_{y}/(2\pi)\simeq\nu_{x}/(2\pi)=3.95$ MHz. The Rabi frequency is kept constant at $\Omega_{i}/2\pi=\Omega_{j}/2\pi=50$ kHz and the wave numbers are $k_{1}=k_{2}=2\pi/369.75$ nm*-1*. Thus for 171Yb+ ions we get $\eta_{im}\simeq 0.01-0.03$ giving $|\eta_{im}\Omega_{i}/2\pi|\simeq 500$ Hz. We have $\Delta_{s}/2\pi=3.5-4.5$ MHz and $\delta/2\pi\simeq 100$ kHz, for which the above conditions are easily satisfied.

Usually the spin models are generated by choosing the detuning far away from the entire bunch of mode frequencies to obtain spin couplings of the same sign. However by tuning close to the center of mass frequency along the $x$ axis, we introduce disorder in a controlled manner. To implement disordered spin interactions, we choose the bichromatic beat frequency to be slightly below and above $\nu_{x}$ which gives the random spin couplings that we need for our model. A spatially dependent random intensity profile constructed using spatial light modulation methods Nogrette et al. (2014) or speckled intensity techniques McGehee et al. (2013) can also be used to create disordered spin couplings. To obtain different power laws in $J_{ij}$ as a function of distance, we scan the detuning $\Delta_{s}$ over a certain range. Fig. 6 characterizes the power law behaviour and the disorder in our spin couplings given by Eq.(9). The power law scaling in Fig. 6(a) is obtained after averaging over the couplings over different ions. The power law scaling over distance in Eq.(9) is usually an empirical approximation and our fittings are limited by finite size effects. Particularly for small system sizes, the power law scaling is highly sensitive to ion equilibrium positions. A pair of ions near the center of the chain are closer to each other as opposed to those at the edge. The typical inter-ionic distances for our chosen trap frequencies ranges from $1.4$ $\mu$m to $2.0$ $\mu$m. Large number of ions in a linear Paul trap can lead to unstable transverse motion and unwanted zig-zag patterns Schiffer (1993). Including a larger number of ions requires having lower values of $\epsilon=\nu_{z}/\nu_{x}$ which are determined by experimental constraints. The scenario involving a two dimensional Penning trap is more relaxed in this regard and good power law fits are obtained for hundreds of ions Britton et al. (2012) but have yet to be shown for disordered couplings. Fig. 6(b) and (c) depict the disorder in the couplings for different values of $\alpha$. One clearly sees the random fluctuation in sign and magnitude of $J_{ij}$ for different values of $i,j$. Although, one can couple to axial modes as well which are typically at a much lower frequency and are more widely separated compared to the transverse modes, our calculations seem to suggest that there is no particular advantage in choosing either type of mode for small systems.

V Conclusion

To summarize, we have studied the entanglement properties of the ground state of the long-range one dimensional quantum $XY$ spin glass. We find that it exhibits a strong area-law violation, which in contrast with the disordered long-range fermionic hopping model, extends into the short-range regime. Such strong area-law violation spells doom for numerical techniques to be effective beyond extremely small system sizes. This motivates the implementation of a quantum simulation approach to access such models. We thus consider trapped ions to realize our spin Hamiltonian which can be generalized to larger systems in future experiments. Our work also indicates the possibility of achieving strong area-law violating ground states using disordered long ranged couplings. Quantum spin models with disordered long range couplings have not been explored much so far either theoretically or experimentally. We hope that our work will help boost more interseting research engaging such spin models.

We are grateful to the High Performance Computing (HPC) facility at IISER Bhopal, where large-scale calculations in this project were run. N.R is grateful to University Grants Commission (UGC), India for support via his Ph.D fellowship. R.M would like to acknowledge J. Lishman for his useful comments and the Max-Planck society for funding under the MPG-IISER partner group program and . A.S acknowledges SERB for the startup grant (File Number: YSS/2015/001696).

Appendix A: Methodology to calculate entanglement entropy for quantum spin glass

The spin glass Hamiltonian is given by

[TABLE]

where $N$ is the total number of spins with spin-$\frac{1}{2}$. The Hamiltonian matrix is written in the $\sigma_{z}$ basis ($\ket{\uparrow}$ and $\ket{\downarrow}$) where total magnetization is zero i.e. number of up spins is equal to the number of down spins. $H_{SG}$ is exactly diagonalized in the basis and we obtain the normalized many-body ground state $\ket{\Psi_{g}}=\sum_{l}\psi_{l}\ket{l}$. The corresponding density matrix is given by $\rho=\ket{\Psi_{g}}\bra{\Psi_{g}}$, which can be written as

[TABLE]

We consider the first $L$ (where $L\leq N/2$) spins as the subsystem A and the remaining $N-L$ spins as the subsystem B. Each element in the basis then can be written as $\ket{l}=\ket{l_{A}}\ket{l_{B}}$, where $\ket{l_{A}}$ and $\ket{l_{B}}$ are the orthonormal sets in subsystem A and B respectively. The reduced density matrix defined as $\rho_{A}=Tr_{B}[\rho]$, can be written as

[TABLE]

where $\rho_{A}^{l_{A},l_{A}^{\prime}}$ is the $(l_{A},l_{A}^{\prime})\textsuperscript{th}$ element of matrix $\rho_{A}$. One then diagonalizes the matrix $\rho_{A}$ exactly and the von Neumann entanglement entropy $S^{A}=-Tr[\rho_{A}\ln\rho_{A}]$ is computed with

[TABLE]

where $\zeta_{k}$’s are the eigenvalues of $\rho_{A}$.

Appendix B: Methodology to calculate entanglement entropy for spinless fermionic model

The fermionic Hamiltonian is given by

[TABLE]

whose diagonal form is

[TABLE]

where $b_{k}=\sum\limits_{j=1}^{N}\psi_{j}(k)c_{j}$. The many-body fermionic ground state is

[TABLE]

Due to the Slater determinant structure of $\ket{\Psi_{0}}$, all higher correlations can be obtained by two point correlation $C_{ij}=\langle c_{i}^{\dagger}c_{j}\rangle$Peschel (2003); Peschel and Eisler (2009); Peschel (2012). The density matrix of the full system is $\rho=\ket{\Psi_{0}}\bra{\Psi_{0}}$ and the reduced density matrix of subsystem A is $\rho_{A}=Tr_{B}(\rho)$. By definition a one particle function, in this case two-point correlation in the subsystem, can be written as

[TABLE]

where $i,j\in A$. The subsystem density matrix of an eigenstate of a quadratic fermionic Hamiltonian is thermalPeschel (2003)

[TABLE]

where $H_{A}=\sum\limits_{i,j=1}^{L}H_{ij}^{A}c_{i}^{\dagger}c_{j}$ is called the entanglement Hamiltonian, and $Z$ is obtained to satisfy the condition $Tr[\rho_{A}]=1$. The entanglement Hamiltonian can thus be written in the diagonal form as,

[TABLE]

where $a_{k}=\sum_{j=1}^{L}\phi_{j}(k)c_{j}$. The reduced density matrix is then given by

[TABLE]

Using Eq.21, we can write Eq.18 as,

[TABLE]

This shows that the matrices $C$ and $H_{A}$ share the eigenstate $\ket{\phi_{k}}$ and their eigenvalues are related by

[TABLE]

where $\lambda_{k}$’s are eigenvalues of matrix $C$ in the subsystem.

The von Neumann entanglement entropy can be simplified Sharma and Rabani (2015) using Eq.21 and Eq.23 as

[TABLE]

Appendix C: Effective spin model using time dependent perturbation theory

Specific spin-spin interactions are induced by the excitation lasers and trapping conditions. Our bichromatic beat frequency is tuned symmetrically to the red and blue sidebands simultaneously and thus we have $\Delta_{s}=\Delta_{b}=-\Delta_{r}$. In the Lamb-Dicke regime, we expand $\mathbf{k}_{L}\cdot\hat{x}_{i}$ in terms of the motional operators to obtain the well known spin dependent Hamiltonian Foss-Feig et al. (2013); Arrazola et al. (2016); Kim et al. (2009); Britton et al. (2012); Islam et al. (2013)

[TABLE]

Here we have assumed that the relative phase for the pair lasers is $-\pi/2$. As discussed in the text, we can make our detuning asymmetric by shifting the levels by $\delta$, i.e. $\Delta_{s}+\delta=-\Delta_{r}+\delta=\Delta_{b}+\delta$, which essentially implies that our ladder operators rotate in the manner: $\sigma^{\pm}\rightarrow\sigma^{\pm}e^{\mp i\delta t}$ and thus our Hamiltonian gets modified to

[TABLE]

where we define

[TABLE]

and $\Delta^{\pm}_{m}=\Delta_{s}\pm\nu_{m}$. For the parameters chosen in our article, typical values of $\Delta^{\pm}_{m}$ are slightly less than one MHz. Using the Magnus expansion, we can write the time evolution operator as follows

[TABLE]

where $H^{(1)}(t)=-i\int^{t}_{0}H(t_{1})dt_{1}$ and $H^{(2)}(t)=-\frac{1}{2}\int^{t}_{0}\int^{t_{1}}_{0}[H(t_{1}),H(t_{2})]dt_{1}dt_{2}$. $H^{(1)}(t)$ corresponds to the spin-phonon coupling and can be ignored if we have $\delta\pm\Delta^{\pm}_{m}\gg|\eta_{im}\Omega_{i}|$. For our parameters, $\eta_{im}\simeq 0.01-0.03$ giving $|\eta_{im}\Omega_{i}|\simeq 500$ Hz. The smallest value of $\Delta^{\pm}_{m}$ is around $880$ kHz and $\delta$ can be chosen between 10 to 100 kHz for the spin-phonon coupling to be ignored. The second order term gives our effective spin model which can be split into two parts as follows

[TABLE]

where

[TABLE]

Using the commutator relationships

[TABLE]

and solving the integrals over time gives

[TABLE]

where the definition for the spin-spin coupling is

[TABLE]

and the longitudinal field is

[TABLE]

The approximated expressions of $J_{ij}$ and $B_{i}$ are obtained by assuming $|\Delta_{m}^{\pm}|>>\delta$.