Dynamical quantum correlations of Ising models on an arbitrary lattice and their resilience to decoherence

It is now straightforward to calculate the dynamics of the Bloch vectors

$$\begin{equation} {\boldsymbol{S}}_j(t)=\textstyle\frac{1}{2}\{\langle\sigma_{j}^{x}(t)\rangle,\langle\sigma_{j}^{y}(t)\rangle,\langle\skew2\hat{\sigma}_j^z(t)\rangle\}. \end{equation} \tag{ 18 }$$

Because $\skew2\hat {\sigma }^{z}_j$ commutes with the Ising interaction the z component of spin is time independent, and given by $S_j^z=\frac {1}{2}g^{-}_j(0)=\frac {1}{2}\cos \theta _j$. The transverse spin components S^x_j(t) and S^y_j(t) can be obtained from the real and imaginary parts, respectively, of $\langle \skew2\hat {\sigma }^{+}_j(t)\rangle$. A straightforward application of equation (13) gives

$$\begin{eqnarray} \langle\skew2\hat{\sigma}^{+}_j(t)\rangle&=&\skew2\bar{f}_j(1)f^{\phantom *}_j(-1)\prod_{k\neq j}g^{+}_{k}(2J_{jk} t)\\ &=&\frac{1}{2}\mathrm{e}^{\mathrm{i}\phi_j}\sin\theta_j\prod_{k\neq j}\left(\cos 2J_{jk}t-\mathrm{i}\sin 2J_{jk}t\cos\theta_k\right). \nonumber \end{eqnarray} \tag{ 19 }$$

For the special case where all spins point along θ = π/2 at t = 0 there is pure decay of the Bloch vector without any rotation. In figure 3 we show the projection into the xy plane of the Bloch vector S₀(t) (where j = 0 labels the central site of a 55-site triangular lattice) for an initial state in which all spins point in a single direction lying outside of the xy plane (θ = π/4, ϕ = 0). The Bloch vector spirals inwards: the precession can be understood as a mean-field effect [33], with the spin rotating due to the average magnetization of the other spins, while the decay is due to the development of quantum correlations. Note that, in a finite system, the length S(t) = |S(t)| of the total Bloch vector ${\boldsymbol{S}}(t)=\sum _j{\boldsymbol{S}}_j(t)$ decays even at the mean-field level for any ζ ≠ 0. This decay, however, is due to the existence of a spatially inhomogeneous mean-field, and cannot be associated with the development of correlations.

Zoom In Zoom Out Reset image size

Spin–spin correlation functions can be calculated just as easily from equations (13) and (16). All two-point correlation functions can be calculated from the four quantities

$$\begin{equation} \mathcal{C}^{zz}_{jk}(t)\equiv\langle \skew2\hat{\sigma}_{j}^{z}(t)\skew2\hat{\sigma}_k^{z}(t)\rangle,\end{equation} \tag{ 20 }$$

$$\begin{equation} \mathcal{C}^{+z}_{jk}(t)\equiv\langle \skew2\hat{\sigma}_{j}^{+}(t)\skew2\hat{\sigma}_k^{z}(t)\rangle,\end{equation} \tag{ 21 }$$

$$\begin{equation} \mathcal{C}^{++}_{jk}(t)\equiv\langle \skew2\hat{\sigma}_{j}^{+}(t)\skew2\hat{\sigma}_k^{+}(t)\rangle,\end{equation} \tag{ 22 }$$

$$\begin{equation} \mathcal{C}^{+-}_{jk}(t)\equiv\langle \skew2\hat{\sigma}_{j}^{+}(t)\skew2\hat{\sigma}_k^{-}(t)\rangle \end{equation} \tag{ 23 }$$

and their complex conjugates. Since the Hamiltonian commutes with all $\skew2\hat {\sigma }_j^z$, the first one is given trivially by $\mathcal {C}^{zz}_{jk}=g^-_j(0)g^-_k(0)=\cos \theta _j\cos \theta _k$. The second one can be obtained from equations (16) and (19) as

$$\begin{eqnarray} \mathcal{C}_{jk}^{+z}&=&\skew2\bar{f}_j(1)f_j(-1)g^{-}_{k}(2J_{jk}t)\prod_{l\neq j,k}g^{+}_{l}(2J_{jl} t)\nonumber\\ &=&\frac{1}{2}\mathrm{e}^{\mathrm{i}\phi_j}\sin\theta_jg^{-}_{k}(2J_{jk}t)\prod_{l\neq j,k}g^{+}_{l}(2J_{jl} t), \end{eqnarray} \tag{ 24 }$$

and the third and fourth are

$$\begin{eqnarray*} \mathcal{C}_{jk}^{++}&=&\frac{1}{4}\mathrm{e}^{\mathrm{i}(\phi_j+\phi_k)}\sin\theta_j\sin\theta_k\prod_{l\neq k,j}g^{+}_{l}(2J_{jl}t+2J_{kl}t),\\ \mathcal{C}_{jk}^{+-}&=&\frac{1}{4}\mathrm{e}^{\mathrm{i}(\phi_j-\phi_k)}\sin\theta_j\sin\theta_k\prod_{l\neq k,j}g^{+}_{l}(2J_{jl}-2J_{kl}t). \end{eqnarray*}$$

These correlation functions can be used, for example, to calculate the time dependence of the spin squeezing parameter

$$\begin{equation} \xi(t)=\sqrt{\mathcal{N}}\frac{\Delta S_{\mathrm{min}}(t)}{S(t)}. \end{equation} \tag{ 25 }$$

Here ΔS_min(t) is defined to be the minimum uncertainty along a direction perpendicular to S(t). If we choose our x-axis to be in the direction of S(t), and define $\skew2\hat {S}_{\psi }=\frac {1}{2}\sum _{j} (\cos (\psi )\skew2\hat {\sigma }^z_j+\sin (\psi )\skew2\hat {\sigma }_j^y)$, then ΔS_min(t) is obtained by minimizing $(\langle \skew2\hat {S}_{\psi }^2\rangle -\langle \skew2\hat {S}_{\psi }\rangle ^2)^{1/2}$ over ψ. The squeezing parameter determines the phase sensitivity in a suitably performed Ramsey experiment, which is enhanced over the standard quantum limit whenever ξ < 1 [34]. For ζ = 0 (infinite-range interactions) the calculation was first performed in [34]. However, in many experimentally relevant situations the interactions have some finite range and the maximum achievable squeezing is diminished (figure 4).

Zoom In Zoom Out Reset image size

It is also possible to calculate correlation functions involving the application of spin operators at different times, which describe the propagation in time of a perturbation to the system. As an example, we can easily calculate dynamical response functions of the form

$$\begin{eqnarray} S^{ab}_{ij}(t_1,t_2)&=&\frac{1}{4}\left\langle\skew2\hat{\sigma}^a_i(t_1)\skew2\hat{\sigma}^b_j(t_2)\right\rangle\nonumber\\ &=&\frac{1}{4}\sin\theta_{\mathrm{i}}\sin\theta_j\mathrm{e}^{-2\mathrm{i}abJ_{ij}(t_2-t_1)}\prod_{k\neq i,j}g^{+}\left(2a t_1 J_{ik}+2b t_2 J_{jk}\right). \end{eqnarray} \tag{ 26 }$$

