Short-time expansion of Heisenberg operators in open collective quantum spin systems

TL;DR

This paper introduces a short-time expansion method for Heisenberg operators in large open collective spin systems, enabling accurate calculations of quantum properties under decoherence relevant for quantum metrology and information scrambling.

Contribution

The authors develop a Taylor expansion approach for Heisenberg operators in collective spins, allowing efficient short-time analysis in large, decohering quantum systems, including extensions to spin-boson models.

Findings

Accurately computed spin squeezing and correlation functions for 10^4 spins.

Demonstrated effectiveness in strong-decoherence regimes inaccessible to previous methods.

Applicable to systems relevant for trapped ion and cavity QED experiments.

Abstract

We present a new method to compute short-time expectation values in large collective spin systems with generic Markovian decoherence. Our method is based on a Taylor expansion of a formal solution to the equations of motion for Heisenberg operators. This expansion can be truncated at finite order to obtain virtually exact results at short times that are relevant for metrological applications such as spin squeezing. In order to evaluate the expansion for Heisenberg operators, we compute the relevant structure constants of a collective spin operator algebra. We demonstrate the utility of our method by computing spin squeezing, two-time correlation functions, and out-of-time-ordered correlators for $10^4$ spins in strong-decoherence regimes that are otherwise inaccessible via existing numerical methods. Our method can be straightforwardly generalized to the case of a collective spin coupled to bosonic modes, relevant for trapped ion and cavity QED experiments, and may be used to investigate short-time signatures of quantum chaos and information scrambling.

Equations

$$\hat{\mathcal{O}}=\sum _{\mathbf{m}}{\mathcal{O}}_{\mathbf{m}}{\hat{\mathcal{S}}}_{\mathbf{m}}$$

$$\frac{d}{dt}⟨\hat{\mathcal{O}}\left(t\right)⟩=⟨\check{T}\hat{\mathcal{O}}\left(t\right)⟩=\sum _{\mathbf{m},\mathbf{n}}⟨{\hat{\mathcal{S}}}_{\mathbf{m}}⟩T_{\mathbf{m}\mathbf{n}}{\mathcal{O}}_{\mathbf{n}}\left(t\right)$$

$$\check{T}{\hat{\mathcal{S}}}_{\mathbf{n}}\equiv i{\left[\hat{H},{\hat{\mathcal{S}}}_{\mathbf{n}}\right]}_{-}+\sum _{\mathcal{J}}{\gamma }_{\mathcal{J}}\check{\mathcal{D}}\left(\mathcal{J}\right){\hat{\mathcal{S}}}_{\mathbf{n}}=\sum _{\mathbf{m}}{\hat{\mathcal{S}}}_{\mathbf{m}}{T}_{\mathbf{m}\mathbf{n}},$$

$$\check{\mathcal{D}}\left(\mathcal{J}\right)\hat{\mathcal{O}}\equiv \sum _{\hat{J}\in \mathcal{J}}\left({\hat{J}}^{†}\hat{\mathcal{O}}\hat{J}-\frac{1}{2}{\left[{\hat{J}}^{†}\hat{J},\hat{\mathcal{O}}\right]}_{+}\right).$$

$$\frac{d}{dt}⟨\hat{\mathcal{O}}\left(t\right)⟩\to \sum _{w\left(\mathbf{m}\right)<W}⟨{\hat{\mathcal{S}}}_{\mathbf{m}}⟩\sum _{\mathbf{n}}{T}_{\mathbf{m}\mathbf{n}}{\mathcal{O}}_{\mathbf{n}}\left(t\right)$$

$$⟨\hat{\mathcal{O}}\left(t\right)⟩=⟨e^{t\check{T}}\hat{\mathcal{O}}\left(0\right)⟩=\sum _{k\ge 0}\frac{t^k}{k!}\sum _{\mathbf{m},\mathbf{n}}⟨{\hat{\mathcal{S}}}_{\mathbf{m}}⟩T_{\mathbf{m}\mathbf{n}}^k{\mathcal{O}}_{\mathbf{n}}\left(0\right),$$

$$T_{\mathbf{m}\mathbf{n}}^0$$

$$T_{\mathbf{m}\mathbf{n}}^1$$

$$T_{\mathbf{m}\mathbf{n}}^{k>1}$$

$$⟨\hat{\mathcal{O}}\left(t\right)⟩\to \sum _{k=0}^M\frac{t^k}{k!}\sum _{\mathbf{m},\mathbf{n}}⟨{\hat{\mathcal{S}}}_{\mathbf{m}}⟩T_{\mathbf{m}\mathbf{n}}^k{\mathcal{O}}_{\mathbf{n}}\left(0\right).$$

$$H_{\text{OAT}}$$

$$H_{\text{TAT}}$$

$$H_{\text{TNT}}$$

$${\xi }^2\equiv \frac{{\left(\Delta {\varphi }_{\text{min}}\right)}^2}{{\left(\Delta {\varphi }_{\text{polarized}}\right)}^2}=\frac{N}{|⟨\hat{\mathbf{S}}⟩{|}^2}\underset{\begin{array}{c}\mathbf{v}\perp ⟨\hat{\mathbf{S}}⟩\\ \mathbf{v}\cdot \mathbf{v}=1\end{array}}{\min }⟨{\left(\hat{\mathbf{S}}\cdot \mathbf{v}\right)}^2⟩,$$

$$C\left(t\right)\equiv \frac{1}{S}\left(⟨{\hat{S}}_{+}\left(t\right){\hat{S}}_{-}\left(0\right)⟩-⟨{\hat{S}}_{+}\left(t\right)⟩⟨{\hat{S}}_{-}\left(0\right)⟩\right),$$

$$D\left(t\right)\equiv \frac{1}{S^2}{⟨{\left[{\hat{S}}_{+}\left(t\right),{\hat{S}}_{-}\left(0\right)\right]}_{-}^{†}{\left[{\hat{S}}_{+}\left(t\right),{\hat{S}}_{-}\left(0\right)\right]}_{-}⟩}_{\text{nn}},$$

$${\left[s_{\text{z}}^{(j)},s_{\mu }^{(k)}\right]}_{-}$$

$${\left[s_{\mu }^{(j)},s_{\bar{\mu }}^{(k)}\right]}_{-}$$

$${\left(\mu S_{\text{z}}\right)}^m{s}_{\mu }^{(j)}={\left(\mu S_{\text{z}}\right)}^{m-1}{s}_{\mu }^{(j)}\left(1+\mu S_{\text{z}}\right)={\left(\mu S_{\text{z}}\right)}^{m-2}{s}_{\mu }^{(j)}{\left(1+\mu S_{\text{z}}\right)}^2=\cdots =s_{\mu }^{(j)}{\left(1+\mu S_{\text{z}}\right)}^m,$$

$$\mu s_{\text{z}}^{(j)}{S}_{\mu }^m=S_{\mu }\mu s_{\text{z}}^{(j)}{S}_{\mu }^{m-1}+s_{\mu }^{(j)}{S}_{\mu }^{m-1}=\cdots =S_{\mu }^m\mu s_{\text{z}}^{(j)}+m{s}_{\mu }^{(j)}{S}_{\mu }^{m-1},$$

$${\left(\mu S_{\text{z}}\right)}^m{S}_{\mu }=S_{\mu }{\left(1+\mu S_{\text{z}}\right)}^m,$$

$${\left(\mu S_{\text{z}}\right)}^{\ell }{S}_{\mu }^m={\left(\mu S_{\text{z}}\right)}^{\ell -1}{S}_{\mu }^m\left(m+\mu S_{\text{z}}\right)={\left(\mu S_{\text{z}}\right)}^{\ell -2}{S}_{\mu }^m{\left(m+\mu S_{\text{z}}\right)}^2=\cdots =S_{\mu }^m{\left(m+\mu S_{\text{z}}\right)}^{\ell }.$$

$${\left(\nu S_{\text{z}}\right)}^{\ell }{S}_{\mu }^m=S_{\mu }^m{\left(\mu \nu m+\nu S_{\text{z}}\right)}^{\ell },$$

$$S_{\mu }^m{s}_{\bar{\mu }}^{(j)}$$

$$=S_{\mu }^{m-2}{s}_{\bar{\mu }}^{(j)}{S}_{\mu }^2+S_{\mu }^{m-2}2\mu s_{\text{z}}^{(j)}{S}_{\mu }+S_{\mu }^{m-1}2\mu s_{\text{z}}^{(j)}$$

$$=s_{\bar{\mu }}^{(j)}{S}_{\mu }^m+\sum _{k=0}^{m-1}{S}_{\mu }^{k}2\mu s_{\text{z}}^{(j)}{S}_{\mu }^{m-k-1}.$$

$$\sum _{k=0}^{m-1}{S}_{\mu }^{k}2\mu s_{\text{z}}^{(j)}{S}_{\mu }^{m-k-1}$$

$$=m\left(m-1\right)s_{\mu }^{(j)}{S}_{\mu }^{m-2}+m{S}_{\mu }^{m-1}2\mu s_{\text{z}}^{(j)},$$

$$S_{\mu }^m{s}_{\bar{\mu }}^{(j)}=s_{\bar{\mu }}^{(j)}{S}_{\mu }^m+m\left(m-1\right)s_{\mu }^{(j)}{S}_{\mu }^{m-2}+m{S}_{\mu }^{m-1}2\mu s_{\text{z}}^{(j)},$$

$$S_{\mu }^m{S}_{\bar{\mu }}=S_{\bar{\mu }}{S}_{\mu }^m+m{S}_{\mu }^{m-1}\left(m-1+2\mu S_{\text{z}}\right).$$

