Relating out-of-time-order correlations to entanglement via multiple-quantum coherences

TL;DR

This paper demonstrates that multiple-quantum coherence spectra, a type of out-of-time-order correlation, can serve as an experimentally accessible entanglement witness, linking quantum chaos, entanglement, and fundamental physics.

Contribution

It establishes a formal connection between multiple-quantum coherences and entanglement, enabling experimental probing of multiparticle entanglement via OTOCs.

Findings

Multiple-quantum coherences act as entanglement witnesses.

OTOCs can directly probe multiparticle entanglement.

Results suggest a link between entanglement, thermodynamics, and quantum gravity.

Abstract

Out-of-time-order correlations (OTOCs) characterize the scrambling, or delocalization, of quantum information over all the degrees of freedom of a system and thus have been proposed as a proxy for chaos in quantum systems. Recent experimental progress in measuring OTOCs calls for a more thorough understanding of how these quantities characterize complex quantum systems, most importantly in terms of the buildup of entanglement. Although a connection between OTOCs and entanglement entropy has been derived, the latter only quantifies entanglement in pure systems and is hard to access experimentally. In this work, we formally demonstrate that the multiple-quantum coherence spectra, a specific family of OTOCs well known in NMR, can be used as an entanglement witness and as a direct probe of multiparticle entanglement. Our results open a path to experimentally testing the fascinating idea that entanglement is the underlying glue that links thermodynamics, statistical mechanics, and quantum gravity.

Figures

Equations

$$F_t(\varphi )\equiv \text{tr}[{\hat{\rho }}_0{\hat{\rho }}_f]=\text{tr}[{\hat{\rho }}_t{\hat{\rho }}_t(\varphi )]=\sum _m{I}_m({\hat{\rho }}_t)e^{-im\varphi },$$

$$F_t(\varphi )\equiv \text{tr}[{\hat{W}}_t^{†}(\varphi ){\hat{V}}_0^{†}{\hat{W}}_t(\varphi ){\hat{V}}_0{\rho }_0]=⟨W_t^{†}(\varphi ){\hat{V}}_0^{†}{\hat{W}}_t(\varphi ){\hat{V}}_0⟩$$

$$F_I({\hat{\rho }}_t,\hat{A})$$

$$={-2\frac{{\partial }^2{F}_t(\varphi )}{\partial {\varphi }^2}|}_{\varphi =0}\le F_Q({\hat{\rho }}_t,\hat{A})$$

$$I_m^{\max ,\mathrm{sep}}=\underset{N_{+}\in \{0,\dots ,N\}}{\max }\frac{(2N_{+})!}{4^{N_{+}}(N_{+}-m)!(N_{+}+m)!}.$$

$${\hat{H}}_{\text{int}}=-J/N{\hat{S}}_x^2-\Omega {\hat{S}}_z$$

$$\text{tr}[{\rho }_0{\rho }_f]$$

$$=\text{tr}[e^{-i{H}_{\text{int}}t}{\rho }_0{e}^{i{H}_{\text{int}}t}U(\varphi )e^{-i{H}_{\text{int}}t}{\rho }_0{e}^{i{H}_{\text{int}}t}{U}^{†}(\varphi )]$$

$$=\text{tr}[\rho U(\varphi )\rho U^{†}(\varphi )]=\text{tr}[\rho \rho (\varphi )]$$

$$=\text{tr}[\sum _{m^{\prime }}{\rho }_{m^{\prime }}\sum _m{\rho }_m{e}^{-im\varphi }]$$

$$=\sum _m\text{tr}[{\rho }_{-m}{\rho }_m]e^{-im\varphi }\equiv \sum _m{I}_m(\rho )e^{-im\varphi }$$

$$\dot{\rho }=-i[H_{\text{int}},\rho ]+\sum _n\mathcal{L}({\hat{\Gamma }}_n)\rho ,$$

$$\mathcal{L}({\hat{\Gamma }}_n)\rho =\sum _j{\hat{\Gamma }}_{n,j}\rho {\hat{\Gamma }}_{n,j}^{†}-\frac{1}{2}({\hat{\Gamma }}_{n,j}^{†}{\hat{\Gamma }}_{n,j}\rho +\rho {\hat{\Gamma }}_{n,j}^{†}{\hat{\Gamma }}_{n,j})$$

$$\text{tr}[{\rho }_0{\rho }_f]=\text{tr}[{\rho }_t{\rho }_t(\varphi )]$$

$$\text{tr}[{\rho }_0(\exp [{\mathcal{L}}_{\text{back}}t]{\rho }_t(\varphi ))]=\text{tr}[(\exp [\mathcal{L}t]{\rho }_0){\rho }_t(\varphi )]$$

$$\text{tr}[{\rho }_0\exp [{\mathcal{L}}_{\text{back}}\Delta t]{\rho }_t(\varphi )]$$

$$\iff \text{tr}[{\rho }_0({\mathcal{L}}_{\text{back}}{\rho }_t(\varphi ))]$$

$$\text{tr}[{\rho }_0({\mathcal{L}}_{\text{diss}}{\rho }_t(\varphi ))]=\text{tr}[({\mathcal{L}}_{\text{diss}}{\rho }_0){\rho }_t(\varphi )].$$

$$\text{tr}[\rho \rho (\varphi )]$$

$$=\text{tr}[\rho (1-iA\varphi -\frac{1}{2}{A}^2{\varphi }^2+\mathcal{O}({\varphi }^3))\rho (1+iA\varphi -\frac{1}{2}{A}^2{\varphi }^2+\mathcal{O}({\varphi }^3))]$$

$$=\text{tr}[{\rho }^2]-{\varphi }^2\text{tr}[{\rho }^2{A}^2-(\rho A{)}^2]+\mathcal{O}({\varphi }^3)$$

$$2\text{tr}[{\rho }^2{A}^2-(\rho A{)}^2]$$

$$={-\frac{d^2}{d{\varphi }^2}\sum _m{I}_m(\rho )e^{-im\varphi }|}_{\varphi =0}$$

$$=\sum _m{m}^2{I}_m(\rho ).$$

$$F_Q(\rho ,A)$$

$$={-2\frac{d^2}{d{\varphi }^2}{\left(\text{tr}\left[\sqrt{\sqrt{\rho }\rho (\varphi )\sqrt{\rho }}\right]\right)}^2|}_{\varphi =0}$$

$$f(\rho ,\rho (\varphi ))$$

$$=\text{tr}[\rho \rho (\varphi )]+1-\text{tr}[{\rho }^2].$$

$$I_m\left(\sum _k{p}_k{\rho }_k\right)\le \sum _k{p}_k{I}_m({\rho }_k).$$

$$I_m(p{\rho }^{(1)}+(1-p){\rho }^{(2)})=\text{tr}[(p{\rho }_{-m}^{(1)}+(1-p){\rho }_{-m}^{(2)})(p{\rho }_m^{(1)}+(1-p){\rho }_m^{(2)})]$$

$$\iff p^2\text{tr}[{\rho }_{-m}^{(1)}{\rho }_m^{(1)}]+(1-p{)}^2\text{tr}[{\rho }_{-m}^{(2)}{\rho }_m^{(2)}]+p(1-p)\text{tr}[{\rho }_{-m}^{(1)}{\rho }_m^{(2)}+{\rho }_{-m}^{(2)}{\rho }_m^{(1)}]$$