These can be combined to calculate dynamical response functions involving arbitrary Pauli matrices, some examples of which are shown in figure 5.

Zoom In Zoom Out Reset image size

As we will see in section 3, such dynamical correlation functions allow us to calculate the effect of spontaneous relaxation and excitation on the dynamics of the system. For instance, for all spins initially polarized along the x-axis, the effect of the sudden relaxation of a spin on site k at time t_s on the time dependence of $\skew2\hat {\sigma }^{x}_j$ (for j ≠ k) is given by

$$\begin{eqnarray} \left\langle\skew2\hat{\sigma}^{+}_k(t^{*}_{\mathrm{s}})\skew2\hat{\sigma}^{x}_j(t)\skew2\hat{\sigma}^{-}_k(t_{\mathrm{s}})\right\rangle&=&\mathrm{Re}\left[\mathrm{e}^{2\mathrm{i}J_{jk}(t-2t_{\mathrm{s}})}\prod_{l\neq j,k}\cos(2J_{jl}t)\right]\nonumber\\ &=&\cos(2J_{jk}t[1-2t_{\mathrm{s}}/t])\langle\skew2\hat{\sigma}^{x}_j(t)\rangle. \end{eqnarray} \tag{ 27 }$$

Note that this is the same time evolution we would obtain if no spontaneous relaxation had occurred, the k'th spin were simply absent, and the j'th spin were coupled to a longitudinal magnetic field of strength 2J_jk(2t_s/t − 1). The reason for this behavior is straightforward when considering the time evolution in the Schrödinger picture. The application of $\skew2\hat {\sigma }^{-}_k$ to the wave function at time t_s not only forces the k'th spin to point down between times t_s and t, but also destroys the piece of the initial wave function having weight into states with σ_k = −1. Hence it is as if the k'th spin pointed up for a time t_up = t_s, and down for a time t_down = t − t_s, thus contributing an inhomogeneous longitudinal magnetic field of strength 2J_jk(t_up − t_down)/t = 2J_jk(2t_s/t − 1) (the factor of 2 arises because the J_jk couple to Pauli matrices rather than the spin-1/2 matrices).

Given a density matrix ϱ describing a system coupled to a reservoir, it is always possible to express the expectation value of a system operator $\skew2\hat {\mathcal{A}}$ in terms of the system reduced density matrix ${\rho}={\rm tr}_{\mathrm {R}}\left [\varrho \right ]$ as $\langle \skew6\hat {\mathcal {A}}\rangle ={\mathrm { tr}}_{\mathrm {S}}[\rho \skew6\hat {\mathcal {A}}]$. Here tr_R(S) denotes a trace over the reservoir (system) degrees of freedom—we will drop these subscripts from now on, since all future instances of tr refer to a trace over system degrees of freedom only. In this language, the effect of a finite system-reservoir coupling is that an initially pure system density matrix ρ(0) ≡ |Ψ(0)〉〈Ψ(0)| will evolve into a mixed state (reflecting entanglement between the system and reservoir degrees of freedom). When the system–environment coupling is weak (i.e. small compared to the inverse of relevant system time scales) and the reservoir correlation time is small, the Born–Markov approximation is justified and the reduced system density matrix obeys a Markovian master equation of Lindblad form [55]. We choose a very general master equation appropriate for describing the various types of decoherence relevant to trapped ions [53], Rydberg atoms [52] and condensed matter systems such as quantum dots [5] and nitrogen vacancy centers [6]:

$$\begin{equation} \dot{\rho}=-\mathrm{i}\mathscr{H}(\rho)-\mathscr{L}_{\mathrm{ud}}(\rho)-\mathscr{L}_{\mathrm{du}}(\rho)-\mathscr{L}_{\mathrm{el}}(\rho), \end{equation} \tag{ 28 }$$

where

$$\begin{equation} \mathscr{H}(\rho)=\left[\mathcal{H},\rho\right],\end{equation} \tag{ 29 }$$

$$\begin{equation} \mathscr{L}_{\mathrm{ud}}(\rho)=\frac{\Gamma_{\mathrm{ud}}}{2}\sum_{j}\left(\skew2\hat{\sigma}^+_j\skew2\hat{\sigma}^-_j\rho+\rho \skew2\hat{\sigma}^+_j\skew2\hat{\sigma}^-_j-2\skew2\hat{\sigma}^{-}_j\rho\skew2\hat{\sigma}^+_j\right),\end{equation} \tag{ 30 }$$

$$\begin{equation} \mathscr{L}_{\mathrm{du}}(\rho)=\frac{\Gamma_{\mathrm{du}}}{2}\sum_{j}\left(\skew2\hat{\sigma}^-_j\skew2\hat{\sigma}^+_j\rho+\rho \skew2\hat{\sigma}^-_j\skew2\hat{\sigma}^+_j-2\skew2\hat{\sigma}^{+}_j\rho\skew2\hat{\sigma}^-_j\right),\end{equation} \tag{ 31 }$$

$$\begin{equation} \mathscr{L}_{\mathrm{el}}(\rho)=\frac{\Gamma_{\mathrm{el}}}{8}\sum_{j}\left(2\rho-2\skew2\hat{\sigma}^{z}_j\rho\skew2\hat{\sigma}^z_j\right). \end{equation} \tag{ 32 }$$

The first term involving a commutator describes coherent evolution due to the Ising interaction, and the various terms having subscripts 'ud', 'du' and 'el' correspond respectively to spontaneous relaxation, spontaneous excitation and dephasing⁹. Equation (28) has the formal solution $\rho (t)=\mathscr {U}(t)\rho (0)$, with

$$\begin{equation} \mathscr{U}(t)=\exp\left[-t\left(\mathrm{i}\mathscr{H}+\mathscr{L}_{\mathrm{ud}}+\mathscr{L}_{\mathrm{du}}+\mathscr{L}_{\mathrm{el}}\right)\right]. \end{equation} \tag{ 33 }$$

The exponential of super-operators is meant to be understood via its series expansion, in which the multiplication of two Lindblad super-operators implies composition $(\mathscr {L}_1\times \mathscr {L}_2)\left (\rho \right )=\mathscr {L}_1\left (\mathscr {L}_2\left (\rho \right )\right )$. Our goal in what follows is to compute the time dependence of an arbitrary operator $\skew6\hat {\mathcal {A}}$ at time t, given in the Schrödinger picture by

$$\begin{equation} \mathcal{A}(t)={\rm tr}[\skew2\hat{\mathcal{A}}\rho(t)]. \end{equation} \tag{ 34 }$$

An immediate simplification follows from the observation that

$$\begin{equation} \left[\mathscr{L}_{\mathrm{el}},\mathscr{H}\right]=\left[\mathscr{L}_{\mathrm{el}},\mathscr{L}_{\mathrm{du}}\right]=\left[\mathscr{L}_{\mathrm{el}},\mathscr{L}_{\mathrm{ud}}\right]=0. \end{equation} \tag{ 35 }$$