$$S_{\mu }^m{S}_{\bar{\mu }}^n=S_{\mu }^{m-1}{S}_{\bar{\mu }}^n{S}_{\mu }+n\left[S_{\mu }^{m-1}{S}_{\bar{\mu }}^{n-1}\left(1-n+2\mu S_{\text{z}}\right)\right]=S_{\bar{\mu }}^n{S}_{\mu }^m+n\sum _{k=0}^{m-1}{S}_{\mu }^{m-k-1}{S}_{\bar{\mu }}^{n-1}\left(1-n+2\mu S_{\text{z}}\right)S_{\mu }^k,$$

$${\left[S_{\mu }^m,S_{\bar{\mu }}^n\right]}_{-}=C_{mn;\mu }\equiv n\sum _{k=0}^{m-1}{S}_{\mu }^{m-k-1}{S}_{\bar{\mu }}^{n-1}\left(1-n+2\mu S_{\text{z}}\right)S_{\mu }^k.$$

$$C_{mn;\mu }$$

$$=D_{mn;\mu }+n\sum _{k=0}^{m-2}{C}_{m-k-1,n-1;\mu }{S}_{\mu }^k\left(2k+1-n+2\mu S_{\text{z}}\right),$$

$$D_{mn;\mu }\equiv mn{S}_{\bar{\mu }}^{n-1}{S}_{\mu }^{m-1}\left(m-n+2\mu S_{\text{z}}\right).$$

Full text

I Introduction

Collective spin systems are a versatile resource in quantum science for a range of applications including quantum-enhanced metrology and quantum simulation. The study of such systems dates back to the mid-twentieth century with the introduction of the Dicke modelDicke (1954) that describes atoms cooperatively interacting with a single mode of a radiation field, and the Lipkin-Meshkov-Glick (LMG) model, a toy model for testing many-body approximation methods in contemporary nuclear physicsLipkin et al. (1965); Meshkov et al. (1965); Glick et al. (1965). On the experimental side, the development of advanced trapping, cooling, and control techniques in atomic, molecular, and optical (AMO) systems have enabled the realization of collective spin models in a broad range of platforms, including cold atomic gassesTakano et al. (2009); Appel et al. (2009), Bose-Einstein condensatesKlinder et al. (2015); Estève et al. (2008); Riedel et al. (2010); Gross et al. (2010), ultracold Fermi gassesMartin et al. (2013); Bromley et al. (2018); Smale et al. (2019), trapped ionsBohnet et al. (2016), and optical cavitiesSchleier-Smith et al. (2010a); Chen et al. (2011); Baumann et al. (2010); Leroux et al. (2010); Bohnet et al. (2014); Cox et al. (2016); Hosten et al. (2016a, b); Norcia et al. (2018); Ritsch et al. (2013), among others. These implementations compliment innumerable theoretical studies in a variety of rich subjects, including quantum phase transitions and criticalityLatorre et al. (2005); Alcalde et al. (2007); Wang et al. (2012); Majd et al. (2014), non-equilibrium phenomenaWalls et al. (1978); Morrison and Parkins (2008a, b, c); Kessler et al. (2012); Bhattacherjee (2014); Zhiqiang et al. (2017); Lang et al. (2018), and precision mentrologyWineland et al. (1992); Kitagawa and Ueda (1993); Zhong (2010); Schleier-Smith et al. (2010b); Ma et al. (2011); Huang et al. (2015a); Muessel et al. (2015); Huang et al. (2015b); Hu et al. (2017); Mirkhalaf et al. (2018); Lewis-Swan et al. (2018); He et al. (2019).

One of the primary motivations for studying collective spin systems is their application to quantum-enhanced metrology. Quantum projection noise limits the error $\Delta\phi$ in the measurement of a phase angle $\phi$ with $N$ independent spins to $\Delta\phi\sim 1/\sqrt{N}$Wineland et al. (1992); Itano et al. (1993). Collective spin systems provide a means to break through this limit via the preparation of many-body entangled states such as spin-cat statesAgarwal et al. (1997); Lau et al. (2014); Huang et al. (2015b) and most notably spin-squeezed statesWineland et al. (1992); Kitagawa and Ueda (1993); Ma et al. (2011) that allow for measurement errors $\Delta\phi\sim 1/N^{\varepsilon}$ with $1/2<\varepsilon\leq 1$, where $\varepsilon=1$ saturates the Heisenberg limitZwierz et al. (2010). Such entangled states can be prepared either via heralded methods such as quantum nondemolition measurementsTakano et al. (2009); Appel et al. (2009); Schleier-Smith et al. (2010a); Chen et al. (2011), or via deterministic methods that require nonlinear dynamics, typically realized with phonon-mediatedBohnet et al. (2016), photon-mediatedKlinder et al. (2015); Baumann et al. (2010); Leroux et al. (2010); Bohnet et al. (2014); Cox et al. (2016); Hosten et al. (2016a, b); Norcia et al. (2018); Ritsch et al. (2013) or collisionalEstève et al. (2008); Riedel et al. (2010); Gross et al. (2010); Martin et al. (2013); Bromley et al. (2018); Smale et al. (2019) interactions. Although a truly collective spin model requires uniform, all-to-all interactions, as long as measurements do not distinguish between constituent particles, even non-uniform systems may be effectively described by a uniform model with renormalized parametersHu et al. (2015).

In the absence of decoherence, permutation symmetry and total spin conservation divide the total Hilbert space of a collective spin system into superselection sectors that grow only linearly with system size $N$, thereby admitting efficient classical simulation of its dynamics. Decoherence generally violates total spin conservation and requires the use of density operators, increasing the dimension of accessible state space to $O\left(N^{3}\right)$Hartmann (2016); Xu et al. (2013). In this case, exact simulations can be carried out for $N\lesssim 100$ particles. If decoherence is sufficiently weak, dynamics can be numerically solvable for $N\lesssim 10^{5}$ particles via “quantum trajectory” Monte Carlo methodsPlenio and Knight (1998); Zhang et al. (2018) (also known as “quantum jump” or “Monte Carlo wavefunction” methods) that can reproduce all expectation values of interest. When decoherence is strong, however, these Monte Carlo methods can take a prohibitively long time to converge, as simulations become dominated by incoherent jumps that generate large numbers of distinct quantum trajectories that need to be averaged in order to accurately compute expectation values. Even with strong decoherence, dynamics are sometimes solvable through the cumulant expansionMeiser and Holland (2010) that neglects all $n$-body connected correlators for $n>2$. The growth of genuinely multi-body correlations, however, eventually causes the cumulant expansion to yield incorrect results with no clear signature of failure. In the absence of other means to compute correlators, it can therefore be difficult to identify the point at which correlators computed via the cumulant expansion can no longer be trusted.

In this work, we present an efficient method to compute short-time dynamics of collective spin systems with typical forms of Markovian decoherence. The only restriction on decoherence (beyond Markovianity) is that, like the coherent collective dynamics, it must act identically on all constituent particles. Our method is based on a formal solution to the equations of motion for Heisenberg operators, thereby bearing some resemblance to the Mori formalismMori (1965) and related workAnnett et al. (1994). Specifically, we expand a formal solution for a Heisenberg operator into a Taylor series whose truncation can yield negligible error at sufficiently short times. Evaluating the resulting expansion requires knowing the structure constants of a collective spin operator algebra; the calculation of these structure constants (in Appendices A–C) is one of the main technical results of this work, which we hope will empower both analytical and numerical studies of collective spin systems in the future. We benchmark our method against exact results from both analytical calculations and quantum trajectory Monte Carlo computations of spin squeezing in accessible parameter regimes, highlighting both advantages and limitations of the short-time expansion. Finally, we showcase applications of our method by computing quantities that are inaccessible to other numerical methods.

II Theory

In this section we provide the basic theory for our method to compute expectation values of collective spin operators, deferring lengthy derivations to the appendices. We consider a system of $N$ distinct spin-1/2 particles. Defining individual spin-1/2 operators $\hat{s}_{\alpha=\text{x},\text{y},\text{z}}\equiv\hat{\sigma}_{\alpha}/2$ and $\hat{s}_{\pm}\equiv\hat{s}_{\text{x}}\pm i\hat{s}_{\text{y}}=\hat{\sigma}_{\pm}$ with Pauli operators $\hat{\sigma}_{\alpha}$, we denote an operator that acts with $\hat{s}_{\alpha}$ on the spin indexed by $j$ and trivially (i.e. with the identity $\hat{\mathds{1}}$) on all other spins by $\hat{s}_{\alpha}^{(j)}$. We then define the collective spin operators $\hat{S}_{\alpha}\equiv\sum_{j=1}^{N}\hat{s}_{\alpha}^{(j)}$ for $\alpha\in\set{\text{x},\text{y},\text{z},+,-}$. Identifying the set $\set{\hat{\mathcal{S}}_{\bm{m}}}$ as a basis for all collective spin operators, with $\bm{m}\equiv\left(m_{+},m_{\text{z}},m_{-}\right)\in\mathbb{N}_{0}^{3}$ and $\hat{\mathcal{S}}_{\bm{m}}\equiv\hat{S}_{+}^{m_{+}}\hat{S}_{\text{z}}^{m_{\text{z}}}\hat{S}_{-}^{m_{-}}$, we can expand any collective spin operator $\hat{\mathcal{O}}$ in the form

[TABLE]

