Simulating generic spin-boson models with matrix product states

The global coupling of few-level quantum systems (“spins”) to a discrete set of bosonic modes is a key ingredient for many applications in quantum science, including large-scale entanglement generation, quantum simulation of the dynamics of long-range interacting spin models, and hybrid platforms for force and spin sensing. We present a general numerical framework for treating the out-of-equilibrium dynamics of such models based on matrix product states. Our approach applies for generic spin-boson systems: it treats any spatial and operator dependence of the two-body spin-boson coupling and places no restrictions on relative energy scales. We show that the full counting statistics of collective spin measurements and infidelity of quantum simulation due to spin-boson entanglement, both of which are difficult to obtain by other techniques, are readily calculable in our approach. We benchmark our method using a recently developed exact solution for a particular spin-boson coupling relevant to trapped ion quantum simulators. Finally, we show how decoherence can be incorporated within our framework using the method of quantum trajectories, and study the dynamics of an open-system spin-boson model with spatially nonuniform spin-boson coupling relevant for trapped atomic ion crystals in the presence of molecular ion impurities.

Figure 1

Spins globally coupled to bosonic modes. A schematic description of ${\mathcal{N}}_{\text{s}}=7$ spins coupled globally to ${\mathcal{N}}_{\text{b}}=3$ bosonic modes (colored lines). Position in the array corresponds to, e.g., position of the spins in real space. The amplitude of the boson line denotes the value of the spin-boson coupling at that position. Our approach applies for any spatial dependence of the bosonic mode coupling.

Figure 2

(a) Tensor network diagram representation of an MPS with open boundary conditions. Lines extending to the top of the page correspond to physical states $|j_i\rangle$ of local system $i$. (b) Tensor network diagram representation of an MPO with open boundary conditions. Each tensor has two physical indices, corresponding to a linear transformation (operator).

Figure 3

Application of spin-boson propagator through swap gates. The tensor network diagram for the application of $e^{-i\delta t{\widehat{H}}_{\mathrm{s}-\mathrm{b}}}$ on a state of ${\mathcal{N}}_{\text{b}}=2$ bosons (circles) and ${\mathcal{N}}_{\text{s}}=3$ spins (squares) through the Trotter-Suzuki decomposition with swap gates (crossed lines) demonstrates that only nearest-neighbor operations are required. For clarity, circles containing the boson mode number are drawn next to boson index lines, and the boson lines have been color coded.

Figure 4

MPS topology for the long-range time evolution method with ${\mathcal{N}}_{\text{s}}={\mathcal{N}}_{\text{b}}=\text{3}$. To avoid variational metastability when directly simulating time evolution with long-ranged spin-boson interactions, we use an MPS representation in which spins are always neighbored by bosons, allowing for strong spin-boson entanglement to be built up locally. The coherence between spins and bosons which are not neighboring is built up both by repeated variational sweeping and by mixing in a small amount of a state which contains long-range spin-boson coherence (see the text for details).

Figure 5

Tensor network representation of the characteristic function. The tensor network description of the characteristic function for a given Fourier variable $k$ and spin quadrature $\mathbf{n}$ is especially simple due to the fact it is the expectation of a product operator. Here, ${\widehat{O}}_j=exp\left[\pi ik{\widehat{\mathbf{\sigma }}}_j·\mathbf{n}/({\mathcal{N}}_{\mathrm{s}}+1)\right]$.

Figure 6

Tensor network representation of the fidelity between a traced spin-boson state and a pure spin state. An example tensor network diagram giving the (squared) fidelity of the spin density matrix ${\widehat{\rho }}_{\text{b}}$ of a system with three bosons (circles) and three spins (squares) with the bosons traced out with the pure density matrix $|\mathrm{\Psi }\rangle \langle \mathrm{\Psi }|$ of a system comprised of only spins (triangles). Assuming the boson tensors are in left canonical form, the upper diagram is transformed to the lower diagram, which is the inner product of a vector $\mathcal{L}$.

Figure 7

Equilibrium structures and effective spin-spin couplings for 1D and 2D ion traps. (a) The equilibrium crystal structure for a 1D Paul trap with ${\mathcal{N}}_{\text{s}}=11$ ions and a trap anisotropy of ${\omega }_z/{\omega }_x=0.184$. The spacing between ions is nonuniform, with the smallest spacing at the center of the chain. (b) Equilibrium positions for ${\mathcal{N}}_{\text{s}}=61$ ions in a Penning trap, using the trap parameters of Ref. [23]. Here, ${\ell }_0={(2e^2/4\pi {\varepsilon }_{0}M{\omega }_z^2)}^{1/3}$. The ions roughly form a triangular lattice. (c) Spin-spin couplings, Eq. (21), for the 1D trap with detuning $\delta /h=80\text{kHz}$ above the center-of-mass mode (blue circles) and 150 kHz above the center of mass (red squares), with the force adjusted so that the maximum value of $J_{i,j}$ is roughly 400 Hz. Best power law fits are shown with solid lines. Significant deviations from translational invariance are seen, especially for the larger detuning. (d) Spin-spin couplings, Eq. (21), for the 2D trap with a detuning of $\delta /h=2\text{kHz}$ above the center of mass and a spin-dependent force similar to Ref. [23]. While there is some nonuniformity, the interactions are well modeled by an all-to-all interaction, justifying the use of only the center-of-mass mode in the main text.