That the last two commutators vanish is less obvious than the first, but physically it has a very clear meaning: spontaneous relaxation/excitation on a site j causes the j'th spin to have a well defined value of σ^z_j, and thus to be unentangled with the rest of the system. Since the dephasing jump operator $\skew2\hat {\sigma }^z_j$ changes the relative phase between the states |σ^z_j = ± 1〉, whether spontaneous relaxation/excitation occurs before or after a dephasing event only affects the sign of the overall wave function, which is irrelevant. As a result, we can write

$$\begin{equation} \mathscr{U}(t)=\mathrm{e}^{-t\mathscr{L}_{\mathrm{el}}}\mathrm{e}^{-t\left(\mathrm{i}\mathscr{H}+\mathscr{L}_{\mathrm{ud}}+\mathscr{L}_{\mathrm{du}}\right)}, \end{equation} \tag{ 36 }$$

and the time dependence of an arbitrary observable $\mathcal {A}(t)=\mathrm {Tr} [\rho (t)\skew6\hat {\mathcal {A}}]$ can be written¹⁰

$$\begin{equation} \mathcal{A}(t)=\mathrm{Tr}\left[\mathrm{e}^{-t\left(\mathrm{i}\mathscr{H}+\mathscr{L}_{\mathrm{ud}}+\mathscr{L}_{\mathrm{du}}\right)}\rho(0)\mathrm{e}^{-t\mathscr{L}_{\mathrm{el}}}\skew2\hat{\mathcal{A}}\right]. \end{equation} \tag{ 37 }$$

Note that application of a superoperator does not commute with operator products ( $\mathscr {L}(\mathcal {O}_1)\mathcal {O}_2\neq \mathscr {L}(\mathcal {O}_1\mathcal {O}_2)$), and we use the convention that a superoperator should act on the operator immediately to its right. The effect of the time evolution due to $\mathscr {L}_{\mathrm {el}}$ can be understood by considering its effect on the Pauli operators:

$$\begin{equation} \mathrm{e}^{-t\mathscr{L}_{\mathrm{el}}}\skew2\hat{\sigma}^{x,y}_{j}=\mathrm{e}^{-\Gamma_{\mathrm{el}}t/2}\skew2\hat{\sigma}^{x,y}_j,\quad \mathrm{e}^{-t\mathscr{L}_{\mathrm{el}}}\skew2\hat{\sigma}^{z}_{j}=\skew2\hat{\sigma}^{z}_j. \end{equation} \tag{ 38 }$$

In light of equations (37) and (38), we are free to ignore the dephasing terms in the master equation at the expense of attaching a factor of e^−Γ_elt/2 to every operator σ^x,y_j occurring inside an expectation value:

$$\begin{equation} \mathrm{Tr}\left[\rho(t)\skew2\hat{\mathcal{A}}(\skew2\hat{\sigma}^x_j,\skew2\hat{\sigma}^y_j)\right]\rightarrow \mathrm{Tr}\left[\rho(t)\skew2\hat{\mathcal{A}}(\mathrm{e}^{-\Gamma_{\mathrm{el}}t/2}\skew2\hat{\sigma}^x_j,\mathrm{e}^{-\Gamma_{\mathrm{el}}t/2}\skew2\hat{\sigma}^y_j)\right], \end{equation} \tag{ 39 }$$

where ρ(t) on the right-hand side evolves under the master equation without the dephasing term.

Because the effects of dephasing are fully included by equation (39), the remaining problem is to compute the time dependence of operators whose expectation values are taken in a density matrix evolving simultaneously under $\mathscr {H}$, $\mathscr {L}_{\mathrm {ud}}$ and $\mathscr {L}_{\mathrm {du}}$:

$$\begin{equation} \mathcal{A}(t)=\mathrm{Tr}\left[\mathrm{e}^{-t\left(\mathrm{i}\mathscr{H}+\mathscr{L}_{\mathrm{ud}}+\mathscr{L}_{\mathrm{du}}\right)}\rho(0)\skew2\hat{\mathcal{A}}\right]. \end{equation} \tag{ 40 }$$

Formally the challenge of including the effects of spontaneous relaxation and excitation is related to the nontrivial commutation relation

$$\begin{equation} \left[\mathscr{H},\mathscr{L}_{\mathrm{ud(du)}}\right]\neq0. \end{equation} \tag{ 41 }$$

Physically, the obstacle is that spontaneous relaxation and excitation change the value of σ^z for the spin which they affect, and this change feeds back on the system through the Ising couplings. From the results on coherent dynamics presented in section 2, we know that time evolution under $\mathscr {H}$ alone is tractable. This suggests that we attempt to solve equation (40) by rewriting

$$\begin{equation} \mathrm{i}\mathscr{H}+\mathscr{L}_{\mathrm{ud}}+\mathscr{L}_{\mathrm{du}}=\mathrm{i}\mathscr{H}_{\mathrm{eff}}-\mathscr{R}, \end{equation} \tag{ 42 }$$

where $\mathscr {R}$ contains all and only terms that do not commute with $\mathscr {H}$, and then doing perturbation theory in $\mathscr {R}$. The above separation is accomplished by defining an effective (non-Hermitian) Hamiltonian and its corresponding superoperator

$$\begin{equation} \mathcal{H}_{\mathrm{eff}}=\mathcal{H}-\mathrm{i}\,\gamma\sum_j\skew2\hat{\sigma}_j^z-\mathrm{i}\,\mathcal{N}\frac{\Gamma_{\mathrm{r}}}{4}\quad \mathrm{and}\quad \mathscr{H}_{\mathrm{eff}}(\rho)=\mathcal{H}_{\mathrm{eff}}\rho-\rho\mathcal{H}_{\mathrm{eff}}^{\dagger}, \end{equation} \tag{ 43 }$$

and the recycling term

$$\begin{equation} \mathscr{R}(\rho)=\Gamma_{\mathrm{ud}}\sum_{j}\skew2\hat{\sigma}^{-}_j\rho\skew2\hat{\sigma}^{+}_{j}+\Gamma_{\mathrm{du}}\sum_{j}\skew2\hat{\sigma}^{+}_j\rho\skew2\hat{\sigma}^{-}_{j}, \end{equation} \tag{ 44 }$$

in terms of which the time evolution operator is $\mathscr {U}(t)=\mathrm {e}^{-t(\mathrm {i}\mathscr {H}_{\mathrm {eff}}-\mathscr {R})}$. In equation (43) we have defined $\gamma =\frac {1}{4}(\Gamma _{\mathrm {ud}}-\Gamma _{\mathrm {du}})$ and Γ_r = Γ_ud + Γ_du. Defining $\mathscr {U}_0(t)=\mathrm {e}^{-\mathrm {i}t\mathscr {H}_{\mathrm {eff}}}$, we can now expand the time evolution operator as a power series in $\mathscr {R}$ in order to obtain the time-dependent expectation value $\mathcal {A}(t)$:

$$\begin{eqnarray} &&\fl\mathcal{A}(t)= \sum_{n}\!\int_0^t \!\!\!\mathrm{d}t_{n}\cdots\!\!\int_0^{t_2}\!\!\! \!\mathrm{d}t_1{\rm tr}\left[\skew2\hat{\mathcal{A}}\mathscr{U}_0(t-t_{n})\mathscr{R}\mathscr{U}_0(t_{n}-t_{n-1})\cdots\mathscr{U}_0(t_2-t_1)\mathscr{R}\mathscr{U}_0(t_1)\rho(0)\right]\!. \end{eqnarray} \tag{ 45 }$$

This expansion is the underlying object being evaluated when Monte Carlo wave function methods [54] (quantum trajectories) are used to approximate the density matrix. In appendix A, we show in detail how the series in equation (45) leads to an expression for $\mathcal {A}(t)$ in terms of the closed-time path ordered correlation functions obtained in section 2, and the summation of that series is carried out in B. Here we will simply summarize the calculation, and explain in physical terms the essential structure of the Hamiltonian and decoherence that allows the result to be cast in closed form.

We begin by noting that when writing $\mathcal {A}(t)$ as a sum over closed-time path ordered correlation functions, each operator inserted along the forward leg of the time-contour is accompanied by its Hermitian conjugate appearing at the same time along the backward part of the time contour. If we had explicitly included an environment and attempted to trace over it (rather than starting with a Markovian master equation), this feature of the problem would emerge as a direct consequence of the Markov approximation. As a result, it suffices to describe any term in the series expansion by specifying the occurrence of operators on the forward time contour. To facilitate this description, we introduce notation describing the occurrence of operators belonging to a particular site j (see figure 6 for a summary). We take $\mathcal {R}^{\pm }_j$ to be the number of times the operator $\skew2\hat {\sigma }_j^{\pm }$ occurs along the forward time path, $\{t_{1}^j,\ldots ,t_{\mathcal{R}_j}^j\}$ to be the set of times at which jump operators are applied to site j, $\mathcal {R}_j=\mathcal {R}^{+}_j+\mathcal {R}^-_j$, κ_j = ± 1 depending on whether the operator at the latest time along the forward path is $\skew2\hat {\sigma }_j^{\pm }$, and

$$\begin{equation} \tau_j=(1-\alpha_j)\int_{0}^{t}\sigma_j^{z}(t). \end{equation} \tag{ 46 }$$

We will also use bold symbols ${\mathcalb{R}}$, κ and τ to specify the complete set of these variables on all lattice sites. Note that specifying $\mathcal {R}_j$ and κ_j determines both $\mathcal {R}^+_j$ and $\mathcal {R}_j^{-}$, since two consecutive (in the closed-time-path-ordering) applications of an operator $\skew2\hat {\sigma }_j^{\pm }$ gives zero. Therefore, we will only include the $\mathcal {R}_j$ and κ_j as explicit arguments in the correlation functions below. These variables are sufficient to determine the value of any term in the series expansion of $\mathcal {A}(t)$, so we do not need to keep track of the individual times at which each jump operator is applied.

Zoom In Zoom Out Reset image size

All nonzero terms in equation (45) are captured by summing over the $\mathcal {R}_j$ and κ_j, and integrating over the times $t^j_1,\ldots ,t^j_{\mathcal{R}_j}$, denoted

$$\begin{equation} \sum\!\!\int=\prod_{j=1}^{\mathcal{N}}\sum\!\!\int_{\!\!j}=\prod_{j=1}^{\mathcal{N}}\left(\delta_{\mathcal{R}_j,0}+\sum_{\kappa_j}\sum_{\mathcal{R}_j=1}^{\infty}\int_{0}^{t} \mathrm{d}t_{\mathcal{R}_j}^{j}\ldots\int_{0}^{t_2^j} \mathrm{d}t_1^{j}\right). \end{equation} \tag{ 47 }$$

If $\skew6\hat {\mathcal {A}}$ can be written as a product of operators $\skew2\hat {\sigma }^{b_j}_j$ (b_j = ±) on sites j contained in a set η, then $\mathcal {A}(t)$ can be compactly expressed (see appendix A) as

$$\begin{equation} \mathcal{A}(t)=\sum\!\!\int \mathcal{P}({\mathcalb{R}},{\boldsymbol{\kappa}},{\boldsymbol{\tau}})\mathcal{G}({\boldsymbol{s}},{\mathcalb{R}},{\boldsymbol{\tau}}). \end{equation} \tag{ 48 }$$

The prefactor

$$\begin{equation} \mathcal{P}({\mathcalb{R}},{\boldsymbol{\kappa}},{{\boldsymbol{\tau}}})=\prod_j\mathcal{P}_j(\mathcal{R}_j,\kappa_j,\tau_j)=\prod_j\left(\mathrm{e}^{-\Gamma_{\mathrm{r}}t/2}\left(\Gamma_{\mathrm{ud}}\right)^{\mathcal{R}^{-}_j}\left(\Gamma_{\mathrm{du}}\right)^{\mathcal{R}^{+}_j}\mathrm{e}^{-2\gamma\tau_j}\right), \end{equation} \tag{ 49 }$$

is closely related to the probability that a series of jumps described by the variables ${\mathcalb{R}}$, κ, and τ has occurred. Defining the symbol s as the vector of quantities s_j = φ_j/t − 2iγ, the correlation functions

$$\begin{equation} \mathcal{G}({\boldsymbol{s}},{\mathcalb{R}},{\boldsymbol{\tau}})=\prod_j\mathcal{G}_j(s_j,\mathcal{R}_j,\tau_j)=\prod_{j}\left(\mathrm{e}^{-\mathrm{i}\vartheta_j(\tau_j)}\left[p_j+g_j^{+}(s_j t)\right]\right), \end{equation} \tag{ 50 }$$

are of the general form presented in section 2. Careful bookkeeping reveals that the variables φ and ϑ defined in section 2 are given by

$$\begin{equation} \varphi_j=2t\sum_{k\in\eta}b_kJ_{jk},\end{equation} \tag{ 51 }$$

$$\begin{equation} \vartheta({\boldsymbol{\tau}})=\sum_{j}\vartheta_j(\tau_j), \end{equation} \tag{ 52 }$$