with scalar coefficients $\mathcal{O}_{\bm{m}}\in\mathbb{C}$. If $\hat{\mathcal{O}}$ is self-adjoint, for example, then $\mathcal{O}_{\bm{m}}^{*}=\mathcal{O}_{\bm{m}^{*}}$ with $\bm{m}^{*}\equiv\left(m_{-},m_{\text{z}},m_{+}\right)$. The corresponding Heisenberg operator is then $\hat{\mathcal{O}}\left(t\right)=\sum_{\bm{m}}\mathcal{O}_{\bm{m}}\left(t\right)\hat{\mathcal{S}}_{\bm{m}}+\hat{\mathcal{E}}_{\mathcal{O}}\left(t\right)$, with time-dependent coefficients $\mathcal{O}_{\bm{m}}\left(t\right)$ for time-independent Schrödinger operators $\hat{\mathcal{S}}_{\bm{m}}$, and mean-zero “noise” operators $\hat{\mathcal{E}}_{\mathcal{O}}\left(t\right)$ that result from interactions between the spin system and its environment, initially $\hat{\mathcal{E}}_{\mathcal{O}}\left(0\right)=0$. These noise operators will essentially play no role in the present work, but are necessary to include for a consistent formalism of Heisenberg operators in an open quantum system; see Appendix N for further discussion. The expectation values of Heisenberg operators evolve according to

[TABLE]

with a Heisenberg-picture time derivative operator $\check{T}=d/dt$ whose matrix elements $T_{\bm{m}\bm{n}}\in\mathbb{C}$ are defined by

[TABLE]

where $\left[X,Y\right]_{\pm}\equiv XY\pm YX$; $\hat{H}$ is the collective spin Hamiltonian; $\mathcal{J}$ is a set of jump operators with a corresponding decoherence rate $\gamma_{\mathcal{J}}$; and $\check{\mathcal{D}}$ is a Heisenberg-picture dissipator, or Lindblad superoperator, defined by

[TABLE]

Decoherence via uncorrelated decay of individual spins, for example, would be described by the set of jump operators \mathcal{J}_{-}\equiv\set{\hat{s}_{-}^{(j)}\mathrel{\mathop{\mathchar 58\relax}}j=1,2,\cdots,N}. The commutator in Eq. (3) can be computed by expanding the product $\hat{\mathcal{S}}_{\bm{\ell}}\hat{\mathcal{S}}_{\bm{m}}=\sum_{\bm{n}}f_{\bm{\ell}\bm{m}\bm{n}}\hat{\mathcal{S}}_{\bm{n}}$ with structure constants $f_{\bm{\ell}\bm{m}\bm{n}}\in\mathbb{R}$ that we work out in Appendices A–C, and the effects of decoherence from jump operators (i.e. elements of $\mathcal{J}$) of the form $\hat{g}^{(j)}=\sum_{\alpha}g_{\alpha}\hat{s}_{\alpha}^{(j)}$ and $\hat{G}=\sum_{\alpha}G_{\alpha}\hat{S}_{\alpha}$ are worked out in Appendices D–G. We consider these calculations to be some of the main technical contributions of this work, with potential applications beyond the short-time simulation method presented here. These ingredients are sufficient to compute matrix elements $T_{\bm{m}\bm{n}}$ of the time derivative operator $\check{T}$ in Eq. (3) in most cases of practical interest.

We note that particle loss is an important decoherence mechanism in many experimental realizations of collective spin modelsMa et al. (2011). In principle, a spin model has no notion of the particle annihilation operators that generate particle loss, and therefore cannot capture this effect directly. Nonetheless, for a system initially composed of $N$ particles, the effect of particle loss can be emulated with $O(1/N)$ error by the dissipator $\check{\mathcal{D}}_{\text{loss}}$ defined by $\check{\mathcal{D}}_{\text{loss}}\hat{\mathcal{S}}_{\bm{m}}=-\lvert\bm{m}\rvert\hat{\mathcal{S}}_{\bm{m}}$, where $\lvert\bm{m}\rvert\equiv\sum_{\alpha}m_{\alpha}$ (see Appendix H). Furthermore, the effect particle loss can be accounted for exactly by

1. introducing an additional index on spin operators, $\hat{\mathcal{S}}_{\bm{m}}\to\hat{\mathcal{S}}_{N\bm{m}}$, to keep track of different sectors of fixed particle number within a multi-particle Fock space, and

2. constructing jump operators that appropriately couple spin operators within different particle-number sectors.

We defer a detailed exact accounting of particle loss to future work.

The time derivative operator $\check{T}$ will generally couple spin operators $\hat{\mathcal{S}}_{\bm{n}}$ to spin operators $\hat{\mathcal{S}}_{\bm{m}}$ with higher “weight”, i.e. with $\lvert\bm{m}\rvert>\lvert\bm{n}\rvert$. The growth of operator weight signifies the growth of many-body correlations. Keeping track of this growth eventually becomes intractable, requiring us to truncate our equations of motion somehow. The simplest truncation strategy would be to take

[TABLE]

for some weight measure $w$, e.g. $w\left(\bm{m}\right)=\lvert\bm{m}\rvert$, and a high-weight cutoff $W$. The truncation in Eq. (5) closes the system of differential equations defined by Eq. (2), and allows us to solve it using standard numerical methods. Some initial conditions for this system of differential equations, namely expectation values of collective spin operators with respect to spin-polarized (Gaussian) states that are generally simple to prepare experimentally, are provided in Appendix I.

The truncation strategy in Eq. (5) has a few limitations:

1. simulating a system of differential equations for a large number of operators can be time-consuming,

2. the weight measure $w$ may need to be chosen carefully, as the optimal measure is generally system-dependent, and

3. simulation results can only be trusted up to the time at which the initial values of operators $\hat{\mathcal{S}}_{\bm{m}}$ with weight $w\left(\bm{m}\right)\geq W$ have a non-negligible contribution to expectation values of interest.

The last limitation in particular unavoidably applies in some form to any method tracking only a subset of all relevant operators. We therefore devise an alternate truncation strategy built around limitation 3.

We can formally expand Heisenberg operators $\hat{\mathcal{O}}\left(t\right)$ in a Taylor series about the time $t=0$ to write

[TABLE]

where the matrix elements $T^{k}_{\bm{m}\bm{n}}$ of the $k$-th time derivative operator $\check{T}^{k}$ are

[TABLE]

with $\delta_{\bm{m}\bm{n}}=1$ if $\bm{m}=\bm{n}$ and zero otherwise. For sufficiently short times, we can truncate the series in Eq. (6) by taking

[TABLE]

We refer to Eq. (10) as the truncated short-time (TST) expansion of Heisenberg operators. Note that when computing an expectation value $\braket{\hat{\mathcal{O}}\left(t\right)}$, the relation $\hat{\mathcal{S}}_{\bm{m}}^{\dagger}=\hat{\mathcal{S}}_{\bm{m}^{*}}$, which by Hermitian conjugation of Eq. (2) also implies that $T_{\bm{m}^{*}\bm{n}^{*}}=T_{\bm{m}\bm{n}}^{*}$, cuts both the number of initial-time expectation values $\braket{\hat{\mathcal{S}}_{\bm{m}}}$ and the number of matrix elements $T_{\bm{m}\bm{n}}$ that we may need to explicitly compute roughly in half.

Unlike the weight-based truncation in Eq. (5), the nonzero matrix elements $T^{k}_{\bm{m}\bm{n}}$ for $k=0,1,\cdots,M$ in Eq. (10) tell us which operators $\hat{\mathcal{S}}_{\bm{m}}$ are relevant for computing the expectation value $\braket{\hat{\mathcal{O}}\left(t\right)}$ to a fixed order $M$. The TST expansion thereby avoids the introduction of a weight measure $w$ that chooses which operators to keep track of, and trades the cost of solving a system of differential equations for the cost of computing expectation values $\braket{\hat{\mathcal{S}}_{\bm{m}}}$ and matrix elements $T_{\bm{m}\bm{n}}^{k}$. In all cases considered in this work, we find that the TST expansion is both faster to evaluate and provides accurate correlators $\braket{\hat{\mathcal{O}}\left(t\right)}$ until later times $t$ than the weight-based expansion in (5) with weight measure $w\left(\bm{m}\right)=\lvert\bm{m}\rvert$ and cutoff $W\approx M$. We therefore restrict the remainder of our discussions to the TST expansion in Eq. (10), and provide a pedagogical tutorial for computing correlators using the TST expansion in Appendix J.

Three primary considerations limit the maximum time $t$ to which we can accurately compute a correlator $\braket{\hat{S}_{\bm{n}}\left(t\right)}$ using the TST expansion. First, maintaining accuracy at larger times $t$ requires going to higher orders $M$ in the TST expansion. An order-$M$ TST expansion of the correlator $\braket{\hat{S}_{\bm{n}}\left(t\right)}$ can involve a significant fraction of operators $\hat{S}_{\bm{m}}$ with weight $\lvert\bm{m}\rvert\lesssim M$, which implies the need to compute $O\left(M^{3}\right)$ initial-time expectation values $\braket{\hat{S}_{\bm{m}}}$ and $O\left(M^{4}\right)$ matrix elements $T^{k}_{\bm{m}\bm{n}}$. In practice, with a straightforward implementation of the TST expansion we find that these requirements generally restrict $M\lesssim 50$ – $70$ with $8$ – $50$ gigabytes of random access memory (RAM). Second, individual terms at high orders of the TST expansion in Eq. (10) can grow excessively large, greatly amplifying any numerical errors and thereby spoiling cancellations that are necessary to arrive at a physical value of a correlator, i.e. with $\lvert\braket{\hat{S}_{\bm{n}}\left(t\right)}\rvert\lesssim S^{\lvert\bm{n}\rvert}$ (where $S\equiv N/2$). Finally, the TST expansion is essentially perturbative in the time $t$, which implies that its validity as a formal expansion eventually breaks down. Precisely characterizing the implications of these last two considerations for the TST expansion requires additional analysis that we defer to future work. An investigation of connections between the TST expansion and past work related to the Mori formalismMori (1965); Annett et al. (1994), for example, might answer questions about the breakdown and convergence of the TST expansion. As we show from benchmarks of the TST expansion in Section III, however, a detailed understanding of breakdown is not necessary to diagnose the breakdown time $t_{\text{break}}^{(M)}$ beyond which the TST expansion yields inaccurate results. Empirically, we find that going beyond order $M\approx 35$ yields no significant gains in all cases considered in this work.