$$\iff p(1-p)\text{tr}[{\rho }_{-m}^{(1)}{\rho }_m^{(1)}+{\rho }_{-m}^{(2)}{\rho }_m^{(2)}-{\rho }_{-m}^{(1)}{\rho }_m^{(2)}-{\rho }_{-m}^{(2)}{\rho }_m^{(1)}]$$

$$\iff \text{tr}[({\rho }_{-m}^{(1)}-{\rho }_{-m}^{(2)})({\rho }_m^{(1)}-{\rho }_m^{(2)})]$$

$$\iff \text{tr}[(\Delta {\rho }_m^{†}\Delta {\rho }_m]$$

$$I_m(\rho )$$

$$=\text{tr}[(\sum _k{\rho }_{-m-k}^A\otimes {\rho }_{-m+k}^B)(\sum _{k^{\prime }}{\rho }_{m-k^{\prime }}^A\otimes {\rho }_{m+k^{\prime }}^B)]$$

$$=\sum _{k{k}^{\prime }}\text{tr}[({\rho }_{-m-k}^A{\rho }_{m-k^{\prime }}^A)\otimes ({\rho }_{-m+k}^B{\rho }_{m+k^{\prime }}^B)]$$

$$=\sum _k\text{tr}[({\rho }_{-m-k}^A{\rho }_{m+k}^A)\otimes ({\rho }_{-m+k}^B{\rho }_{m-k}^B)]$$

$$=\sum _k{I}_{m+k}({\rho }_A)I_{m-k}({\rho }_B),$$

$$|\text{CSS}(\theta ,\varphi )⟩={\left(\sin (\theta /2)|\uparrow ⟩+e^{i\phi }\cos (\theta /2)|\downarrow ⟩\right)}^{\otimes N}=\sum _{k=0}^N{\alpha }_{N,k}|N,k⟩,$$

Full text

Entanglement in quantum systems is a resource for quantum computation and communication and has been called the characteristic trait of quantum mechanics Schrödinger (1935). Recently, it has also been proposed D’Alessio et al. (2016) and experimentally tested in proof-of-principle experiments Kaufman et al. (2016); Neill et al. (2016) that quantum entanglement is in fact the key concept behind thermalization in isolated quantum systems. Essentially, the approach to equilibrium can be understood as the spreading of entanglement through the system’s degrees of freedom. In parallel, the concept of “scrambling” in many-body systems, which refers to the delocalization of quantum information over all of a system’s degrees of freedom, has gained great attention Hosur et al. (2016); Yao et al. (2016); Swingle et al. (2016); Rozenbaum et al. (2017); Zhu et al. (2016); Shen et al. (2017); Li et al. (2017); Bohrdt et al. (2017); Swingle and Chowdhury (2017), motivated by the finding that special models with thermal states “holographically dual” to black holes can thermalize and scramble quantum information at the fastest rate allowed by nature Maldacena et al. (2016); Sekino and Susskind (2008). The scrambling rate can be quantified through out-of-time-ordered correlators (OTOCs), which have been connected to entanglement via the Rényi entropy Hosur et al. (2016); Fan et al. (2017). However, the Rényi entropy is a strict entanglement monotone only for pure systems and hard to access experimentally, requiring resources that scale exponentially with the subsystem size as well as single-particle addressing. Therefore, it is desirable to establish experimentally accessible entanglement witnesses applicable to open as well as isolated quantum systems which can be used to quantify scrambling.

In this Letter, we formally show that a specific family of OTOCs, first developed in NMR under the name of the multiple-quantum coherence (MQC) spectra, are useful entanglement witnesses. The MQC protocol has been known for many years to be a suitable method to quantify the development of many-body quantum coherences Baum et al. (1985); Cappellaro (2014). Recently, it has been applied to describe the spreading of correlations Baum et al. (1985); Baum and Pines (1986); Munowitz et al. (1987); Sánchez et al. (2014) and as a signature of localization effects Álvarez and Suter (2010); Álvarez et al. (2013, 2015); Wei et al. (2016). While connections between MQCs and entanglement have been pointed out in Refs. Doronin (2003); Furman et al. (2008, 2009) and witnesses of two-particle entanglement have been constructed in Refs. Fel’dman and Pyrkov (2008); Fel’dman et al. (2012), to date a formal relation between the MQC spectrum and multiparticle entanglement generally applicable to mixed states does not exist. Here, we formally establish such a relation by deriving entanglement witnesses from the MQC intensities, as well as a relationship between MQCs and the quantum Fisher information (QFI) Tóth and Apellaniz (2014); Pezzé et al. (2016), a well-known witness of multiparticle entanglement.

To illustrate the power of these connections, we use the specific example of a long-range Ising model in a transverse field. We start the dynamics from a pure initial state, but show the applicability of the witness to mixed states by including decoherence arising from light scattering during the dynamics. This type of decoherence is relevant for a broad class of quantum systems. Our results demonstrate the existence of an experimentally accessible link between scrambling measured by OTOCs and entanglement, provided by the MQCs.

MQCs have a long tradition in NMR systems, which typically operate at high temperature. Measuring MQCs in pure and almost zero temperature initial states is now becoming feasible in cold-atom experiments, including Bose-Einstein condensates, ultracold atoms in cavities, or trapped ions Widera et al. (2008); Cucchietti (2010); Leroux et al. (2010); Douglas et al. (2015); Linnemann et al. (2016); Swingle et al. (2016); Macrì et al. (2016); Gärttner et al. (2017). Such experiments open the possibility to probe the rich information contained in an entangled state via MQCs.

We start by introducing the MQCs, which have been used as a means for quantifying quantum coherence Cappellaro (2014); Baum et al. (1985); Munowitz et al. (1987). Let $\left|\psi_{i}\right>$ be the eigenstates of a Hermitian operator $\hat{A}$ and $\lambda_{i}$ the corresponding discrete eigenvalues. We divide the density matrix of an arbitrary state $\hat{\rho}$ into blocks as $\hat{\rho}=\sum_{m}\sum_{\lambda_{i}-\lambda_{j}=m}\rho_{ij}\left|\psi_{i}\right>\left<\psi_{j}\right|=\sum_{m}\hat{\rho}_{m}$. Thus, $\hat{\rho}_{m}$ contains all coherences between states with eigenvalues of $\hat{A}$ that differ by $m$. An experimentally accessible quantifier of these MQCs is the Frobenius norm $I_{m}(\hat{\rho})=(\|\hat{\rho}_{m}\|_{2})^{2}=\text{tr}[\hat{\rho}_{m}^{\dagger}\hat{\rho}_{m}]$ called multiple-quantum intensity. The key idea is that $I_{m}$ can be directly accessed in an experiment that has the ability to reverse the dynamics that created the state of interest $\hat{\rho}$ from an initial fiducial state $\hat{\rho}_{0}$. In this context, the time reversal can be connected with the concept of many-body Loschmidt echoes, well-known probes of irreversibility and chaos Rhim et al. (1970, 1971); Pastawski et al. (2000); Zangara et al. (2016); Zangara and Pastawski (2017).

The protocol to measure $I_{m}$ is as follows Baum et al. (1985); Gärttner et al. (2017) (see Fig. 1): evolve $\hat{\rho}_{0}$ into $\hat{\rho}_{t}$ under a nontrivial unitary evolution $e^{-i\hat{H}_{\text{int}}t}$, apply $\hat{W}(\phi)=e^{-i\hat{A}\phi}$, evolve backward with $e^{i\hat{H}_{\text{int}}t}$ to $\hat{\rho}_{f}$, and finally measure the probability to find the system in the initial state $\text{tr}[\hat{\rho}_{0}\hat{\rho}_{f}]$ (if $\hat{\rho}_{0}$ is pure, this is the fidelity). Noting that $\hat{W}(\phi)\hat{\rho}_{m}\hat{W}^{\dagger}(\phi)=e^{im\phi}\hat{\rho}_{m}$ and using cyclic permutations under the trace, one finds

[TABLE]

where $\hat{\rho}_{t}(\phi)=\hat{W}(\phi)\hat{\rho}_{t}\hat{W}^{\dagger}(\phi)$. Thus, by Fourier transforming the signal with respect to $\phi$, one obtains the MQC spectrum $\left\{I_{m}(\hat{\rho}_{t})\right\}$ (see SM for details). For NMR systems typically operating at infinite temperature, this overlap measurement reduces to a magnetization measurement, making it possible to observe coherences as high as $m\sim 7000$ Álvarez et al. (2013, 2015). Nevertheless, the perturbative nature of the coherences present in highly mixed states, which facilitates experimental access of the MQCs, also implies that the underlying quantum complexity and entanglement content in those states are small in comparison to pure states. For pure states, measuring MQCs requires a fidelity measurement that encodes information about $N$-body correlations in an $N$-particle system 111While in NMR the overlap measurement consists in measuring the magnetization $\sum_{i}\hat{\sigma}_{i}^{z}$, for pure states the projector $\bigotimes_{i}(1-\hat{\sigma}_{i}^{z})$ needs to be measured, which involves spin correlations of any order.. Despite the fact that, in general, the resources required for measuring fidelity scale unfavorably with the system size, the feasibility of such a measurement has been demonstrated for up to 50 particles Gärttner et al. (2017), much beyond what is possible with schemes involving measuring entanglement entropy.

The connection between MQCs and OTOCs becomes apparent from $\hat{\rho}_{t}=e^{-i\hat{H}_{\text{int}}t}\hat{\rho}_{0}e^{i\hat{H}_{\text{int}}t}$. By defining $\hat{V}_{0}=\hat{\rho}_{0}$, if $\hat{V}_{0}\hat{\rho}_{0}=\hat{\rho}_{0}$ 222This can be achieved easily for pure states and for strongly mixed states encountered in NMR experiments., the above expression can be recast as Gärttner et al. (2017); Wei et al. (2016)

[TABLE]

where $\hat{W}_{t}(\phi)=e^{i\hat{H}_{\text{int}}t}\hat{W}(\phi)e^{-i\hat{H}_{\text{int}}t}$. $F_{t}(\phi)$ is therefore an OTOC function, a specific product of Heisenberg operators not acting in normal order. When $\hat{W}(\phi)$ and $\hat{V}_{0}$ are chosen to be initially commuting operators, then $F_{t}(\phi)=1-\langle|[\hat{W}_{t}(\phi),\hat{V}_{0}]|^{2}\rangle$. The growth of the norm of the commutator, i.e., the degree by which the initially commuting operators fail to commute at later times due to the many-body interactions generated by $\hat{H}_{\text{int}}$, is commonly used as an operational definition of the scrambling rate Hosur et al. (2016); Swingle et al. (2016); Yao et al. (2016). Scrambling can be interpreted as the process by which the information encoded in the initial state, through the interactions, is distributed over the other degrees of freedom of the system. This process makes it no longer possible to retrieve the initial information by local operations and measurements.

We are now in the position to state the main results of the Letter.

First, the second moment of the MQC spectrum [$F_{I}/2$, defined in Eq. (3)] provides a lower bound on the quantum Fisher information $F_{Q}$

[TABLE]

This expression becomes an equality for pure states $\hat{\rho}_{t}$.

The QFI has been introduced to quantify the maximal precision with which a parameter $\phi$ in the unitary $\hat{W}(\phi)=e^{-i\hat{A}\phi}$ can be estimated using the quantum state $\hat{\rho}$ as an input to an interferometer. It bounds the minimal variance of $\phi$ as $\Delta\phi\geq 1/\sqrt{F_{Q}(\hat{\rho},\hat{A})}$ (Cramér-Rao bound) (Braunstein and Caves, 1994). It has been shown that if $F_{Q}(\hat{\rho},\hat{A})>b_{k}\equiv nk^{2}+(N-nk)^{2}$ (note that $b_{k}\geq Nk$), where $n$ is the integer part of $N/k$, then $\hat{\rho}$ is $(k+1)$-particle entangled Pezzé and Smerzi (2009); Hyllus et al. (2012); Tóth (2012). To derive expression (3), we used the relation $F_{I}(\hat{\rho},\hat{A})=4\text{tr}[\hat{\rho}^{2}\hat{A}^{2}-(\hat{\rho}\hat{A})^{2}]$, which is a lower bound on the QFI Yadin and Vedral (2016) (see also SM ). The choice of the generator $\hat{A}$ can be optimized for detecting the entanglement of a given state using intuition from quantum metrology.

The relation (3) has a number of implications, the most direct one being that $F_{I}$ inherits the property of $F_{Q}$ of being a witness for multiparticle entanglement; i.e., $F_{I}>b_{k}$ implies $F_{Q}>b_{k}$ and thus $(k+1)$-particle entanglement. This allows us to establish an intimate connection between scrambling of quantum information and the buildup of entanglement. Namely, the $\phi$ dependence of the OTOC $F_{t}(\phi)$ encodes information about the entanglement content of the state $\hat{\rho}_{t}$. We also note that, for thermal states, QFI can be directly related to dynamic susceptibilities, as demonstrated in Refs. Hauke et al. (2016) [see explicitly Eq. (4)] and Pappalardi et al. (2017), which are well-known signatures of quantum critical behavior and phase transitions. Moreover, the QFI is a measure of macroscopic coherences, such as appear in ”cat states” Yadin and Vedral (2016).

Second, each individual $I_{m}$ by itself can be used as an entanglement witness. The quantity $F_{I}$ only characterizes the second moment of the MQC spectrum or, equivalently, only depends on the small-$\phi$ behavior of the measured observable $F_{t}(\phi)$, while the MQC spectrum, i.e., each individual $I_{m}$ contains much more detailed information about the state $\hat{\rho}$. To show that individual $I_{m}$ can witness entanglement, we use two properties SM : First, the $I_{m}$ are convex, or nonincreasing under mixing [$I_{m}[p\hat{\rho}_{1}+(1-p)\hat{\rho}_{2}]\leq pI_{m}(\hat{\rho}_{1})+(1-p)I_{m}(\hat{\rho}_{2})$ for $m\neq 0$]. Second, coherences of product states can be obtained from those of the constituent subensembles by $I_{m}(\hat{\rho}_{A}\otimes\hat{\rho}_{B})=\sum_{k}I_{m+k}(\hat{\rho}_{A})I_{m-k}(\hat{\rho}_{B})$. With these two properties, one can, for a given $m$, bound the maximal $I_{m}$ achievable on the set of separable states.

In the following, we outline how to derive such bounds for systems of spin $1/2$ particles. The detailed proof can be found in the Supplemental Material SM . The spins are described by Pauli operators $\hat{\sigma}_{j}^{\alpha}$, $\alpha=x,y,z$, $j=1,\dots,N$, with the eigenstates of $\hat{\sigma}_{j}^{z}$ denoted by $\left|\uparrow\right>_{j}$ and $\left|\downarrow\right>_{j}$. We calculate the maximal $I_{m}$ achievable with a separable state. Without loss of generality, we choose $\hat{A}=\hat{S}_{z}=\sum_{j}\hat{\sigma}_{j}^{z}/2$ 333For a general local operator $\hat{A}=\sum_{j}\mathbf{n}_{j}\cdot\mathbf{\hat{s}}_{j}$, we can just define the basis of the Hilbert space as a product of the eigenstates of the local spin operators $\mathbf{n}_{i}\cdot\mathbf{\hat{s}}_{j}$, which formally maps the problem to the case we consider explicitly.. It follows from the convexity that the maximum $I_{m}$ is assumed for pure states, which for separable states take the most general form $\bigotimes_{j}\left(\sqrt{p_{j}}\left|\uparrow\right>_{j}+e^{i\varphi_{j}}\sqrt{1-p_{j}}\left|\downarrow\right>_{j}\right)$. From the rule for building tensor products, it follows that $I_{m}$ is independent of $\varphi_{j}$ and is a quadratic polynomial in the $p_{j}$. Noting that $I_{m}$ is invariant under $p_{j}\rightarrow 1-p_{j}$, the maximum is assumed when all $p_{j}$ are either extremal (zero or one) or equal to $1/2$. For such a state with $N_{+}$ spins in the equal superposition state ($p=1/2$), $I_{m}$ can be calculated analytically and optimized numerically with respect to $N_{+}$, which yields

[TABLE]

Thus, if for a given state $\hat{\rho}$ and rotation generated by $\hat{A}$, one has $I_{m}>I_{m}^{\rm max,sep}$ for some $m$, then $\hat{\rho}$ must be entangled. Note also that $I_{N}$ is a witness of genuine $N$-partite entanglement Leibfried et al. (2005).

We now illustrate these results by applying them to the specific case of collective spin models. We consider a system of $N$ spin $1/2$ particles and the coherences with respect to the collective spin operator $\hat{A}=\hat{S}_{\mathbf{n}}=\sum_{j}\mathbf{\hat{s}}_{j}\cdot\mathbf{n}$, with $\mathbf{\hat{s}}_{j}=(\hat{\sigma}_{j}^{x},\hat{\sigma}_{j}^{y},\hat{\sigma}_{j}^{z})/2$ and a unit vector $\mathbf{n}=(n^{x},n^{y},n^{z})$. Thus, the spectrum of $\hat{A}$ consists of the (half) integers $M=-N/2,\ldots,N/2$, and we define the $m$th order coherence $\hat{\rho}_{m}$ as the block of the density matrix spanned by $\left|\phi_{M}\right>\left<\phi_{M+m}\right|$, where $\left|\phi_{M}\right>$ are the eigenstates of $\hat{A}$ with eigenvalue $M$. We study an all-to-all transverse-field Ising model

[TABLE]

where the spins are initially prepared in $\left|\psi_{0}\right>=\left|\uparrow\right>^{\otimes N}$. In the absence of decoherence, the dynamics is restricted to the symmetric Dicke manifold, which makes it very easy to numerically simulate the dynamics of large numbers of spins.

In Fig. 2, we illustrate the time evolution of the coherence spectrum $I_{m}$ for zero and nonzero transverse field. The QFI per particle, shown as a black line, is proportional to the variance of the coherence spectrum. The figure shows that the $I_{m}$ surpass the bounds for separable states in large parts of the spectrum. A complex pattern of self-interference emerges as soon as the coherences become distributed across the entire spectrum and the initially Gaussian state completely delocalizes in spin space. The two snapshots on the right show a relatively short evolution time, where $\hat{\rho}_{t}$ is a spin-squeezed near-Gaussian state, and a longer time, where the state becomes clearly non-Gaussian and the $I_{m}$ develop an intricate structure for both the pure Ising and the transverse-field Ising case. This snapshot corresponds to the longest time that has been measured experimentally for these parameters in Gärttner et al. (2017). At this time, the $I_{m}$ fall off at most linearly with $m$, while the bound decreases exponentially [cf. Eq. (4)]. This means that the degree ($I_{m}/I_{m}^{\rm max,sep}$) to which the entanglement bound is violated increases exponentially with $m$.

Next, we discuss the impact of decoherence for an example relevant to recent trapped-ion experiments Bohnet et al. (2016); Gärttner et al. (2017). We find that decoherence can substantially reduce the state overlap $F_{t}(\phi)$. However, for the parameters of Ref. Gärttner et al. (2017), detecting entanglement should be feasible. The main source of decoherence in these experiments is off-resonant light scattering, which can be captured by including Lindblad terms in the master equation SM . Specifically, we consider elastic Rayleigh scattering, which leads to coherence decay with rate $\Gamma_{el}$, and Raman scattering, i.e., incoherent transitions from $\left|\downarrow\right>$ to $\left|\uparrow\right>$ ($\Gamma_{du}$) and vice versa ($\Gamma_{ud}$). We emphasize that if $\Gamma_{du}=\Gamma_{ud}$, which is typically the case in the trapped-ion experiments, $\text{tr}[\hat{\rho}_{0}\hat{\rho}_{f}]=\text{tr}[\hat{\rho}_{t}\hat{W}(\phi)\hat{\rho}_{t}\hat{W}^{\dagger}(\phi)]$ in Eq. (1) still holds, and thus the $I_{m}$ can still be detected using the time reversal scheme SM .

The role of decoherence is illustrated in Fig. 3, where the choice of parameters is motivated by the experimental capabilities demonstrated in Ref. Gärttner et al. (2017). Typical experimental parameters are $J\lesssim 5\,$kHz, $t\lesssim 1\,$ms, and a total decoherence rate $\Gamma\approx 60\,$s*-1*, dominated by $\Gamma_{el}$. Numerical simulations were performed using an efficient density matrix symmetrization approach SM .

Comparing $F_{Q}(\hat{\rho}_{t},\hat{A})/N$ (black line) with the bound $F_{I}(\hat{\rho}_{t},\hat{A})/N$ (red dashed), one recovers $F_{Q}(\hat{\rho}_{t},\hat{A})=F_{I}(\hat{\rho}_{t},\hat{A})$ for pure states ($\Gamma=0$), but as decoherence rates are increased, the bound quickly becomes less tight. While the QFI decays slowly at small $\Gamma t$, the decay of the bound $F_{I}\sim e^{-N\Gamma t}$ is $N$-fold enhanced compared to the single-particle decay rate $\Gamma$ because the global state overlap $\text{tr}[\hat{\rho}_{0}\hat{\rho}_{f}]$ decays with this rate. The inverse spin-squeezing parameter Kitagawa and Ueda (1993), which also provides a lower bound on QFI, does not witness any entanglement for the case of Fig. 3(b), as the state is already strongly oversqueezed.

Figures 3(c) and 3(d) show the coherence spectra for two values of $\Gamma$. The main effect of dephasing is a global decay of the $I_{m}$ with $e^{-N\Gamma t}$, approximately independent of $m$, as expected at short times in an initially pure system. Nevertheless, even for strong dephasing, the $I_{m}$ still violate the entanglement bound for sufficiently large $m$, since the bound decreases exponentially with $m$, while the $I_{m}$ decay much more slowly. Therefore, even in the presence of single-particle decoherence processes, we observe that the $I_{m}$ remain useful entanglement witnesses in the considered scenario. Nevertheless, one needs to deal with the experimental challenge of detecting a small signal, especially for large $N$. We note, however, that in Ref. Gärttner et al. (2017), MQCs below $10^{-2}$ have been resolved.

In summary, we have derived inseparability criteria from the MQCs and a formal connection between MQCs and the QFI. Our results demonstrate that MQCs, a specific type of OTOCs, can serve as an experimentally accessible probe for detecting scrambling of quantum information and multiparticle entanglement in mixed states.

We thank Arghavan Safavi-Naini, Michael Wall, John Bollinger, Justin Bohnet, Graeme Smith, and Felix Leditzky for discussions. Supported by Defense Advanced Research Projects Agency (DARPA, W911NF-16-1-0576 through ARO), NSF Grant No. PHY 1521080, JILA-NSF Grant No. NSF-PFC-PHY-1734006, AFOSR-MURI, NIST, the DFG Collaborative Research Center SFB1225 (ISOQUANT), the Austrian Science Fund (FWF), through SFB FoQuS (No.  F4016-N23), the ERC Synergy Grant UQUAM, and the ERC Advanced Grant EntangleGen (Project-ID 694561).

I Extracting MQCs from the many-body echo sequence

In the time reversal sequence described in the main text, an initial state $\rho_{0}=\rho(t=0)$ evolves unitarily for a time $t$ into $\rho=e^{-iH_{\text{int}}t}\rho_{0}e^{iH_{\text{int}}t}$. Then, the rotation $U(\phi)=\exp(iA\phi)$ is applied transforming $\rho$ into $\rho(\phi)=U(\phi)\rho U^{\dagger}(\phi)$, and subsequently the time evolution is reversed to give the final state $\rho_{f}=e^{iH_{\text{int}}t}\rho(\phi)e^{-iH_{\text{int}}t}$ Measuring the projector on $\rho(0)$ amounts to

[TABLE]

Thus, the measured observable is the Fourier transform of the multiple-quantum coherence spectrum of the state $\rho$ and the MQCs $I_{m}$ can be extracted from it.

In the simulations presented in the main text, we focus on ensembles of $N$ (pseudo-) spin $1/2$ particles and consider coherences with respect to the generator of a global rotation $A=S_{\mathbf{n}}=n_{x}S_{x}+n_{y}S_{y}+n_{z}S_{z}$. In this case, the spectrum of $A$ consists of (half) integer values $-N/2\ldots N/2$, and thus $m=-N\ldots N$ and all coherences can be extracted by scanning the phase in $U(\phi)$ over the range $[0,2\pi]$. Generalizations to global rotations in higher spin systems are straight forward. Also, for $A=1/2\sum_{j}(n_{j}^{x}\sigma_{j}^{x}+n_{j}^{y}\sigma_{j}^{y}+n_{j}^{z}\sigma_{j}^{z})$, i.e. rotating each spin about an individual axis, none of our conclusions change.

II Effect of decoherence on the MQC detection scheme

The above equations assumed a unitary evolution. However, realistic experiments will suffer from diverse sources of decoherence. While these may suppress the MQCs, as discussed in the main text, we now show that the scheme for detecting the MQCs works equally well in a range of situations where decoherence is present.

Assuming decoherence due to a featureless Markovian bath, as is usually the case in trapped-ion or ultracold-atoms experiments, we can take its effect into account through a Lindblad master equation

[TABLE]

where

[TABLE]

is a Lindblad operator ($j$ is the particle index). Note, that here we use hats on the jump operators $\hat{\Gamma}$ to distinguish them from their corresponding jump rates $\Gamma$. The most relevant types of decoherence are spontaneous emission up ($\hat{\Gamma}_{du}=\sqrt{\Gamma_{du}}\left|\uparrow\right>\left<\downarrow\right|$) and down ($\hat{\Gamma}_{ud}=\sqrt{\Gamma_{ud}}\left|\downarrow\right>\left<\uparrow\right|$) as well as elastic dephasing ($\hat{\Gamma}_{el}=\sqrt{\Gamma_{el}}\left|\uparrow\right>\left<\uparrow\right|$, or, more generally, $\hat{\Gamma}_{el}^{\mathbf{n}}=\sqrt{\Gamma_{el}}\left|\uparrow\right>_{\mathbf{n}}\left<\uparrow\right|_{\mathbf{n}}$, where $\left|\uparrow\right>_{\mathbf{n}}$ is the up-eigenstate in the $\mathbf{n}$ direction on the Bloch sphere). Often, one faces situations where the dominant contribution is either only from elastic dephasing Uys et al. (2010) or where $\Gamma_{du}\approx\Gamma_{ud}$ Schindler et al. (2013). In both cases, the detection scheme remains valid, as we show now.

What we need to prove is that

[TABLE]

or

[TABLE]

where $\mathcal{L}\rho=-i[H,\rho]+\mathcal{L}_{\text{diss}}[\rho]$ and $\mathcal{L}_{\text{back}}\rho=-i[-H,\rho]+\mathcal{L}_{\text{diss}}[\rho]$. By writing the time evolution as a Trotter expansion, we only have to prove the above for an infinitesimal time step, which can then be applied iteratively. Expanding the exponential for an infinitesimal step $\Delta t$, we have

[TABLE]

Using $\text{tr}(\rho_{0}[H,\rho_{t}(\phi)])=-\text{tr}([H,\rho_{0}]\rho_{t}(\phi))$ due to cyclic permutation under the trace, what remains to be shown is

[TABLE]

For this, we note $\text{tr}[\rho_{1}(\mathcal{L}(\hat{\Gamma})\rho_{2})]=\text{tr}[(\mathcal{L}(\hat{\Gamma}^{\dagger})\rho_{1})\rho_{2}]$. Using this relation together with $\mathcal{L}_{\text{diss}}\rho=(\mathcal{L}[\hat{\Gamma}_{el}]+\mathcal{L}[\hat{\Gamma}_{ud}]+\mathcal{L}[\hat{\Gamma}_{ud}^{\dagger}])\rho$, valid for $\Gamma_{ud}=\Gamma_{du}$, and the fact that $\hat{\Gamma}_{el}=\hat{\Gamma}_{el}^{\dagger}$, demonstrates Eq. (12) and in consequence also Eq. (10).

III Derivation of the lower bound on quantum Fisher information

The statement $4\text{tr}[\rho^{2}A^{2}-(\rho A)^{2}]=2\sum_{m}I_{m}(\rho)m^{2}$ used in the main text (time arguments have been dropped) can be proven by expanding $\text{tr}[\rho\rho(\phi)]$ in $\phi$ around $\phi=0$

[TABLE]

Taking the second derivative with respect to $\phi$ this gives

[TABLE]

The inequality $F_{Q}(\rho,A)\geq 4\text{tr}[\rho^{2}A^{2}-(\rho A)^{2}]$ Girolami and Yadin (2017) follows from the relation of the QFI with the Uhlmann fidelity $f(\rho,\rho(\phi))$ Jozsa (1994); Liu et al. (2014)

[TABLE]

and the inequality Miszczak et al. (2009)

[TABLE]

IV Properties of the MQCs

Classification of MQCs in common categories of coherence measures: In the main text, we define the multiple quantum coherences of a quantum state $\rho$ with respect to a hermitian operator $A$ as $\rho_{m}=\sum_{M-M^{\prime}=m}\rho_{M^{\prime}M}\left|M^{\prime}\right>\left<M\right|$, where $\left|M\right>$ are the eigenstates of $A$ with eigenvalue $M$. The $\rho_{m}$ are identical with the $\delta$-coherences defined in Ref. Yadin and Vedral (2016), which quantify the macroscopicity of a state Kwon et al. (2018). The function $F_{I}(\rho,A)=\sum_{m}m^{2}I_{m}(\rho)$ defined in the main text as well as the quantum Fisher information are examples for measures of $\delta$-coherence. The $\delta$-coherence is a special case of the translationally covariant (TC) coherence defined in Ref. Marvian and Spekkens (2016). In that definition, a state $\rho=\sum_{m}\rho_{m}$ is called incoherent with respect to $A$ if it is invariant under the action of the unitary generated by $A$, i.e., if $e^{-iA\phi}\rho e^{iA\phi}=\sum_{m}\rho_{m}e^{-im\phi}=\rho$. This is the case if and only if $\rho_{m}=0\quad\forall\,m\neq 0$. The generalized rotation $U(\phi)$ generated by $A$ is called a translation, whence the TC coherence takes its name. The trace norm $\lvert\rho_{m}\rvert_{1}=\text{tr}[\sqrt{\rho_{m}^{\dagger}\rho_{m}}]$, Marvian and Spekkens (2016) as well as the quantum Fisher information Girolami (2014) are examples for measures of TC-coherence.

These definitions quantify the type of coherence that is relevant for tasks appearing, e.g., in quantum metrology where a distinguished operator $A$ exists. For example, in quantum phase estimation $A$ generates translations of a phase that one seeks to determine. For such tasks, off-diagonal entries in the density matrix that connect states with the same eigenvalue of $A$ do not constitute relevant coherences, as they are insensitive towards the translations generated by $A$ Marvian and Spekkens (2016). Except in the case where the spectrum of $A$ is non-degenerate, these definitions deviate from the one of Ref. Baumgratz et al. (2014), which considers all off-diagonal entries of the density matrix as coherent.

Convexity: The MQCs of a mixture of two (or more) states ($\rho=\sum_{k}p_{k}\rho_{k}$, where $\sum_{k}p_{k}=1$) cannot become larger than the weighted sum of the MQCs of the components of the mixture:

[TABLE]

Proof: We prove the inequality for the special case $\rho=p\rho^{(1)}+(1-p)\rho^{(2)}$, which can be immediately generalized to arbitrary $\rho=\sum_{k}p_{k}\rho_{k}$. For $p=0$ or $p=1$ the inequality follows trivially, so we need only prove the case $p\neq 0,1$.

[TABLE]

where the last statement is true since the left side is just the (non-negative) Frobenius norm of $\Delta\rho_{m}$.

Direct product: We want to express the MQCs of a state $\rho$ that is a direct product of two subsystems $\rho=\rho_{A}\otimes\rho_{B}$ by the MQCs of $\rho_{A}$ and $\rho_{B}$. Collecting the $m$th order coherences of $\rho=\rho_{A}\otimes\rho_{B}$ (writing $\rho_{X}=\sum_{m}\rho_{m}^{X}$) we have to sum over all the tensor products of sectors such that the sum of the coherenes is $m$: $\rho_{m}=\sum_{k}\rho_{m-k}^{A}\otimes\rho_{m+k}^{B}$. In the case of a state of $N_{A}+N_{B}$ spin $1/2$ particles and $U(\phi)$ a global rotation the sum runs over indices such that $m-k\in[-N_{A},N_{A}]$ and $m+k\in[-N_{B},N_{B}]$. Thus for the $I_{m}$ we obtain

[TABLE]

where we used the distributive law of inner and outer products, $\text{tr}[A\otimes B]=\text{tr}[A]\text{tr}[B]$, and the fact that the $\text{tr}[\rho_{m}\rho_{m^{\prime}}]\neq 0$ only if $m=-m^{\prime}$.

V Deriving entanglement bounds for the individual $I_{m}$

Next, we sketch the proof for the entanglement bounds on the individual coherences $I_{m}$ stated in the main text. We first calculate the MQCs of a coherent spin state (CSS). A general CSS can be expressed in terms of fully symmetric Dicke states.

[TABLE]

where

[TABLE]

and

[TABLE]

are the symmetric Dicke states.

Without loss of generality (see also Sec. I), we restrict ourselves to the coherences with respect to the basis $\{\left|\uparrow\right>,\left|\downarrow\right>\}$, which usually denotes the eigenstates of $\sigma_{z}$, but can be the eigenstates of any spin operator $\mathbf{n}\cdot\boldsymbol{\sigma}$. Thus, we obtain for the CSS

[TABLE]

and thus the coherence spectrum of the CSS with $\theta=\pi/2$ is

[TABLE]

We now proceed to derive entanglement bounds for the $I_{m}$. For this, we maximize $I_{m}$ on the set of fully separable states

[TABLE]

where $\sum_{k}\beta_{k}=1$ and $\rho_{k}^{(j)}$ general single particle densities. The convexity of the $I_{m}$ yields

[TABLE]

So the supremum of $I_{m}$ can be upper bounded by optimizing over pure product states

[TABLE]

Using their behavior under tensoring, the coherences $I_{m}$ can be calculated recursively from the coherences of the individual spin states $\rho_{j}=\left|\psi_{j}\right>\left<\psi_{j}\right|$:

[TABLE]

The single-particle MQCs are

[TABLE]

Thus $I_{m}$ is a second order polynomial in $p_{j}$. We also notice that $I_{0}(\rho_{j})$ and $I_{1}(\rho_{j})$ (and thus $I_{m}$) are invariant under the substitution $p_{j}\rightarrow 1-p_{j}$. Thus the maximum with respect to $p_{j}$ is assumed either at $p_{j}=1/2$ (if $d^{2}I_{m}/dp_{j}^{2}<0$) or at the edges $p_{j}=0,1$. Hence, in order to maximize $I_{m}$, each spin must either be in the state $(\left|\uparrow\right>+\left|\downarrow\right>)/\sqrt{2}$ or in one of $\left|\uparrow\right>$ and $\left|\downarrow\right>$. According to equation (28), $I_{m}$ does not change if a particle in $\left|\uparrow\right>$ or $\left|\downarrow\right>$ is added. Thus, the maximal $I_{m}$ for a separable state of $N$ spins equals the coherence of a coherent spin state of $N_{\text{opt}}$ particles with $\theta=\pi/2$, where $N_{\text{opt}}\in\{0,1\ldots N\}$ maximizes $I_{m}$.

VI Numerical method

VI.1 Symmetrized Liouville space

To study the MQCs numerically, we implemented the master equation (7) with $H$ the one-axis twisting Hamiltonian and with Lindblad operators $\hat{\Gamma}_{du}$, $\hat{\Gamma}_{ud}$, and $\hat{\Gamma}_{el}$. As long as the interactions are all-to-all and the decoherence processes affect all particles in the same way, Eq. (7) is invariant under exchange of particles. Since the global rotations involved in the sequence also preserve this symmetry and since the initial state is fully symmetric under particle exchange, the dynamics is restricted to the space of fully symmetric density matrices Sarkar and Satchell (1987); Hartmann (2016); Xu et al. (2013), allowing for an efficient numerical implementation.

In general, any density matrix can be represented as

[TABLE]

where

[TABLE]

and $\sigma_{\alpha_{j}}\in\{1,\sigma_{z},\sigma_{+},\sigma_{-}\}$. Here, $\alpha$ represents the vector of the single particle states for a given basis state. The symmetry constraint now means that this Liouville space basis can be restricted to symmetrized states

[TABLE]

where $(n_{z},n_{+},n_{-})$ are the number of occurrences of $(z,+,-)$ in $\alpha$ (the index $1$ occurs $n_{1}=N-n_{z}-n_{+}-n_{-}$ times) and $\chi(\alpha)$ is a permutations of the indices. There are $N!$ permutations but many of them generate identical states. We choose the normalization factor $\mathcal{N}=n_{z}!n_{+}!n_{-}!n_{1}!$ such that all (different) states occur with unit weight. Thus the possible symmetric states are effectively given by the Fock states of $N$ bosons on a 4-site lattice, which gives a dimension of

[TABLE]

which is a tremendous reduction of dimensionality compared to the $4^{N}$ basis states needed to represent a general $N$-particle density matrix. We will later make use of other basis choices, like $\sigma_{\alpha_{j}}\in\{1,\sigma_{x},\sigma_{y},\sigma_{z}\}$, but in the context of the MQC spectrum the $+/-$ basis is very intuitive since the block structure in terms of coherence blocks $m=n_{+}-n_{-}$ is already built in.

In order to carry out calculations on the symmetrized subspace we only have to write the master equation as well as the initial state and the observables in terms of the coefficients $c_{\alpha}$.

• •

Initial state:

[TABLE]

• •

Liouvillian: In order to calculate the matrix elements of the relevant Liouvillian operator we calculate its action on each basis state $(n_{z},n_{+},n_{-})$, i.e., we determine the coefficients $a_{\alpha\beta}$ in $\mathcal{L}[\rho_{\alpha}]=\sum_{\beta}a_{\alpha\beta}\rho_{\beta}$

Interactions:

[TABLE]

This term conserves $n_{+}$ and $n_{-}$ and is thus block diagonal if written in blocks of fixed $(n_{+},n_{-})$. These blocks are tri-diagonal.

Incoherent terms: We consider the Linbald terms, e.g.,

[TABLE]

which gives

[TABLE]

[TABLE]

[TABLE]

Here, we can see that $n_{+}$ and $n_{-}$ are conserved and if $\Gamma_{\text{ud}}=\Gamma_{\text{du}}$, the whole dissipative term is diagonal, i.e., just leads to an exponential decay of coherences. The joint action of one-axis twisting and dissipation is block diagonal with block size at most $N+1$, which greatly simplifies the problem.

With the above, the master equation can be written in terms of the coefficients $c_{\alpha}$:

[TABLE]

Thus the matrix-elements calculated above have to be transposed to yield the master equation for coefficient vector $c_{\alpha}$. We used that the basis states $\rho_{\alpha}$ are orthogonal, i.e. $\text{tr}[\rho_{\alpha}\rho_{\beta}^{\dagger}]\propto\delta_{\alpha\beta}$.

• •

Observables are expressed in terms of the $c_{\alpha}$ by expressing them in terms of the basis states $\rho_{\alpha}$.

We will repeatedly make use of

[TABLE]

The factor $2^{n_{1}+n_{z}}$ is due the fact that $\text{tr}[1^{2}]=\text{tr}[\sigma_{z}^{2}]=2$, and could be removed by defining the basis in terms of $1/2$ and $\sigma_{z}/2$. Similarly, one could normalize the basis states by the number of permutations, but then those factors would appear in the initial state.

Components of total spin (or populations $N_{\uparrow}=N/2+\langle S_{z}\rangle$):

[TABLE]

[TABLE]

[TABLE]

The all up probability:

[TABLE]

Muliple quantum coherences:

[TABLE]

[TABLE]

Second moment of the collective spin:

[TABLE]

$m$th moment:

[TABLE]

where $P_{n}$ is the probability to find $n$ spins in state $\left|\uparrow\right>$. $P_{n}$ is easier to directly compute by combinatoric arguments than $\langle S_{z}^{m}\rangle$. To calculate the moments of the spin along other axes $S_{\mathbf{n}}^{m}$, one can simply rotate the state into the corresponding direction and then calculate $S_{z}^{m}$.

Probability $P_{n}$ for $n$ spins in $\left|\uparrow\right>$:

[TABLE]

where binomial coefficients with $n<k$ are defined to be zero and ${{}_{2}}F_{1}$ denotes a hypergeometric function.

VI.2 Efficient implementation

We noticed that the the evolution under one-axis twisting and dissipation has a particularly simple form in the $(n_{z},n_{+},n_{-})$ basis. The Liouvillian $\mathcal{L}_{\text{int}}$ is block-diagonal with blocks of size $\leq N+1$ and the blocks are tri-diagonal matrices. The map $\exp[\mathcal{L}_{\text{int}}t]$ can thus be evaluated efficiently, scaling as $N^{2}$ for each block, giving an overall scaling of $N^{4}$ (there are $\mathcal{O}[N^{2}]$ blocks). Similarly, rotations are block-diagonal with blocks of size $\leq N+1$ in the $(n_{x},n_{y},n_{z})$ basis because, e.g., for a rotation about $y$, $n_{y}$ as well as $n_{x}+n_{z}$ is conserved. Moreover, the basis transformation between $(z+-)$ and $(xyz)$ basis also has such block structure since $n_{z}$ is conserved and $n_{x}+n_{y}=n_{+}+n_{-}$. We can thus write a rotation in the $(z+-)$ basis as a sequence of simpler operations: Transform to $(xyz)$, rotate, transform back. Each of these steps is a matrix-vector multiplication where each matrix has $\mathcal{O}[N^{4}]$ non-zero elements (recall that a rotation in the $(z+-)$ basis has $\mathcal{O}[N^{5}]$ non-zero elements).

We compute the matrix elements for rotations in the $(n_{x},n_{y},n_{z})$ basis by working out the action of $\exp[\mathcal{L}_{\text{rot}}\phi]$ on the basis states. Here, we give the matrix elements for $y$-rotations ($S_{y}=\sum_{i}\sigma_{y}^{(i)}/2$),

[TABLE]

with

[TABLE]

The transformations between the two different basis sets is given by

[TABLE]

with

[TABLE]

and for the reverse transformation

[TABLE]

with

[TABLE]

VI.3 More complex observables: QFI and entanglement entropies

Calculating more complex observables such as quantum Fisher information and Rényi- and von-Neumann entanglement entropies requires to determine the eigenvalues and eigenvectors of the density matrix. It is not obvious how this can be done given the coefficient vector in the symmetrized basis. In the following, we outline how to adopt the procedure described in Ref. Xu et al. (2013) to accomplish this. The main idea is that any permutation-symmetric state is block-diagonal if written in the Dicke basis of angular momentum states $\left|J,M,\beta\right>$, where $J=N/2,N/2-1\ldots 0$ (or $1/2$), $M=-J/2\ldots J/2$, and $\beta=1\ldots n_{N,J}$, where $n_{N,J}$ is the degeneracy factor of each $J,M$ pair. By block-diagonal we mean that $\left<J^{\prime}M^{\prime}\beta^{\prime}\right|\rho\left|JM\beta\right>=\delta_{J^{\prime}J}\delta_{\beta^{\prime}\beta}\left<JM^{\prime}\beta\right|\rho\left|JM\beta\right>$. Also, blocks with different $\beta$ but same $J$ are identical, and thus one only has to calculate one representative element of each block for each $J$ and keep track of the degeneracy factor $n_{N,J}$. This means that the number of matrix-elements that have to be determined is again $\sim N^{3}$, in fact the number of non-redundant and non-zero matrix elements is exactly the same as in any other basis, namely $\binom{N+3}{3}$.

In Ref. Xu et al. (2013), Xu et al. derive a method to recursively construct the matrix elements in the $\left|JM\right>$ basis from a symmetrized basis that is constructed with the basis operators $\left|\uparrow\right>\left<\uparrow\right|=(1+\sigma_{z})/2$, $\left|sd\right>\left<sd\right|=(1-\sigma_{z})/2$, $\left|\uparrow\right>\left<sd\right|=\sigma_{+}$, $\left|\downarrow\right>\left<\uparrow\right|=\sigma_{-}$. Thus, before we can apply the recursion, we have to transform the coefficient vector from the $(n_{z},n_{+},n_{-})$ basis to this $(n_{\uparrow},n_{\downarrow},n_{+},n_{-})$ basis. It is straight forward to show that this is accomplished by

[TABLE]

and for the reverse transformation:

[TABLE]

Note that the coefficient vector obtained in this way corresponds to symmetrized basis states that are not normalized in the same way as those in Ref. Xu et al. (2013). In order to obtain the same normalization, we have to multiply each coefficient $c_{(n_{\uparrow},n_{\downarrow},n_{+},n_{-})}$ by the number possible permutations $\frac{N!}{n_{\uparrow}!n_{\downarrow}!n_{+}!n_{-}!}$.

Now, finding the representation of $\rho$ in the basis $\left|JM\right>\left<JM^{\prime}\right|$, means that we want to express the coefficients $d_{JMM^{\prime}}=\left<JM\right|\rho\left|JM^{\prime}\right>$ in terms of the coefficients $c_{(n_{\uparrow},n_{\downarrow},n_{+},n_{-})}$. Before we outline the recursive method, we note that the basis states of both bases are eigenstates of the superoperators $\mathcal{S}_{L}[\cdot]=S_{z}\cdot$ and $\mathcal{S}_{R}[\cdot]=\cdot S_{z}$, where $S_{z}=\sum_{i}\sigma_{z}^{(i)}/2$

[TABLE]

Since two eigenstates with different eigenvalues are orthogonal to each other, this means that $c_{(n_{\uparrow},n_{\downarrow},n_{+},n_{-})}$ can only contribute to $d_{JMM^{\prime}}$ if the two constraints $M+M^{\prime}=n_{\uparrow}-n_{\downarrow}$ and $M-M^{\prime}=n_{+}-n_{-}$ are fulfilled. This reduces the number of summands that contribute to $d_{JMM^{\prime}}$ to $\leq N$. The second constraint reflects the fact that in both bases one can uniquely assign a coherence order $m=M-M^{\prime}=n_{+}-n_{-}$ to each basis state, i.e., in both representation there is a natural division of the density matrix into blocks of different coherence order.

We now outline the recursive procedure for determining the basis transformation matrix. First, we notice that $\left<N/2\,N/2\right|\rho\left|N/2\,N/2\right>=\left<\uparrow\ldots\uparrow\right|\rho\left|\uparrow\ldots\uparrow\right>=c_{(N000)}$, so we already know the first row of our basis transformation matrix. To get the next element $\left<N/2\,N/2\right|\rho\left|N/2\,N/2-1\right>$ we notice that using the properties of the angular momentum eigenstates $\left|JM\right>$

[TABLE]

where $S_{-}=\sum_{i}\sigma_{-}^{(i)}$. We have applied $S_{-}$ form the right to each basis state and then made use of our knowledge of the previously calculated row of the matrix, i.e., we match the shifted indices of the coefficients to the non-zero entries of the previous row. This way, we can recursively calculate all other elements $d_{N/2\,N/2\,M}$. Then, exploiting that the density matrix is Hermitian, we get $d_{N/2\,M\,N/2}=d_{N/2\,N/2\,M}$ from which we can calculate $d_{N/2\,M\,M^{\prime}}$ by applying the same recursion as before starting with $d_{N/2\,M\,N/2}$.

The next step is to calculate the matrix elements for other blocks with $J<N/2$. For this we take the trace on both sides of

[TABLE]

where the sum on left is constrained by $N=n_{\uparrow}+n_{\downarrow}+n_{+}+n_{-}$. The basis states with non-zero trace are the ones with $n_{+}=n_{-}=0$ on the left and $M=M^{\prime}$ on the right. Using in addition that coefficients of basis states with different $\mathcal{S}_{L}$ (or $\mathcal{S}_{R}$) eigenvalues do not depend on each other, we find

[TABLE]

With this relation we can calculate $d_{JJJ}$ using that we already know $d_{J^{\prime}JJ}$ with $J^{\prime}>J$. For example $d_{N/2-1\,N/2-1\,N/2-1}=(c_{(N-1,1,0,0)}-d_{N/2,N/2,N/2})/(N-1)=(c_{(N-1,1,0,0)}-c_{(N,0,0,0)})/(N-1)$. Form here, we can then again recursively calculate all $d_{N/2-1\,M,M^{\prime}}$ as described above. The degeneracy factor $n_{N,J}$ can be obtained from a Young tableau. A closed form expression is given by

[TABLE]

Given the density matrix of $\rho$ in the $\left|JM\right>\left<JM^{\prime}\right|$ basis, we can now diagonalize each block giving the eigenvalues $\lambda_{J,i}$ and eigenstates $\left|\psi_{J,i}\right>$ and calculate the QFI Braunstein and Caves (1994)

[TABLE]

where in the last expression the indices $k$ and $l$ run over the size $2J+1$ of the respective block.

In order to calculate entanglement entropies, we have to calculate partial traces. This is almost trivial in the basis $(n_{z},n_{+},n_{-})$. If we trace over $n$ out of $N$ particles, then a basis state will only contribute if all the $n$ particles are in state $1$. Taking the trace gives a factor $2^{n}$ for the normalization we use. Thus $c^{(N-n)}_{n_{z},n_{+},n_{-}}=2^{n}c^{(N)}_{n_{z},n_{+},n_{-}}$, where the particle number of the system is represented as an upper index. To calculate for example the von Neumann entropy of the reduced system, we then transform to the $\left|JM\right>\left<JM^{\prime}\right|$ basis, diagonalize and calculate $S_{E}=-\sum_{J}n_{N,J}\sum_{k}\lambda_{J,k}\log\lambda_{J,k}$. We proceed similarly for Rényi entropy and mutual information.