$$\begin{equation} \vartheta_j(\tau_j)= \left\{\begin{array}{@{}ll@{}} 0 & \mbox{if}~ j\in\eta,\\ -2\tau_j\sum_{k\in\eta}b_kJ_{jk} & \mbox{if}~ j\notin\eta.\end{array} \right. \end{equation} \tag{ 53 }$$

Note that the $\mathcal {R}_j$ dependence in $\mathcal {G}_j$ is hidden in the implicit dependence of p_j and g⁺_j on α_j (which for j∉η satisfies $\alpha _j=\delta _{\mathcal{R}_j,0}$). When evaluating the sums and integrals in equation (48), one must keep in mind that, whenever j∈η, terms with $\mathcal {R}_j\neq 0$ vanish because they contain the consecutive application of either σ⁺_j or σ⁻_j. Physically, the vanishing of such terms reflects the lack of coherence for any spin that has undergone even a single spontaneous spin flip.

The factorization of $\mathcal {P}$ into functions $\mathcal {P}_j$ of local site variables $\{\mathcal {R}_j,\kappa _j,\tau _j\}$ is a direct consequence of the single particle nature of the anti-Hermitian part of $\mathcal {H}_{\mathrm {eff}}$. Physically, this factorization occurs because the dissipation we are considering is uncorrelated from site to site (in contrast to the collective relaxation processes that arise, e.g. in the context of Dicke superradience [56]). The factorization of the correlation function $\mathcal {G}$ into functions $\mathcal {G}_j$ of local site variables $\{\mathcal {R}_j,\tau _j\}$ is a more surprising result, and depends crucially on the occurrence of jump operators at the same times on the forward and backward time evolution appendix A). The γ appearing in the argument of g⁺_j(s_jt = φ_j − 2iγt) affects the value of any term in equation (48) in which no jump operators are applied to site j (such that α_j = 1 and therefore g⁺_j(s_jt) ≠ 0). It arises from the term $-\mathrm {i}\,\gamma \skew2\hat {\sigma }^z$ in $\mathcal {H}_{\mathrm {eff}}$ (equation (43)), and decreases the expectation value of $\skew2\hat {\sigma }^z_j$ for γ > 0 (when spontaneous relaxation outweighs spontaneous excitation). This effect is often referred to as null measurement state reduction: gaining knowledge that the spin on site j has not spontaneously flipped affects the expected value when measuring $\skew2\hat {\sigma }^z_j$.

Because $\mathcal {G}$ and $\mathcal {P}$ both factorize, we need only to evaluate the quantities

$$\begin{equation} \Phi_{j}(s_j,t)=\sum\!\!\int_j\mathcal{P}_j(\mathcal{R}_j,\kappa_j,\tau_j)\mathcal{G}_j(s_j,\mathcal{R}_j,\tau_j) \end{equation} \tag{ 54 }$$

in terms of which

$$\begin{equation} \mathcal{A}(t)=\prod_j\Phi_j(s_j,t). \end{equation} \tag{ 55 }$$

The explicit dependence on s_j is included to remind the reader that, after the sums and integrals (over $\mathcal {R}_j$, κ_j, and $t_1^j\cdots t_{\mathcal{R}_l}^{j}$) have been carried out, s_j is the only site-dependent quantity on which Φ_j(s_j,t) depends. We evaluate these sums and integrals in B, obtaining

$$\begin{equation} \fl \Phi_j(s_j,t)\!=\!\left\{ \begin{array}{*2{l}} \!\!\! \mathrm{e}^{-\Gamma_{\mathrm{r}}t/2}p_j & \mbox{if}~j\in\eta,\vspace{0.2cm}\\ \!\!\! \mathrm{e}^{-\Gamma_{\mathrm{r}}/2}\left[\cos\!\left(\!t\sqrt{\!s_j^2-r}\,\right)+\left(\frac{\Gamma_{\mathrm{r}}}{2}+\mathrm{i}\varphi_j\cos\theta_j\right)t{\rm sinc}\!\left(\!t\sqrt{\!s_j^2-r}\,\right)\!\right] & \mbox{if}~j\notin\eta,\end{array} \right. \end{equation} \tag{ 56 }$$

where r = Γ_udΓ_du.

If $\skew6\hat {\mathcal {A}}$ also contains an operator $\skew2\hat {\sigma }_{l}^z(t)$ (with l∉η), we denote its expectation value by $\mathcal {A}^z_l(t)$. The insertion of $\skew2\hat {\sigma }_l^z$ must be dealt with at the point of equation (54). Keeping in mind the discussion surrounding equation (16), and remembering that s_l = φ_l/t − 2iγ, we must replace Φ_l(s_l,t) with

$$\begin{equation} \Psi_l(s_l,t)=\sum\!\!\int_l\mathcal{P}_l(\mathcal{R}_l,\kappa_l,\tau_l)\left((1-\alpha_l)\kappa_l+\alpha_l\frac{\mathrm{i}}{t}\frac{\partial}{\partial s_l}\right)\mathcal{G}_l(s_l,\mathcal{R}_l,\tau_l) \end{equation} \tag{ 57 }$$

Therefore we have

$$\begin{equation} \mathcal{A}_l^z(t)=\Psi_l(s_l,t)\prod_{j\neq l}\Phi_j(s_j,t), \end{equation} \tag{ 58 }$$

and in B we find

$$\begin{equation} \fl \Psi_{l}(s_l,t)=\mathrm{e}^{-\Gamma_{\mathrm{r}}/2}\left[\cos\!\left(\!t\sqrt{\!s_j^2-r}\,\right)+\left(\mathrm{i}s_l+2\gamma-\frac{\Gamma_{\mathrm{r}}}{2}\cos\theta_l\right)t\,{\rm sinc}\!\left(\!t\sqrt{\!s_j^2-r}\,\right)\!\right]. \end{equation} \tag{ 59 }$$

These equations reveal that correlation functions will generally undergo a qualitative transition in dynamics—from over-damped to oscillatory—whenever the condition s²_l = r is satisfied. This behavior is the most clearly manifest when the couplings J_ij have a simple structure, such as nearest neighbor or all-to-all. For instance, for nearest-neighbor coupling in 1D, assuming {Γ_el,Γ_ud,Γ_du} = {0,Γ,Γ}, and choosing the initial state to point along the x-axis, we find (ignoring boundary effects)

$$\begin{equation} S^x(t)=\frac{\mathcal{N}}{2}\mathrm{e}^{-3\Gamma t}\left[\cos\left(t\sqrt{4J^2-\Gamma^2}\right)+\Gamma t\,{\rm sinc}\left(t\sqrt{4J^2-\Gamma^2}\right)\right]^2, \end{equation} \tag{ 60 }$$

which becomes critically damped at Γ_c = 2J (see figure 7). It is interesting to note that this solution is similar in structure to the damping of a classical harmonic oscillator or a coherently driven two-level system (cf the weak-coupling limit of the spin boson problem [14, 57]). It is important to contrast this behavior with that of a single spin coupled to a Markovian bath, where the decoherence we consider only causes a damped-to-oscillatory transition in the presence of a transverse magnetic field; the Hamiltonian dynamics must be able to restore coherence in the basis for which the environment induces a measurement. In the present case, there is no transverse field, and it is not a priori obvious that such behavior should emerge. In fact, a simple mean-field estimate of the dynamics fails to capture the oscillatory-to-damped transition. Using a site-factorized ansatz for the density matrix, $\rho =\bigotimes _j\rho _{j}$, it is straightforward to see that [33]

$$\begin{equation} S^x_{\mathrm{MF}}(t)=\frac{\mathcal{N}}{2}\mathrm{e}^{-\Gamma t}\cos(4\,J\cos(\theta)t). \end{equation} \tag{ 61 }$$

Thus mean-field theory, which assumes an unentangled density matrix, always predicts under-damped dynamics; the transition to over-damped behavior captured by the exact solution depends crucially on the competition between decoherence and entanglement.

Zoom In Zoom Out Reset image size

In addition to providing an efficient way to compute arbitrary observables for an open many-body system, the above calculation provides significant insight into the interplay between decoherence and interactions. To keep the notation as simple as possible, we will focus on the case of nearest-neighbor coupling on a lattice with coordination number z, and calculate the time dependence of the transverse spin length of a single spin s^x(z,t) = 〈σ^x_j〉 (for an infinite system it is independent of j) for a state in which all spins are initially polarized along the x-axis. In the absence of decoherence but in the presence of a longitudinal magnetic field of strength h, application of results in section 2 gives

$$\begin{equation} s^x_h(z,t)=\frac{1}{2}\mathrm{Re}\left[\mathrm{e}^{-2\mathrm{i}ht}\cos^{z}(2Jt)\right]. \end{equation} \tag{ 62 }$$

In the absence of a longitudinal field but including equal rates of spontaneous relaxation/excitation (Γ_ud,du ≡ Γ, γ = 0), we find instead (from equations (55) and (56))

$$\begin{equation} s^x(z,t)=\mathrm{e}^{-(z+1)\Gamma_{\mathrm{r}}t/2}\left[\cos\left(t\sqrt{4J^2-\Gamma^2}\right)+\Gamma t\,{\rm sinc}\left(t\sqrt{4J^2-\Gamma^2}\right)\right]^{z}. \end{equation} \tag{ 63 }$$

However, it is instructive to temporarily hold off evaluating the sums and integrals implicit in the Φ_j, and work directly with equations (54) and (55):

$$\begin{equation} s^x(z,t)=\Phi_j(s_j,t)\prod_{k}^{\mathrm{NN}}\Phi_k(s_k,t)\prod_l^{\mathrm{D}}\Phi_l(s_l,t). \end{equation} \tag{ 64 }$$

In equation (64) we have divided the product over lattice sites into three parts: the j'th site, the nearest-neighbor sites (product labeled by NN) and the rest of the lattice (product labeled by D, reminding us that these sites are Disconnected from site j). First, we observe that $\Phi _j(s_j,t)=\frac {1}{2}\mathrm {e}^{-\Gamma _{\mathrm {r}}t/2}$. Next we evaluate $\prod ^{\mathrm {D}}_l\Phi (s_l=0,t)=1$ (the s_l = 0 because these sites are disconnected from the j'th site, and hence J_jl = 0, see equation (51)), which follows from a sum rule $\sum \!\int _{l}\mathcal {P}_l(\mathcal {R}_l,\kappa _l,\tau _l)\mathcal {G}_l(0,\mathcal {R}_l,\tau _l)=1$ (derivable from ${\mathrm { tr}}\left [\rho (t)\right ]=1$). It remains only to evaluate