III Spin squeezing, benchmarking, and breakdown

To benchmark our method for computing collective spin correlators, we consider three collective spin models known to generate spin-squeezed states: the one-axis twisting (OAT), two-axis twisting (TAT), and twist-and-turn (TNT) models described by the collective spin HamiltoniansMa et al. (2011)

[TABLE]

where we include a factor of $1/3$ in the TAT Hamiltonian because it naturally appears in realistic proposals to experimentally implement TATLiu et al. (2011); Huang et al. (2015a). For simplicity, we further fix $\Omega=\chi S$ (with $S\equiv N/2$ throughout this work) to the critical value known to maximize the entanglement generation rate of TNT in the large-$N$ limitMicheli et al. (2003); Sorelli et al. (2019).

Note that the OAT model is a special case of the zero-field Ising model, whose quantum dynamics admits an exact analytic solution even in the presence of decoherenceFoss-Feig et al. (2013). The approximate and numerics-oriented TST expansion is therefore an inappropriate tool for studying the OAT model, which will merely serve as an exactly solvable benchmark of our methods. Wherever applicable, we will provide exact results for the OAT model (see Appendix K, as well as the Supplementary Material of Ref. [14]).

The Hamiltonians in Eqs. (11)–(13) squeeze the initial product state $\ket{\text{X}}\propto\left(\ket{\uparrow}+\ket{\downarrow}\right)^{\otimes N}$ with $\hat{S}_{\text{x}}\ket{\text{X}}=S\ket{\text{X}}$. Our measure of spin squeezing is the directionally-unbiased Ramsey squeezing parameter determined by the maximal gain in resolution $\Delta\phi$ of a phase angle $\phi$ over that achieved by any spin-polarized product state (e.g. $\ket{\text{X}}$)Wineland et al. (1992); Ma et al. (2011),

[TABLE]

where $\bm{\hat{S}}\equiv\left(\hat{S}_{\text{x}},\hat{S}_{\text{y}},\hat{S}_{\text{z}}\right)$ is a collective spin operator-valued vector, the minimization is performed over all unit vectors $\bm{v}$ orthogonal to the mean spin vector $\braket{\bm{\hat{S}}}$, and for brevity we have suppressed the explicit time dependence of operators in Eq. (14). This squeezing parameter is entirely determined by one- and two-spin correlators of the form $\braket{\hat{S}_{\alpha}}$ and $\braket{\hat{S}_{\alpha}\hat{S}_{\beta}}$. For the unitary dynamics discussed in this work, these correlators are obtainable via exact simulations of quantum dynamics in the $\left(N+1\right)$-dimensional Dicke manifold of states $\set{\ket{S,m}}$ with net spin $S$ and spin projection $m$ onto the $z$ axis, i.e. with $\braket{S,m}{\bm{\hat{S}}^{2}}{S,m}=S\left(S+1\right)$ and $\braket{S,m}{\hat{S}_{\text{z}}}{S,m}=m$ for $m\in\set{-S,-S+1,\cdots,S}$. In the presence of single-spin or collective decoherence, meanwhile, these correlators are obtainable with the collective-spin quantum trajectory Monte Carlo method developed in ref. Zhang et al. (2018). In this work, these exact and quantum trajectory simulations will be used to benchmark the TST expansion in Eq. (10).

Figure 1 compares the squeezing parameter $\xi^{2}$ for $N=10^{4}$ spins initially in the state $\ket{\text{X}}$ evolved under the Hamiltonians in Eqs. (11)–(13), as computed via both benchmarking simulations and the TST expansion in Eq. (10) with $M=35$. Squeezing is shown for both unitary dynamics (Figure 1a), as well as non-unitary dynamics in the presence of spontaneous decay, excitation, and dephasing of individual spins at rates $\chi$ (Figure 1b), respectively described by the sets of jump operators $\mathcal{J}_{\alpha}\equiv\set{\hat{s}_{\alpha}^{(j)}}$ with corresponding decoherence rates $\gamma_{\alpha}=\chi$ for $\alpha\in\set{-,+,\text{z}}$. The results shown in Figure 1 were computed in a rotated basis with $\left(\hat{s}_{\text{z}},\hat{s}_{\text{x}}\right)\to\left(\hat{s}_{\text{x}},-\hat{s}_{\text{z}}\right)$ and $\ket{\text{X}}\to\ket{-\text{Z}}\equiv\ket{\downarrow}^{\otimes N}$, as well as appropriate transformations of the Hamiltonian and jump operators. The only effects of this rotation on the results presented in Figure 1 are to

1. reduce the time it takes to compute correlators $\braket{\hat{\mathcal{O}}\left(t\right)}$ with the TST expansion, and

2. prolong the time for which the TST expansion of TNT results agree with benchmarking simulations.

The speedup in a different basis occurs because for the initial state $\ket{-\text{Z}}$, all initial-time correlators $\braket{\hat{S}_{\bm{m}}}$ are zero unless $m_{+}=m_{-}=0$, and all non-zero correlators take $O\left(1\right)$ (i.e. constant in $N$) time to compute, rather than $O\left(N\right)$ time (see Appendix I). In total, the use of a rotated basis reduces the computation time of initial-time correlators from $O\left(M^{3}N\right)$ to $O\left(M\right)$. The reason for prolonged agreement of TNT results in a rotated basis is not entirely understood, and provides a clue into the precise mechanism by which the TST expansion breaks down (discussed below). We defer a detailed study of this breakdown to future work.

The main lesson from Figure 1 is that the TST expansion yields essentially exact results right up until a sudden and drastic departure that can be diagnosed by inspection. The breakdown of the TST expansion in Figure 1 induces an unphysical squeezing parameter $\xi^{2}<0$. In general, however, there is no fundamental relationship between the breakdown of the TST expansion and the conditions for a physical squeezing parameter $\xi^{2}$. A proper diagnosis of breakdown therefore requires inspection of the correlators $\braket{\hat{\mathcal{S}}_{\bm{n}}\left(t\right)}$ used to compute the squeezing parameter $\xi^{2}$, which upon breakdown will rapidly take unphysical values with $\lvert\braket{\hat{\mathcal{S}}_{\bm{n}}\left(t\right)}\rvert\gtrsim S^{\lvert\bm{n}\rvert}$ (see Appendix L for an example). The sudden and drastic departure from virtually exact results is consistent with the limitations of the TST expansion discussed at the end of Section II. Specifically, we identify three possible mechanisms for breakdown:

1. a rapid growth in the order $M$ necessary for the TST expansion to converge,

2. the growth of numerical errors in excessively large terms of the TST expansion, and

3. the formal breakdown of the perturbative expansion in the time $t$.

In all of these cases, a detailed cancellation eventually ceases to occur between large terms at high orders in the TST expansion. These large terms grow with the time $t$ raised to some large power (as high as $M$), and therefore rapidly yield wildly unphysical results. In contrast to other approximate methods such as the cumulant expansionMeiser and Holland (2010), the TST expansion can thus diagnose its own breakdown, which is an important feature when working in parameter regimes that are inaccessible via other means to compute correlators. Note that, due to the breakdown mechanisms of the TST expansion, going up through order $M=70$ does not significantly increase the breakdown time $t_{\text{break}}^{(M)}$ in Figure 1, and in some cases even shortens $t_{\text{break}}^{(M)}$.

Although the TST expansion breaks down at short times, it has two key advantages over the quantum trajectory Monte Carlo method to compute correlators in the presence of decoherence. First, computing spin correlators with the TST expansion is generally faster and requires fewer computing resources. The TST expansion results in Figure 1b, for example, take $\sim 10$ seconds to compute with a single CPU on modern computing hardware. The quantum trajectory Monte Carlo results in the same figure, meanwhile, take $\sim 10^{4}$ CPU hours to compute on similar hardware; the bulk of this time is spent performing sparse matrix-vector multiplication, leaving little room to further optimize serial runtime. Parallelization can reduce actual runtime of the Monte Carlo simulations to $\sim 10$ hours by running all trajectories at once, but at the cost of greatly increasing computing resource requirements. Though it may be possible to further speed up quantum trajectory Monte Carlo simulations by introducing new truncation schemes, any modifications

1. should be made carefully to ensure that simulations still yield correct results, and

2. are unlikely to bridge the orders of magnitude in computing resource requirements.

The second advantage of the TST expansion is the capability to compute spin correlators in strong-decoherence regimes of large systems that are entirely inaccessible to other methods. As an example, Figure 2 shows squeezing of $N=10^{4}$ spins initially in $\ket{\text{X}}$, undergoing spontaneous decay, excitation, and dephasing of individual spins at rates $\gamma_{-}=\gamma_{+}=\gamma_{\text{z}}=100\chi$. The system size in these simulations is too large for straightforward application of exact methods for open quantum systems. Quantum trajectory Monte Carlo simulations, meanwhile, take a prohibitively long time to converge with such strong decoherence due to the multiplicity of quantum trajectories that require averaging.

The results in Figure 2 show that the TNT model can generate more squeezing than the OAT or TAT models in the presence of strong decoherence. The better performance of TNT is in part a consequence of the fact that TNT initially generates squeezing at a faster rate than OAT or TAT, thereby allowing it to produce more squeezing before the degrading effects of decoherence kick in. We corroborate this finding with quantum trajectory simulations of a smaller system in Appendix M. Strong-decoherence computations of the sort used for Figure 2 put lower bounds on theoretically achievable spin squeezing via TAT with decoherence in Ref. [48], exemplifying a concrete and practical application of the TST expansion and the collective-spin structure constants calculated in this work.