Figure 8

Comparison of numerical and analytically computed collective spin observables; multimode case. (a) Interaction-induced depolarization of the collective spin computed numerically (red symbols) and analytically (blue lines) for the $\delta /h=80\text{kHz}$ detuned configuration of Fig. 7, showing excellent agreement. The inset shows small oscillations due to spin-phonon entanglement, which are not captured within spin-only approaches. (b) Collective spin correlations for the $\delta /h=150\text{kHz}$ detuned configuration of Fig. 7. The difference between the numerical values (symbols) and analytical results (lines) is not visible on the scale of the plot.

Figure 9

Fidelity of spin-boson dynamics with pure spin dynamics. The fidelity of the spin-boson dynamics with the bosons traced out with the pure spin dynamics, calculated using the methods of Sec. 2e, is displayed for the $\delta /h=80\text{kHz}$ detuning (red $+$'s) and the $\delta /h=150\text{kHz}$ detuning (blue $\times$'s) configurations of Fig. 7. Rapid oscillations resulting from the development of spin-phonon entanglement with all phonon modes is visible, but the overall fidelity remains quite high.

Figure 10

Squeezing dynamics; single-mode case. Dynamics of the Ramsey squeezing parameter in dB for the coupled spin-boson system (oscillating red points) and the pure spin system (blue points) evolving under Eq. (21), with parameters from Fig. 7. Coupling to phonons as well as time dependence of the Ising spin-spin couplings [see Eq. (20)] cause modulation of the squeezing in the spin-boson dynamics. The spin-boson and pure spin approaches agree at the decoupling points, here given by integer multiples of 0.5 ms.

Figure 11

Full counting statistics of spin-boson and pure spin states. The full counting statistics, obtained using the methods of Sec. 2e, are displayed for the squeezed quadrature minimizing the spin variance (upper panels) as well as the antisqueezed quadrature along ${\widehat{S}}^y$ (lower panels) for the parameters of Fig. 7. The left panels display the full coupled spin-boson dynamics, and the right panels correspond to spins evolving under a time-independent Ising model. Periodic modulations from spin-phonon entanglement broaden the squeezed quadrature and can either broaden or narrow features in the antisqueezed quadrature. Generally, the antisqueezed quadrature is less sensitive to phonon coupling than the squeezed quadrature.

Figure 12

Dynamics of MPS bond dimension in 1D and 2D. Panels (a) and (b) show the dynamics of the spin-boson MPS bond dimension (solid black) and the spin MPS bond dimension (magneta) in the 2D case for the parameters of Fig. 10. Panel (a) uses the spatially homogenous center-of-mass mode and panel (b) one of the inhomogeneous tilt modes, leading to the Ising spin-spin coupling constants $J_{j{j}^{\prime }}$ between sites $j$ and $j^{\prime }$ in the MPS representation shown in the insets. The bond dimensions of the spin-boson and spin states have the same overall trend at longer times. Panel (c) shows the 1D, multimode (${\mathcal{N}}_{\mathrm{b}}={\mathcal{N}}_{\mathrm{s}}$) case for the parameters of Fig. 8 and multiple ${\mathcal{N}}_{\mathrm{s}}$. The addition of more boson modes causes a substantial increase in bond dimension due to spin-boson entanglement compared to the spin-only state. For comparison, the fractional decay of the magnetization due to spin-spin interactions $\sim 0.95$ at 0.03 ms in panel (c) and 1 ms in panel (a).

Figure 13

Equilibrium structures and center-of-mass mode of an ion crystal with impurities. (a) Equilibrium crystal structure for 37 atomic ${\mathrm{Be}}^{+}$ ions (empty circles) and 44 ${\mathrm{BeH}}^{+}$ molecular ions (filled circles) in a Penning trap. The heavier molecular ions move to the outside of the crystal due to centrifugal forces in the rotating frame of the ion crystal. (b) Amplitude of the highest frequency axial normal mode, analogous to the center-of-mass mode in a crystal without impurities. The molecular ions have a reduced amplitude due to their larger mass, and the amplitudes of the atomic ions (inside boxed region) are spatially nonuniform, which in turn makes the spin-phonon coupling spatially inhomogeneous.

Figure 14

Spin correlations in an ion crystal with impurities and decoherence. (a) Normalized collective magnetization $\langle {\sigma }^x\rangle$ and its square $\langle {\left({\sigma }^x\right)}^2\rangle$ computed for the spatially inhomogeneous mode pictured in Fig. 13 (blue dashed) and its closest uniform approximation (red solid). The black dashed line is the mean-field prediction of the decay of magnetization due to decoherence. (b) Spatially resolved magnetization $\langle {\sigma }^x(x,y)\rangle$ at the final time shown in (a), with circle size indicating magnitude. Blue, dashed border (red, solid border) circles are the inhomogeneous mode (uniform approximation). Note that all red circles have the same diameter, and only the blue circles are inhomogeneous. The spread in magnetization values is roughly twice the mean for the inhomogeneous mode.