$$\begin{equation} \sum^{\mathrm{NN}}\!\!\int\left(\prod^{NN}_{j}\mathcal{P}_j(\mathcal{R}_j,\kappa_j\tau_j)\mathcal{G}_j(s_j=2J,\mathcal{R}_j,\tau_j)\right), \end{equation} \tag{ 65 }$$

where we've adopted an abbreviated notation $\sum ^{\mathrm {NN}}\!\!\int =\prod _j^{\mathrm {NN}}\sum \!\int _j$ for the sums and integrals over the nearest-neighbor sites. Utilizing¹¹

$$\begin{equation} \prod_j^{\mathrm{NN}}\mathcal{G}_j(2J,\mathcal{R}_j,\tau_j)=2^{-\mathcal{R}}\mathrm{e}^{-2\mathrm{i}ht}\cos(2Jt)^{z-\mathcal{R}}=\frac{2}{2^{\mathcal{R}}}s^x_h(z-\mathcal{R},t), \end{equation} \tag{ 66 }$$

where $\mathcal {R}=\sum _{j}^{\mathrm {NN}}\mathcal {R}_{j}$, $\tau =\sum ^{\mathrm {NN}}_j\tau _j$, and h = Jτ/t, and noting that $s_{h}^x(z-\mathcal {R},t)$ only depends on τ and $\mathcal {R}$ (and not any other combinations of the local site variables), we can then write

$$\begin{equation} s^x(z,t)=\mathrm{e}^{-\Gamma_{\mathrm{r}}t/2}\sum_{\mathcal{R}=0}^{z}\int \mathrm{d}h\; P(\mathcal{R},h) s^x_h(z-\mathcal{R},t). \end{equation} \tag{ 67 }$$

Here

$$\begin{equation} P(\mathcal{R}',h')=\frac{t}{J}\sum^{\mathrm{NN}}\!\!\int\left(\prod^{\mathrm{NN}}_j\frac{\mathcal{P}_j(\mathcal{R}_j,\kappa_j,\tau_j)}{2^{\mathcal{R}_j}}\right)\delta(\tau-\tau')\delta_{\mathcal{R},\mathcal{R}'} \end{equation} \tag{ 68 }$$

is obtained by carrying out the sums and integrals while holding $\mathcal {R}$ and τ fixed (the factor of t/J arises when changing variables from τ to h in the remaining integral).

Equation (67) has a very suggestive form. The factor of e^−Γ_rt/2 out front is the single-particle contribution of T₁ decoherence, and would be present in the absence of interactions; it reflects the probability that the j'th spin has not spontaneously flipped before the time t (a prerequisite for having any coherence along x). As the z nearest neighbors evolve in time, they can undergo spin flips which cause them to fluctuate in time between pointing along +z and −z (figure 8(a)). Even once flipped, they influence (via the Ising couplings) the j'th spin in a manner formally equivalent to a longitudinal magnetic field of strength h = J(τ/t) (figure 8(b)). The quantity $P(\mathcal {R},h)$ describes the probability that $\mathcal {R}$ nearest neighbors have spontaneously flipped, collectively contributing an effective magnetic field of strength h to the time evolution of the j'th spin. The sum over $\mathcal {R}$ and integral over h then average the resulting dynamics for the j'th spin, $s_h^x(z-\mathcal {R},t)$, over the possible behaviors of its neighbors. For a given $\mathcal {R}$, the integral over h reduces the Bloch vector length by an amount depending on the width in h of the distribution $P(\mathcal {R},h)$. Physically, this integral captures the phase diffusion of the j'th spin due to the stochastic temporal fluctuation of its immediate environment (neighboring spins, as shown in figure 8). Thus we see very clearly that the interplay between spontaneous spin flips and coherent interactions leads to an emergent source of dephasing: flipping spins act as fluctuating magnetic fields (mediated by the Ising interactions) on other spins, even if these latter spins have not been directly affected by decoherence.

Zoom In Zoom Out Reset image size

Given the results in section 3, analytic calculation of the squeezing parameter in the presence of arbitrary decoherence rates Γ_el, Γ_ud and Γ_du is now straightforward. As can be seen in figure 9(a), the effect of T₂ decoherence is much less severe than that of T₁ decoherence. One reason for this behavior is that the minimum variance ΔS²_min occurs at an angle ψ (in the yz plane) only slightly deviating from the z-axis. Therefore dephasing, which can be thought of as random rotations of the individual spins around the z-axis, does not introduce much noise in the squeezed quadrature (see figure 9(b)). To the contrary, spontaneous relaxation/excitation processes introduce noise directly into the squeezed quadrature.