IV Two-time correlation functions and out-of-time-ordered

correlators

As a final example of collective-spin physics that is numerically accessible via the TST expansion of Heisenberg operators, we consider the calculation of two-time correlation functions and out-of-time-ordered correlators (OTOCs). In particular, we consider the effect of decoherence on short-time behavior of the two-time connected correlator

[TABLE]

and the expectation value of a squared commutator,

[TABLE]

in the context of the squeezing models in Section III. The subscript on $\braket{\cdot}_{\text{nn}}$ in Eq. (16) stands for “no noise”, and denotes a correlator computed without the noise contributions $\hat{\mathcal{E}}_{\mathcal{O}}\left(t\right)$ to Heisenberg operators $\hat{\mathcal{O}}\left(t\right)$. While linear contributions from noise operators as e.g. in Eq. (15) always vanish under Markovian decoherence (see Appendix N), quadratic contributions that would otherwise appear in Eq. (16) generally do notBlocher and Mølmer (2019). Determining the effect of these noise terms generally requires making additional assumptions about the environment, which would be a digression for the purposes of the present work. We therefore exclude these noise terms in (16) in order to keep our discussion simple and general; see Ref. [65] for more detailed discussions of noise terms and the quantum regression theorem underlying the calculation of multi-time correlators.

In an equilibrium setting, correlation functions similar to that in Eq. (15) contain information about the linear response of Heisenberg operators to perturbations of a system; in a non-equilibrium setting, they contribute to short-time linear response (see Appendix O). Similar correlators have made appearances as order parameters for diagnosing time-crystalline phases of matterTucker et al. (2018). Squared commutators such as that in Eq. (16), meanwhile, are commonly examined for signatures of quantum chaos and information scramblingMaldacena et al. (2016); Swingle (2018); García-Mata et al. (2018). In typical scenarios, such squared commutators initially vanish by construction through a choice of spatially separated operators. Collective spin systems, however, have no intrinsic notion of locality or spatial separation. In our case, therefore, with the choice of initial state $\ket{\text{X}}\propto\left(\ket{\uparrow}+\ket{\downarrow}\right)^{\otimes N}$ we merely have $D\left(0\right)\sim 1/N$.

Figure 3 shows the behavior of $C\left(t\right)$ and $D\left(t\right)$ for $N=10^{4}$ spins, initially in the state $\ket{\text{X}}$, evolving under the squeezing Hamiltonians in Eqs. (11)–(13) both with and without spontaneous decay, excitation, and dephasing of individual spins at rates $\gamma_{-}=\gamma_{+}=\gamma_{\text{z}}=100\chi$. In the case of unitary evolution under OAT, we find that to an excellent approximation $\lvert C\left(t\right)\rvert$ takes the functional form $f\left(t\right)=f\left(0\right)+aN\chi t+\left(bN\chi t\right)^{2}$ with $a\sim b\sim 1$, and with a virtually perfect fit $D\left(t\right)=D\left(0\right)+\left(\left[N+1\right]\chi t\right)^{2}$. For unitary evolution under TAT and TNT, we find that to an excellent approximation both $\lvert C\left(t\right)\rvert$ and $D\left(t\right)$ take the functional form $f\left(t\right)=f\left(0\right)+a\left[\exp\left(bN\chi t\right)-1\right]$ with $a\sim b\sim 1$. As may be expected, the growth of $C\left(t\right)$ and $D\left(t\right)$ is generally suppressed by decoherence. Figure 3 serves as an example for the type of behavior that is accessible at short times with the TST expansion. These examples are straightforward to extend to equilibrium settings and spin-boson systems.

V Conclusions

We have presented an efficient method for computing correlators at short times in collective spin systems. This method is based on truncating a short-time expansion of Heisenberg operators, and can access correlators on time scales that are relevant to metrological applications such as spin squeezing. In order to evaluate the truncated short-time (TST) expansion of Heisenberg operators, we have computed the structure constants of a collective spin operator algebra, which we hope will empower future analytical and numerical studies of collective spin systems. Even though we considered only non-equilibrium spin-squeezing processes in this work, our method can be applied directly in an equilibrium setting, and is straightforward to generalize to systems such as trapped ions and optical cavities with collective spin-boson interactions. In such contexts, our method may be used to benchmark the short-time effects of decoherence, or study the onset of quantum chaos and information scrambling.

Acknowledgements.

We acknowledge helpful discussions with Robert Lewis-Swan, Kris Tucker, and Colin Kennedy; as well as some technical contributions from Diego Barberena. This work is supported by the Air Force Office of Scientific Research (AFOSR) grant FA9550-18-1-0319; the Defense Advanced Research Projects Agency (DARPA) and Army Research Office (ARO) grant W911NF-16-1-0576 and W911NF-19-1-0210; the National Science Foundation (NSF) grant PHY-1820885; JILA-NSF grant PFC-1734006; and the National Institute of Standards and Technology (NIST).

Appendix A Basic spin operator identities

The appendices in this work make ubiquitous use of various spin operator identities; we collect and derive some basic identities here for reference. Note that despite the working definition of collective spin operators from $S_{\alpha}=\sum_{j}s_{\alpha}^{(j)}$, the identities we will derive involving only collective spin operators apply just as well to large-spin operators that cannot be expressed as the sum of individual spin-1/2 operators. The elementary commutation relations between spin operators are, with ${\bar{\mu}}\equiv-\mu\in\set{+1,-1}$ for brevity,

[TABLE]

These relations can be used to inductively compute identities involving powers of collective spin operators. By pushing through one spin operator at a time, we can find

[TABLE]

and

[TABLE]

where we will generally find it nicer to express results in terms of $\mu s_{\text{z}}^{(j)}$ and $\mu S_{\text{z}}$ rather than $s_{\text{z}}^{(j)}$ and $S_{\text{z}}$. Summing over the single-spin index $j$ in both of the cases above gives us the purely collective-spin versions of these identities:

[TABLE]

where we can repeat the process of pushing through individual $S_{\text{z}}$ operators $\ell$ times to get

[TABLE]

Multiplying (22) through by $\left(\mu\nu\right)^{\ell}$ (for $\nu\in\set{+1,-1}$) and taking its Hermitian conjugate, we can say that more generally

[TABLE]

Finding commutation relations between powers of transverse spin operators, i.e. $S_{\mu}$ and $S_{\bar{\mu}}$, turns out to be considerably more difficult than the cases we have worked out thus far. We therefore save this work for Appendix B.

Appendix B Commutation relations between powers of transverse spin

operators

To find commutation relations between powers of transverse collective spin operators, we first compute

[TABLE]

While (26) gives us the commutator $\left[S_{\mu}^{m},s_{\bar{\mu}}^{(j)}\right]_{-}$, we would like to enforce an ordering on products of spin operators, which will ensure that we only keep track of operators that are linearly independent. We choose (for now) to impose an ordering with all $s_{\bar{\mu}}^{(j)}$ operators on the left, and all $s_{\text{z}}^{(j)}$ operators on the right. Such an ordering will prove convenient for the calculations in this section111In retrospect, it may have been nicer to push all $s_{\mu}^{(j)}$ operators to the right throughout these calculations, due to the enhanced symmetry that expressions would have with respect to Hermitian conjugation. In any case, we provide the final result of this section in both ordering conventions, and therefore feel no need to reproduce these calculations with a different ordering of spin operators.. This choice of ordering compels us to expand

[TABLE]

which implies

[TABLE]

and in turn

[TABLE]

As the next logical step, we take on the task of computing

[TABLE]

which implies

[TABLE]

We now need to rearrange the operators in $C_{mn;\mu}$ into a standard order, which means pushing all $S_{\text{z}}$ operators to the right and, for the purposes of this calculation, all $S_{\bar{\mu}}$ operators to the left. We begin by pushing $S_{\mu}^{k}$ to the left of $S_{\text{z}}$, which takes $2\mu S_{\text{z}}\to 2\mu S_{\text{z}}+2k$, and then push $S_{\mu}^{m-k-1}$ to the right of $S_{\bar{\mu}}^{n-1}$, giving us

[TABLE]

where we have dropped the last ($k=m-1$) term in the remaining sum because $C_{m-k-1,n-1;\mu}=0$ if $k=m-1$, and

[TABLE]

To our despair, we have arrived in (34) at a recursive formula for $C_{mn;\mu}$. Furthermore, we have not even managed to order all spin operators, as $C_{m-k-1,n-1;\mu}$ contains $S_{\text{z}}$ operators that are to the left of $S_{\mu}^{k}$. To sort all spin operators once and for all, we define

[TABLE]

which we can expand as

[TABLE]

and

[TABLE]

While the resulting expression in (42) strongly resembles that in (34), there is one crucial difference: all spin operators in (42) have been sorted into a standard order. We can now repeatedly substitute $C_{mn;\mu}^{(k)}$ into itself, each time decreasing $m$ and $n$ by 1, until one of $m$ or $n$ reaches zero. Such repeated substitution yields the expansion

[TABLE]

where the first two terms in the sum over $p$ are

[TABLE]

and more generally for $p>1$,

[TABLE]

In principle, the expressions in (35), (38), (43), and (46) suffice to evaluate the commutator $\left[S_{\mu}^{m},S_{\bar{\mu}}^{n}\right]_{-}=C_{mn;\mu}$, but this result is – put lightly – quite a mess: the expression for $E_{mn;\mu}^{(p)}$ in (46) involves a sum over $p$ mutually dependent intermediate variables, each term of which additionally contains a product of $p$ factors. We therefore devote the rest of this section to simplifying our result for the commutator $\left[S_{\mu}^{m},S_{\bar{\mu}}^{n}\right]_{-}$.

Observing that in (46) we always have $0\leq k_{1}\leq k_{2}\leq\cdots\leq k_{p}\leq m-p-1$, we can rearrange the order of the sums and relabel $k_{p}\to\ell$ to get

[TABLE]

where for shorthand we define

[TABLE]

We now further define

[TABLE]

and evaluate sums successively over $k_{p-1},k_{p-2},\cdots,k_{1}$, finding

[TABLE]

Substitution of this result together with $D_{m-p,n-p;\mu}^{(\ell)}$ using (38) into (47) then gives us

[TABLE]

with

[TABLE]

We can further simplify

[TABLE]

where ${p\brack q}$ is an unsigned Stirling number of the first kind, and finally

[TABLE]

Putting everything together, we finally have

[TABLE]

with

[TABLE]

where in this final form $E_{mn;\mu}^{(0)}=D_{mn;\mu}$, which together with the expansion for $C_{mn;\mu}$ in (43) implies that

[TABLE]

and

[TABLE]

If we wish to order products of collective spin operators with $S_{\text{z}}$ in between $S_{\bar{\mu}}$ and $S_{\mu}$, then

[TABLE]

where

[TABLE]

with

[TABLE]

Here ${p\brack s}$ is an unsigned Stirling number of the first kind.

Appendix C Product of arbitrary ordered collective spin operators

The most general product of collective spin operators that we need to compute is

[TABLE]

where

[TABLE]

is defined in terms of $Z_{r\ell;\mu}^{(k)}$ and $\zeta_{r\ell}^{ka}$ as respectively given in (65) and (66). The (anti-)commutator of two ordered products of collective spin operators is then simply

[TABLE]

Appendix D Sandwich identities for single-spin decoherence calculations

In this section we derive several identities that will be necessary for computing the effects of single-spin decoherence on ordered products of collective spin operators, i.e. on operators of the form $S_{\mu}^{\ell}\left(\mu S_{\text{z}}\right)^{m}S_{\bar{\mu}}^{n}$. These identities all involve sandwiching a collective spin operator between operators that act on individual spins only, and summing over all individual spin indices. Our general strategy will be to use commutation relations to push single-spin operators together, and then evaluate the sum to arrive at an expression involving only collective spin operators.

We first compute sums of single-spin operators sandwiching $\left(\mu S_{\text{z}}\right)^{m}$, when necessary making use of the identity in (19). The unique cases up to Hermitian conjugation are, for $S\equiv N/2$ and $\mu,\nu\in\set{+1,-1}$,

[TABLE]

We are now equipped to derive similar identities for more general collective spin operators. Making heavy use of identities (20) and (29) to push single-spin operators through transverse collective-spin operators, we again work through all combinations that are unique up to Hermitian conjugation, finding

[TABLE]

Appendix E Uncorrelated, permutationally-symmetric single-spin

decoherence

In this section we work out the effects of permutationally-symmetric decoherence of individual spins on collective spin operators of the form $S_{\mu}^{\ell}\left(\mu S_{\text{z}}\right)^{m}S_{\bar{\mu}}^{n}$. For compactness, we define

[TABLE]

where $g$ is an operator that acts on a single spin, $g^{(j)}$ is an operator that acts with $g$ on spin $j$ and trivially on all other spins, and $N$ is the total number of spins.

E.1 Decay-type decoherence

The effect of decoherence via uncorrelated decay ($\mu=-1$) or excitation ($\mu=1$) of individual spins is described by

[TABLE]

In order to determine the effect of this decoherence on general collective spin operators, we expand the anti-commutator

[TABLE]

which implies, using (79),

[TABLE]

and, using (80),

[TABLE]

Decoherence via jump operators $s_{\bar{\mu}}^{(j)}$ only couples operators $S_{\mu}^{\ell}\left(\mu S_{\text{z}}\right)^{m}S_{\bar{\mu}}^{n}$ to operators $S_{\mu}^{\ell}\left(\mu S_{\text{z}}\right)^{m^{\prime}}S_{\bar{\mu}}^{n}$ with $m^{\prime}\leq m$. Decoherence via jump operators $s_{\mu}^{(j)}$, meanwhile, makes operators $S_{\mu}^{\ell}\left(\mu S_{\text{z}}\right)^{m}S_{\bar{\mu}}^{n}$ “grow” in $m$ through the last term in (85), although the sum $\ell+m+n$ does not grow.

E.2 Dephasing

The effect of decoherence via single-spin dephasing is described by

[TABLE]

From (75), we then have

[TABLE]

Decoherence via single-spin dephasing makes operators $S_{\mu}^{\ell}\left(\mu S_{\text{z}}\right)^{m}S_{\bar{\mu}}^{n}$ “grow” in $m$, although the sum $\ell+m+n$ does not grow.

E.3 The general case

The most general type of single-spin decoherence is described by

[TABLE]

To simplify (88), we expand

[TABLE]

and

[TABLE]

which implies

[TABLE]

In order to compute the effect of this decoherence on general collective spin operators, we expand the anti-commutator

[TABLE]

Recognizing a resemblance between terms in (92) and (76), we collect terms to simplify

[TABLE]

and likewise

[TABLE]

with

[TABLE]

Defining for completion

[TABLE]

and

[TABLE]

we finally have

[TABLE]

Note that the sum $\ell+m+n$ for operators $S_{\mu}^{\ell}\left(\mu S_{\text{z}}\right)^{m}S_{\bar{\mu}}^{n}$ does not grow under this type of decoherence.

Appendix F Sandwich identities for collective-spin decoherence

calculations

In analogy with the work in Appendix D, in this section we work out sandwich identities necessary for collective-spin decoherence calculations. The simplest cases are

[TABLE]

With a bit more work, we can also find

[TABLE]

which implies

[TABLE]

Finally, we compute

[TABLE]

where

[TABLE]

so

[TABLE]

Appendix G Collective spin decoherence

In this section we work out the effects of collective decoherence on general collective spin operators. For shorthand, we define

[TABLE]

where $G$ is a collective spin jump operator.

G.1 Decay-type decoherence and dephasing

Making use of the results in Appendix F, we find that the effects of collective decay-type decoherence on general collective spin operators are given by

[TABLE]

and

[TABLE]

Similarly, the effect of collective dephasing is given by

[TABLE]

G.2 The general case

More generally, we consider jump operators of the form

[TABLE]

whose decoherence effects are determined by

[TABLE]

and

[TABLE]

which implies

[TABLE]

In order to compute the effect of this decoherence on general collective spin operators, we expand the anti-commutators

[TABLE]

Collecting terms and defining

[TABLE]

we then have

[TABLE]

Note that the sum $\ell+m+n$ for operators $S_{\mu}^{\ell}\left(\mu S_{\text{z}}\right)^{m}S_{\bar{\mu}}^{n}$ grows by one if $G_{\mu}\neq 0$ or $G_{\bar{\mu}}\neq 0$, and does not grow otherwise.

Appendix H Emulating particle loss in a spin model

Here we discuss the details of emulating particle loss with $O(1/N)$ error, where $N$ is the initial number of particles in a system that we wish to describe with a spin model. Starting with the full algebra of creation and annihilation operators (whether bosonic or fermionic) in a system, spin models are typically implemented by identifying a subalgebra of relevant “spin” operators that satisfy appropriate commutation relations. Two-state particles on a lattice, for example, are described by annihilation operators $c_{js}$ indexed by a lattice site $j\in\mathbb{Z}$ and an internal state index $s\in\set{\uparrow,\downarrow}$, enabling the straightforward construction of spin operators

[TABLE]

which satisfy the same commutation relations as the standard Pauli operators. These spin operators can be more compactly defined in the form

[TABLE]

where $\sigma_{\alpha}$ for \alpha\in\set{\text{x},\text{y},\text{z},\mathds{1}} is a Pauli operator, with $\sigma_{\mathds{1}}\equiv\mathds{1}$; and $\braket{r}{\sigma_{\alpha}}{s}$ denotes a matrix element of $\sigma_{\alpha}$. This construction exemplifies how the set of jump operators $\mathcal{J}_{\text{loss}}^{\text{bare}}\equiv\set{c_{js}}$ that generate particle loss cannot be constructed from spin operators, which are generally bilinear in particle creation or annihilation operators. When working on the level of a spin model, therefore, we can at best only emulate the effect of particle loss by some indirect means.

To understand the effect of particle loss on collective spin operators, we first define a single multi-body spin operator addressing sites $\bm{j}=\set{j_{1},j_{2},\cdots,j_{\lvert\bm{j}\rvert}}$,

[TABLE]

and expand

[TABLE]

In order to have an actual spin model, fermionic statistics or energetic considerations must forbid multiple occupation of individual lattice sites. In that case, the on-site four-point product $c_{ks}^{\dagger}c_{kq}^{\dagger}c_{kr}c_{ks}=0$ vanishes, and

[TABLE]

Up to $O(1/N)$ corrections, a collective spin operator $\mathcal{S}_{\bm{m}}$ essentially consists of $\lvert\bm{m}\rvert$-body operators of the form $\sigma_{\bm{j}\bm{\alpha}}$ with $\lvert\bm{j}\rvert=\lvert\bm{m}\rvert$, which implies that the dissipator $\mathcal{D}_{\text{loss}}$ defined by $\mathcal{D}_{\text{loss}}\mathcal{S}_{\bm{m}}=-\lvert\bm{m}\rvert\mathcal{S}_{\bm{m}}$ describes particle loss with $O(1/N)$ error. We note that the dissipator $\mathcal{D}_{\text{loss}}$ is essentially the depolarizing channel, i.e. $\mathcal{D}_{\text{loss}}=\mathcal{D}\left(\mathcal{J}_{\text{loss}}\right)$ for $\mathcal{J}_{\text{loss}}=\set{s_{\alpha}^{(j)}}$ with $\alpha\in\set{\text{x},\text{y},\text{z}}$ and $j\in\set{1,2,\cdots,N}$. A direct implementation of $\mathcal{D}_{\text{loss}}$ with $\mathcal{D}_{\text{loss}}\mathcal{S}_{\bm{m}}=-\lvert\bm{m}\rvert\mathcal{S}_{\bm{m}}$, however, is much more efficient than evaluating the depolarizing channel $\mathcal{D}\left(\mathcal{J}_{\text{loss}}\right)$ with the ingredients in Appendices D and E.