Zoom In Zoom Out Reset image size

In the absence of decoherence, the OAT Hamiltonian is known to give rise to $\mathcal {N}$-spin GHZ states¹²

$$\begin{equation} |\mathrm{GHZ}\rangle=\frac{|\Uparrow_x\rangle+\mathrm{i}^{\mathcal{N}+1}|\Downarrow_x\rangle}{\sqrt{2}}, \end{equation} \tag{ 70 }$$

at a time t* = ℏπ/4J, where |⇑_x〉 (|⇓_x〉) denotes the state where all spins point along the positive (negative) x-axis [58]. These entangled states afford Heisenberg-limited phase sensitivity [59], and are a resource for certain types of fault-tolerant quantum computation ([37] and references therein). However, they are also a canonical example of a fragile quantum state, and their usefulness is easily destroyed by decoherence [60]. The effect of dephasing on the production of GHZ state via OAT is well understood [60, 61]. With the results of section 3, however, we can easily calculate the effects of dephasing and spontaneous relaxation/excitation on the production of a GHZ state by OAT in a unified way. In this section we explicitly compare the effects of dephasing to the those of pure spontaneous relaxation ({Γ_el,Γ_ud,Γ_du} = {0,Γ,0}).

We first characterize the GHZ state by its phase coherence, obtained from the expectation value $\mathcal {C}={\mathrm { tr}}\left [\rho (t)\skew2\hat {\mathcal {C}}\right ]$ of the operator

$$\begin{eqnarray} \skew2\hat{\mathcal{C}}&=&|\Uparrow_x\rangle\langle \Downarrow_x\!|\nonumber\\ &=&\frac{1}{2^{\mathcal{N}}}\prod_{j}(\skew2\hat{\sigma}_j^z+\skew2\hat{\sigma}_j^{+}-\skew2\hat{\sigma}_j^{-}). \end{eqnarray} \tag{ 71 }$$

The quantity $\mathcal {C}$ characterizes the extent to which the superposition between the macroscopically distinct states |⇑_x〉 and |⇓_x〉 is quantum mechanical (rather than a classical mixture). Formally, $\mathcal {C}$ serves as a witness to $\mathcal {N}$-particle entanglement of the GHZ type, with entanglement guarantied whenever $|\mathcal {C}|>1/4$ is satisfied [37, 62]. Application of the results in section 3 yields

$$\begin{equation} \mathcal{C}(t)=\frac{1}{2^{\mathcal{N}}}\sum_{m=1}^{\mathcal{N}}\sum_{n=0}^{m}{{\mathcal{N}}\choose{m}}{{m}\choose{n}}\mathrm{e}^{-m\Gamma t}\Psi\left(2J[2n-m]-2\mathrm{i}\gamma,t\right)^{\mathcal{N}-m}, \end{equation} \tag{ 72 }$$

and in the absence of decoherence one finds $|\mathcal {C}(t^{*})|=1/2$. As can be seen in figure 10, the effect of spontaneous relaxation on the coherence $\mathcal {C}$ is comparable to (but worse than) the effect of dephasing. Both types of decoherence cause a loss of phase coherence when $\Gamma \sim J/\mathcal {N}$, with the factor of $\mathcal {N}$ responsible for the fragility of a GHZ state composed of a large number of spins.

Zoom In Zoom Out Reset image size

We can also directly calculate the metrological usefulness of a GHZ state prepared in the presence of spontaneous relaxation. The favorable sensitivity of the state ρ* ≡ ρ(t*) to rotations by angle Ω around the x-axis can be understood as the strong dependence of the expectation value of the parity operator $\skew2\hat {\pi }=\prod _{j}\skew2\hat {\sigma }^z_j$,

$$\begin{equation} P(\Omega)={\rm tr}\left[\rho^{*}(\Omega)\skew2\hat{\pi}\right]\end{equation} \tag{ 73 }$$

$$\begin{equation} ={\rm tr}\left[\rho^{*}\prod_{j}\left(\skew2\hat{\sigma}_j^z\cos\Omega-\skew2\hat{\sigma}_j^y\sin\Omega\right)\right] \end{equation} \tag{ 74 }$$

on the angle Ω (where ρ*(Ω) results from rotating ρ* about the x-axis by angle Ω)  [59]. The phase sensitivity of a GHZ state, denoted $\mathscr {M}$, is given (see [61]) by

$$\begin{equation} \mathscr{M}=\left.\frac{\left|\partial P(\Omega)/\partial\Omega\right|}{\Delta P(\Omega)}\right|_{\Omega=0}\approx\sum_j|{\rm tr}[\rho^{*}\skew2\hat{\sigma}^{y}_j\prod_{k\neq j}\skew2\hat{\sigma}^{z}_k]|, \end{equation} \tag{ 75 }$$

where ΔP(Ω) = 1 − P(Ω)² (taking into consideration that $\skew2\hat {\pi }^2=1$) is the uncertainty of the operator $\skew2\hat {\pi }$ calculated in the state ρ*. The approximation in equation (75) is simply that P(0) ≈ 0 and therefore ΔP(0) ≈ 1. This can be checked explicitly by looking at the large $\mathcal {N}$ limit of

$$\begin{equation} P(0)={\rm tr}\left[\rho^{*}\skew2\hat{\pi}\right]=\left(\frac{2\gamma\left(\mathrm{e}^{-\Gamma_{\mathrm{r}}t}-1\right)}{\Gamma_{\mathrm{r}}}\right)^{\mathcal{N}}. \end{equation} \tag{ 76 }$$

In the absence of decoherence, $\mathscr {M}=\mathcal {N}$ and the Heisenberg limit of phase sensitivity is obtained. By generalizing equation (58) to the case where $\skew2\hat {\sigma }_j^z$ is inserted on $\mathcal {N}-1$ of the sites, we find that in the presence of decoherence

$$\begin{equation} \mathscr{M}=\mathcal{N} \mathrm{e}^{-\Gamma t/2}\;\mathrm{Im}\left[\Psi(2J-2\mathrm{i}\gamma,t^{*})^{\mathcal{N}-1}\right]. \end{equation} \tag{ 77 }$$

This result is plotted in figure 11 for different types of decoherence. The enhancement in $\mathscr {M}$ survives T₁ decoherence only if $\mathcal {N}\Gamma _{\mathrm {r}} t^{*}\lesssim 1$ is satisfied (the scaling by $\mathcal {N}$ is shown in the inset). This result should be contrasted with the effect of T₂ decoherence ({Γ_el,Γ_ud,Γ_du} = {Γ,0,0}). In this case Ψ(2J,t*) = 1, yielding

$$\begin{equation} \mathscr{M}=\mathcal{N} \mathrm{e}^{-\Gamma t/2}, \end{equation} \tag{ 78 }$$