Appendix I Initial conditions

Here we compute the expectation values of collective spin operators with respect to spin-polarized (also Gaussian, or spin-coherent) states. These states are parameterized by polar and azimuthal angles $\theta\in[0,\pi)$, $\phi\in[0,2\pi)$, and lie within the Dicke manifold spanned by states $\ket{k}\propto S_{+}^{S+k}\ket{\downarrow}^{\otimes N}$ with $S\equiv N/2$ and $S_{\text{z}}\ket{k}=k\ket{k}$:

[TABLE]

We can likewise expand, within the Dicke manifold,

[TABLE]

where ${\bar{\mu}}\equiv-\mu\in\set{+1,-1}$ and

[TABLE]

which implies

[TABLE]

This expansion allows us to compute the expectation value

[TABLE]

where

[TABLE]

[TABLE]

Defining the states

[TABLE]

some particular expectation values of interest are

[TABLE]

and

[TABLE]

Appendix J Computing correlators with the truncated short-time (TST)

expansion

Here we provide a pedagogical tutorial for computing correlators using the truncated short-time TST expansion. For concreteness, we nominally consider $N$ spins evolving under the one-axis twisting (OAT) Hamiltonian

[TABLE]

additionally subject to spontaneous single-spin decay at rate $\gamma_{-}$, with jump operators \mathcal{J}_{-}=\set{s_{-}^{(j)}\mathrel{\mathop{\mathchar 58\relax}}j=1,2,\cdots,N}. The equation of motion for a Heisenberg operator $\left(S_{+}^{\ell}S_{\text{z}}^{m}S_{-}^{n}\right)\left(t\right)$ is

[TABLE]

where we have suppressed the explicit time dependence of operators for brevity. Using the results in appendices C and E.1 respectively to evaluate the commutator $\left[S_{\text{z}}^{2},S_{+}^{\ell}S_{\text{z}}^{m}S_{-}^{n}\right]_{-}$ and dissipator $\mathcal{D}\left(\mathcal{J}_{-}\right)\left(S_{+}^{\ell}S_{\text{z}}^{m}S_{-}^{n}\right)$ in (148), we can expand

[TABLE]

In practice, we do not want to keep track of such an expansion by hand, especially in the case of e.g. the two-axis twisting (TAT) and twist-and-turn (TNT) models with more general types of decoherence, for which the analogue of (149) may take several lines just to write out in full. Defining the operators $\mathcal{S}_{\bm{m}}\equiv S_{+}^{m_{+}}S_{\text{z}}^{m_{\text{z}}}S_{-}^{m_{-}}$ with $\bm{m}\equiv\left(m_{+},m_{\text{z}},m_{-}\right)$ for shorthand, we note that the vector space spanned by $\set{\mathcal{S}_{\bm{m}}}$ is closed under time evolution. We therefore expand

[TABLE]

where $T$ is a superoperator that generates time evolution for Heisenberg operators. In the present example, the matrix elements $T_{\bm{m}\bm{n}}\in\mathbb{C}$ of $T$ are defined by (149) and (150). For any Hamiltonian $H$ with decoherence characterized by sets of jump operators $\mathcal{J}$ and decoherence rates $\gamma_{\mathcal{J}}$, the matrix elements $T_{\bm{m}\bm{n}}$ are more generally defined by

[TABLE]

The results in Appendices C, E, and G can be used to write model-agnostic codes that compute matrix elements $T_{\bm{m}\bm{n}}$, taking a particular Hamiltonian $H$ and decoherence processes $\set{\left(\mathcal{J},\gamma_{\mathcal{J}}\right)}$ as inputs.

In order to compute a quantity such as spin squeezing, we need to compute correlators of the form $\braket{\mathcal{S}_{\bm{n}}\left(t\right)}$, where for clarity we will re-introduce the explicit time dependence of Heisenberg operators $\mathcal{S}_{\bm{n}}\left(t\right)$. The order-$M$ truncated short-time (TST) expansion takes

[TABLE]

where $T^{k}_{\bm{m}\bm{n}}$ are matrix elements of the $k$-th time derivative operator $T^{k}$, given by

[TABLE]

Matrix elements $T^{k}_{\bm{m}\bm{n}}$ and initial-time expectation values $\braket{\mathcal{S}_{\bm{m}}\left(0\right)}$ are thus computed as needed for any particular correlator $\braket{\mathcal{S}_{\bm{n}}\left(t\right)}$ of interest, and combined according to (152). Note that initial-time expectation values $\braket{\mathcal{S}_{\bm{m}}\left(0\right)}$ are an input to the TST expansion, and need to be computed separately for any initial state of interest; expectation values with respect to spin-polarized (Gaussian) states are provided in Appendix I. In practice, we further collect terms in (152) to write

[TABLE]

where $c_{\bm{n}k}$ are time-independent coefficients for the expansion of $\braket{\mathcal{S}_{\bm{n}}\left(t\right)}$. After computing the coefficients $c_{\bm{n}k}$, there is only negligible computational overhead to compute the correlator $\braket{\mathcal{S}_{\bm{n}}\left(t\right)}$ for any time $t$.

Appendix K Analytical results for the one-axis twisting model

The one-axis twisting (OAT) Hamiltonian for $N$ spin-1/2 particles takes the form

[TABLE]

where $\sigma_{\text{z}}^{(j)}$ represents a Pauli-$z$ operator acting on spin $j$. This model is a special case of the zero-field Ising Hamiltonian previously solved in Ref. [64] via exact, analytical treatment of the quantum trajectory Monte Carlo method for computing expectation values. The solution therein accounts for coherent evolution in addition to decoherence via uncorrelated single-spin decay, excitation, and dephasing respectively at rates $\gamma_{-}$, $\gamma_{+}$, and $\gamma_{\text{z}}$ (denoted by $\Gamma_{\text{ud}}$, $\Gamma_{\text{du}}$, and $\Gamma_{\text{el}}$ in Ref. [64]). Letting $S\equiv N/2$ and $\mu,\nu\in\set{+1,-1}$, we adapt expectation values computed in Ref. [64] for the initial state $\ket{\text{X}}\propto\left(\ket{\uparrow}+\ket{\downarrow}\right)^{\otimes N}$ with $S_{\text{x}}\ket{\text{X}}=S\ket{\text{X}}$ evolving under $H_{\text{OAT}}$, finding

[TABLE]

where

[TABLE]

for

[TABLE]

In order to compute spin squeezing as measured by the Ramsey squeezing parameter $\xi^{2}$ defined in (14), we additionally need analytical expressions for $\braket{S_{\text{z}}\left(t\right)}$ and $\braket{S_{\text{z}}^{2}\left(t\right)}$. As these operators commute with both the OAT Hamiltonian and the single-spin operators $\sigma_{\text{z}}^{(j)}$, their evolution is governed entirely by decay-type decoherence (see Appendix E.1), which means

[TABLE]

The initial conditions $\braket{S_{\text{z}}\left(0\right)}=0$ and $\braket{S_{\text{z}}^{2}\left(0\right)}=S/2$ then imply

[TABLE]

With appropriate assumptions about the relevant sources of decoherence, the expectation values in (156)–(158) and (163) are sufficient to compute the spin squeezing parameter $\xi^{2}$ in (14) at any time throughout evolution of the initial state $\ket{\text{X}}$ under $H_{\text{OAT}}$.

Appendix L Diagnosing breakdown of the TST expansion

In Figure 1 of the main text, the TST expansion provided nearly exact results for squeezing until a sudden departure that quickly resulted in an unphysical squeezing parameter, $\xi^{2}<0$. In general, however, there is no fundamental relationship between the breakdown of the TST expansion and the conditions for a physical squeezing parameter $\xi^{2}$. A proper diagnosis of breakdown therefore requires inspection of the correlators $\braket{\mathcal{S}_{\bm{n}}\left(t\right)}$ used to compute the squeezing parameter $\xi^{2}$, which upon breakdown will rapidly take unphysical values with $\lvert\braket{\mathcal{S}_{\bm{n}}\left(t\right)}\rvert\gtrsim S^{\lvert\bm{n}\rvert}$. As an example, Figure 4 shows the squeezing parameter $\xi^{2}$ throughout decoherence-free evolution of $N=100$ spins initially in the state $\ket{\text{X}}$. In this example, the squeezing computed by the TST expansion for the TAT model diverges from the exact answer without an immediate and obvious signature of breakdown. Nonetheless, breakdown can still be diagnosed by inspection of individual correlators, as shown in Figure 5, where we plot $\text{Im}\braket{S_{+}^{2}}$ as a function of time for $N=100$ spins evolving under the TAT without decoherence. Figure 5 shows that breakdown clearly occurs around $N\chi t\lesssim 7$, when the correlator $\braket{S_{+}^{2}}$ begins to diverge to values $\gtrsim S^{2}=\left(N/2\right)^{2}=2500$ in magnitude. A joint inspection of figures 4 and 5 suffice to trace the anomalous behavior of $\xi^{2}$ from $N\chi t\approx 7$ back to $N\chi t\approx 6$, when it first took a sudden turn before becoming unphysical at $N\chi t\approx 8$.