and hence the precision enhancement decays on a timescale that is independent of $\mathcal {N}$ (consistent with results in [61]). In contrast, the entanglement witness $\mathcal {C}$ decays at an $\mathcal {N}$-enhanced rate for either type of decoherence.

Zoom In Zoom Out Reset image size

In section 4.2, we showed that spontaneous relaxation significantly degrades the precision enhancement of a GHZ state unless $\Gamma \lesssim J/\mathcal {N}$ is satisfied. In this section, we will show that a time-resolved record of spontaneous relaxation events provides sufficient information to restore the phase enhancement under the much less stringent constraint Γ ≲ J. In particular, we are imagining a situation where spontaneous spin flips are accompanied by the real spontaneous emission of a photon, such that they can be measured by photodetection.

For reasons that will become clear in what follows, we take our initial state to have all spins pointing at an arbitrary angle θ (rather than θ = π/2, as assumed in section 4.2). Our goal is to evaluate the expectation value of the operator $\sum _{j}\skew2\hat {\sigma }^y_j\prod _{l\neq j}\skew2\hat {\sigma }_l^z$ at time t*, which was accomplished above by appealing directly to equations (59) and (56). Pursuing a strategy similar to that employed in section 3.5, we first observe that n spins initialized at θ = π/2, evolving in the absence of decoherence but in the presence of a longitudinal field of strength h, would yield the phase-sensitivity enhancement¹³

$$\begin{equation} \mathscr{M}(n,h)=n\left|\mathrm{Im}\left[\mathrm{e}^{-2\mathrm{i}ht^{*}}(\mathrm{i}\sin 2Jt^{*})^{n-1}\right]\right|\end{equation} \tag{ 79 }$$

$$\begin{equation} =\left|\sin(2ht^{*}+\pi(n-1)/2)\right|. \end{equation} \tag{ 80 }$$

Now we calculate $\mathscr {M}$ in the presence of decoherence, but we hold off evaluating the multiple sums and integrals that yield the functions Ψ_l(s_l,t) in closed form, and instead obtain $\mathscr {M}$ by working directly with equation (48)

$$\begin{equation} \mathscr{M}=\sum\!\!\int \mathcal{P}({\mathcalb{R}},{\boldsymbol{\kappa}},{\boldsymbol{\tau}})\mathcal{G}({\boldsymbol{s}},{\mathcalb{R}},{\boldsymbol{\tau}}). \end{equation} \tag{ 81 }$$

First, we evaluate $\mathcal {G}({\boldsymbol{s}},{\boldsymbol{R}},{\boldsymbol{\tau}})$, obtaining

$$\begin{equation} \mathcal{G}({\boldsymbol{s}},{\mathcalb{R}},{\boldsymbol{\tau}})\!=\!(\mathcal{N}-\mathcal{R})\cos^{2\mathcal{R}}\!\left(\frac{\theta}{2}\right)\!\mathrm{Im}\!\left[\mathrm{e}^{-2\mathrm{i} J\tau}g^{-}(s=2Jt-2\mathrm{i}\gamma t)^{\mathcal{N}-\mathcal{R}-1}\right]\!, \end{equation} \tag{ 82 }$$

which only depends on the $\mathcal {R}_j$ and τ_j only through their sums $\mathcal {R}=\sum _j\mathcal {R}_j$ and $h=(J/t)\sum _j\tau _j$. With the judicious choice $\theta =\pi -2\tan ^{-1}\left (\mathrm {e}^{2\gamma t^{*}}\right )$, we can rewrite

$$\begin{equation} g^{-}(2Jt-2\mathrm{i}\gamma t)=\frac{\mathrm{i}\sin(2Jt)}{\cosh2\gamma t}, \end{equation} \tag{ 83 }$$

and thus we obtain

$$\begin{equation} \mathscr{M}=\sum\!\!\int\mathcal{P}\left(\frac{\cos^{2\mathcal{R}}(\theta/2)}{\cosh^{N-1}(2\gamma t)}\right)\mathscr{M}(\mathcal{N}-\mathcal{R},J\tau/t). \end{equation} \tag{ 84 }$$

The initial value of θ, which places the initial spins slightly above the xy plane, was carefully chosen so that g⁻(2Jt − 2iγ)∝sin(2Jt). This finite tipping angle is required to precisely cancel the null measurement effect, which causes the z-projection of each spin to change in time due to the lack of emission of a photon, and thus ensures that at time t* the unflipped spins are brought down into the xy plane. Finally, we can write

$$\begin{equation} \mathscr{M}=\sum_{\mathcal{R}}\int \mathrm{d}h\; P(\mathcal{R},h)\mathscr{M}(\mathcal{N}-\mathcal{R},h=J\tau/t), \end{equation} \tag{ 85 }$$

where $P(\mathcal {R},h)$ is obtained by carrying out $\sum \!\int$ in equation (84) with $\mathcal {R}$ and τ = ht/J held fixed. As in section 3.5, we can interpret $P(\mathcal {R},h)$ as the probability to have $\mathcal {R}$ flipped spins contributing an effective magnetic field h to the dynamics of the remaining spins. For each particular value of $\mathcal {R}$ and h, the function $\mathscr {M}(n,h)$ (where $n=\mathcal {N}-\mathcal {R}$) yields the precision enhancement of a GHZ state produced from n spins in a longitudinal field of strength h.

If an experiment can record the times ( $t_{1},\ldots ,t_{\mathcal{R}}$) at which photons are emitted, thus gaining access to both $\mathcal {R}$ and

$$\begin{equation} h=\frac{J}{t}\sum_{j=1}^{\mathcal{R}}(2t_j-t), \end{equation} \tag{ 86 }$$

then that experiment produces the conditional density matrix ρ(n,h), corresponding to a GHZ state of n spins produced by OAT in the presence of the longitudinal field h. The effect of this field is simply to rotate the system around the z-axis by a (shot-to-shot random) angle δ = (ht* + π(n − 1)/2). However, rotation by a random angle only causes decoherence if the value of that angle is not known. Because the experiment measures δ indirectly via the photon emission record, it is possible to remove the effect of h in any experimental shot by applying the rotation operator $R(\delta )=\exp \left (-\mathrm {i}\, S^z\delta \right )$ to the conditional density matrix ρ(n,h). In this way we create

$$\begin{equation} \rho(n,0)=R(\delta)\rho(n,h)R^{\dagger}(\delta), \end{equation} \tag{ 87 }$$

which is GHZ state of the form in equation (70) containing n spins, and thus obtain a precision enhancement of n. Because the expected value of n will decay only at the bare rate Γ, there is no longer an enhancement by $\mathcal {N}$ in the decay of precision. This measurement-based coherent feedback can also be applied in the context of spin-squeezing, where once again it can vastly improve the metrological usefulness of a state generated by OAT in the presence of spontaneous relaxation.

Before concluding, we note that this feedback strategy could, in principle, be applied to situations where both spontaneous relaxation and spontaneous excitation are present, and even when the coupling constants J_ij are not uniform. However, in the former case it is necessary to independently record the photon emission record corresponding to excitation and relaxation processes, and in both cases it is necessary to obtain site-resolved (in addition to time-resolved) information about the photon emissions.