Appendix M Spin squeezing with strong decoherence

Here we provide supplementary evidence of our finding in Section III that the TNT model can produce more squeezing than the OAT or TAT models in the presence of strong decoherence. To this end, Figure 6 shows the minimal squeezing parameter $\xi^{2}_{\text{min}}$ achievable with $N=100$ spins through the OAT, TAT, and TNT models as a function of the rate $\gamma_{0}$ at which individual spins undergo spontaneous decay, excitation, and dephasing. These results were computed with quantum trajectory simulations, with $10^{3}$ trajectories per data point. While the OAT and TAT models produce more squeezing than the TNT model with weak decoherence, this squeezing falls off faster with an increasing decoherence rate $\gamma_{0}$. The relative robustness of TNT is in part a consequence of the fact that TNT initially generates squeezing at a faster rate than OAT or TAT, thereby allowing it to produce more squeezing before the degrading effects of decoherence kick in.

Appendix N Heisenberg operators in open quantum systems

Here we explain the origin and character of the mean-zero “noise” operators $\mathcal{E}_{\mathcal{O}}\left(t\right)$ that appear in the expansion of a Heisenberg operator $\mathcal{O}\left(t\right)=\sum_{\bm{m}}\mathcal{O}_{\bm{m}}\left(t\right)\mathcal{S}_{\bm{m}}+\mathcal{E}_{\mathcal{O}}\left(t\right)$ with time-dependent coefficients $\mathcal{O}_{\bm{m}}\left(t\right)$ for time-independent Schrödinger operators $\mathcal{S}_{\bm{m}}$. Our discussion should clarify why noise operators play no role in our calculation of expectation values of the form $\braket{\mathcal{O}\left(t\right)}$ and $\braket{\mathcal{O}\mathcal{Q}\left(t\right)}$, despite the fact that noise operators generally do need to be considered in the calculation of more general multi-time correlators in open quantum systemsBlocher and Mølmer (2019).

In any closed quantum system with initial state $\rho$ and propagator $U\left(t\right)$, such that the state at time $t$ is $\rho\left(t\right)\equiv U\left(t\right)\rho U^{\dagger}\left(t\right)$, time-dependent Heisenberg operators $\mathcal{O}\left(t\right)$ are uniquely defined from time-independent Schrödinger operators $\mathcal{O}$ by

[TABLE]

Enforcing (164) for arbitrary initial states $\rho$ forces $\mathcal{O}\left(t\right)=U^{\dagger}\left(t\right)\mathcal{O}U\left(t\right)$. In an open quantum system, however, the definition of a Heisenberg operator is not so straightforward. Open systems can often be understood as subsystems of a larger closed system. Consider therefore an open system $S$ with environment $E$, a joint initial state $\rho_{SE}$, and propagator $U_{SE}\left(t\right)$. The reduced state $\rho_{S}\left(t\right)$ of $S$ at time $t$ is

[TABLE]

where $\rho_{S}\equiv\rho_{S}\left(0\right)$ is a time-independent state of $S$ in the Heisenberg picture, $\underline{S}$ denotes the space of operators on $S$, and the quantum channel $\mathcal{U}_{\underline{S}}\left(t\right)$ has the decompositionRivas and Huelga (2012)

[TABLE]

with ordinary operators $\mathcal{U}_{S}^{(j)}\left(t\right)$ on $S$. We can therefore expand

[TABLE]

where $\mathcal{U}_{\underline{S}}^{\dagger}\left(t\right)$ is the adjoint map of $\mathcal{U}_{\underline{S}}\left(t\right)$ (with respect to a trace inner product between operators on $S$), and we define the time-dependent operator

[TABLE]

We thus find that substituting $\underline{\mathcal{O}_{S}}\left(t\right)$ in place of $\mathcal{O}_{S}\left(t\right)$ suffices for the calculation of correlators $\braket{\mathcal{O}_{S}\left(t\right)}$, thereby accounting for the validity of the equation of motion in (2). As we show below, this substitution also suffices for the calculation of two-time correlators of the form $\braket{\mathcal{O}_{S}\mathcal{Q}_{S}\left(t\right)}$ when the environment $E$ is Markovian.

The problem with defining Heisenberg operators $\mathcal{O}_{S}\left(t\right)$ by $\underline{\mathcal{O}_{S}}\left(t\right)$ only becomes evident when considering products of Heisenberg operators. One would like for the product of two Heisenberg operators $\mathcal{O}_{S}\left(t\right)$ and $\mathcal{Q}_{S}\left(t\right)$ to satisfy $\mathcal{O}_{S}\left(t\right)\mathcal{Q}_{S}\left(t\right)=\left(\mathcal{O}_{S}\mathcal{Q}_{S}\right)\left(t\right)$. This intuition can be formalized by observing that

[TABLE]

where $\mathds{1}_{E}$ is the identity operator on $E$, expectation values of Heisenberg operators on system $A\in\set{S,E,SE}$ are taken with respect to the state $\rho_{A}$, and

[TABLE]

By expanding Heisenberg operators similarly to (169) and (170), we then find

[TABLE]

The expression in (168), however, makes it clear that generally $\underline{\mathcal{O}_{S}}\left(t\right)\underline{\mathcal{Q}_{S}}\left(t\right)\neq\left(\underline{\mathcal{O}_{S}\mathcal{Q}_{S}}\right)\left(t\right)$. To correct for this discrepancy, we define

[TABLE]

in terms of new “noise” operators $\mathcal{E}_{\mathcal{O}_{S}}\left(t\right)$ that are essentially defined to enforce the consistency of operator products such as $\mathcal{O}_{S}\left(t\right)\mathcal{Q}_{S}\left(t\right)=\left(\mathcal{O}_{S}\mathcal{Q}_{S}\right)\left(t\right)$. Self-consistency forces noise operators to be mean-zero, as

[TABLE]

Furthermore, if the environment $E$ is Markovian, then noise operators are also uncorrelated with initial-time observables, i.e. $\braket{\mathcal{O}_{S}\mathcal{E}_{\mathcal{Q}_{S}}\left(t\right)}=0$, which means that noise operators can be neglected in the calculation of two-time correlators of the form $\braket{\mathcal{O}_{S}\mathcal{Q}_{S}\left(t\right)}$. To see why, we observe that a Markovian environment is essentially defined to satisfy

[TABLE]

with $\rho_{E}$ a time-independent steady state of the environment. If we enforce (174) for all states $\rho_{S}$, e.g. the maximally mixed state $\rho_{S}^{(1)}\propto\mathds{1}_{S}$ and $\rho_{S}^{(2)}\equiv\rho_{S}^{(1)}+\mathcal{O}_{S}$ with $\mathcal{O}_{S}$ any traceless operator on $S$ with operator norm $\mathinner{\!\left\lVert\mathcal{O}_{S}\right\rVert}\leq 1/\tr\mathds{1}_{S}$ (i.e. such that $\rho_{S}^{(2)}$ remains positive semi-definite, or a valid quantum state), then by linearity we find that

[TABLE]

which implies that the Markov approximation (174) holds even if we replace $\rho_{S}$ by any operator on $S$, and in particular

[TABLE]

We can therefore expand

[TABLE]

and invoke the Markov approximation in (176) to find that

[TABLE]

which implies

[TABLE]

Noise operators thus play no role in the calculation of correlators such as $C\left(t\right)$ in (15). In contrast, noise operators generally do play a role in the calculation of multi-time correlators of the form $\braket{\prod_{j}\mathcal{O}_{S}^{(j)}\left(t_{j}\right)}$Blocher and Mølmer (2019). Furthermore, these calculations generally require additional assumptions about the environment. To keep our discussion simple and general, we therefore exclude the effects of noise terms in Section IV.

Appendix O Short-time linear response and two-time correlators

Here we discuss the appearance of two-time correlation functions in the short-time linear response of correlators to perturbations of a Hamiltonian. Consider an initial Hamiltonian $H$ perturbed by an operator $V$ with $\mathinner{\!\left\lVert V\right\rVert}\ll\mathinner{\!\left\lVert H\right\rVert}$, where $\mathinner{\!\left\lVert\mathcal{O}\right\rVert}$ denotes the operator norm of $\mathcal{O}$, such that the net Hamiltonian is $\tilde{H}=H+V$. We denote the generator of Heisenberg time evolution under the perturbed (unperturbed) Hamiltonian by $\tilde{T}$ ($T$). These generators are related by

[TABLE]

where $\underline{V}$ is a superoperator whose action on operators $\mathcal{O}$ is defined by

[TABLE]

Through quadratic order in the time $t$ and linear order in the perturbation $\underline{V}$, we can say that

[TABLE]

Defining perturbed and unperturbed Heisenberg operators $\tilde{\mathcal{O}}\left(t\right)$ and $\mathcal{O}\left(t\right)$ that respectively satisfy $\braket{\tilde{\mathcal{O}}\left(t\right)}=\braket{e^{t\tilde{T}}\mathcal{O}}$ and $\braket{\mathcal{O}\left(t\right)}=\braket{e^{tT}\mathcal{O}}$, we thus find that for sufficiently small times $t$ and weak perturbations $V$,

[TABLE]

Two-time correlators $\braket{V\mathcal{O}\left(t\right)}$ and $\braket{\mathcal{O}\left(t\right)V}$, in addition to the expectation values $\braket{\left(V\mathcal{O}\right)\left(t\right)}$ and $\braket{\left(\mathcal{O}V\right)\left(t\right)}$, thus determine the short-time linear response of correlators $\braket{\mathcal{O}\left(t\right)}$ to perturbations $V$ of a Hamiltonian.