For its fundamental role in quantum physics, quantum entanglement has attracted lasting enthusiasm. Potential applications of entanglement in quantum information and quantum metrology lead to an urgent pursuit of massive entanglement [1–5]. To effectively characterize the massive entanglement in many-body systems, multiple quantum coherences (MQCs), which were originally introduced in nuclear magnetic resonance (NMR) systems [6–11], are employed to quantify the many-body quantum dynamics. With rapid progress in experimental technology, it is now possible to achieve many-body time-reversal Hamiltonian in pioneering trapped-ion experiments where the MQCs are measured [12]. The connection between MQCs and out-of-time-order correlations (OTOCs) is remarkably established. Through the appropriate choice of operators, a special type of fidelity OTOCs become exact MQCs [13]. Therefore, the MQCs are a powerful tool to diagnose massive quantum entanglement and quantum phase transition (QPT) [14, 15].

Spinor Bose–Einstein condensate (BEC) can exhibit rich quantum phases under an external field, and some of them are deeply entangled [16–18]. Recently, great breakthroughs have been made in preparing massively entangled twin-Fock (TF) state and Dicke state in a spinor BEC through quantum phase transitions [19, 20]. The key idea of such experiments is to adiabatically sweep the driving field across different phases and keep the system in its ground state, so that the system evolves slowly to its entangled quantum state [21, 22]. However, it is incredibly challenging to achieve the adiabatic condition in the limited lifetime of the BEC, particularly for up to tens of thousands of atoms in a real experiment. Although Luo et al and Zou et al estimate that thousands of atoms in the system are entangled through entanglement depth measurement [23, 24], the fidelity of the final quantum state to a genuine TF or Dicke state remains unclear. Compared to state tomography and randomized measurements [25–28], the entanglement depth itself is inadequate to fully characterize the final quantum state. In addition, the many-body dynamics of the condensate through the QPTs are incomplete because of the absence of real-time measurement of the entanglement depth. For such a large scale and highly entangled quantum system, it is worth developing new tools to extract more detailed information.

In this work, we study the MQC dynamics in a ferromagnetically interacting spin-1 BEC, which exhibits three quantum phases as the quadratic Zeeman shift is adjusted. We numerically simulate the QPTs in the spin-1 BEC by scanning the quadratic Zeeman shift. When the evolution of the condensate is adiabatic, the MQC spectrum shows dramatic changes around the quantum critical points, catching the QPTs perfectly. Limited to a finite-time evolution, however, the BEC has excitations, leading to oscillatory behavior in low-excitation states. Surprisingly, the MQC spectrum exhibits similar oscillatory behavior, indicating that MQC may act as a useful indicator to many-body dynamics. In fact, a few low-order MQCs, which are favored in experiments, are enough to clearly characterize the complex QPTs in the spinor BEC. To implement experimentally, the MQC spectrum can be obtained approximately with pesuo-echo, as pointed out in [15].

We start by briefly reviewing the MQCs, which were originally introduced in NMR systems and recently adopted to describe the many-body quantum systems [14]. We assume that the given Hermitian operator $\hat{A}$ and its eigenvectors ∣i⟩, satisfy $\hat{A}|i⟩=a_i|i⟩$. A density matrix can be expanded aswhere the block ${\hat{\rho }}_k$ is defined aswith ${\rho }_{i,j}=⟨i|\hat{\rho }|j⟩$.

The property of ${\hat{\rho }}_k$ becomes clear if a rotation ${\hat{R}}_{\varphi }={\mathrm{e}}^{-\mathrm{i}\hat{A}\varphi }$ acts on it, ${\hat{R}}_{\varphi }{\hat{\rho }}_k{\hat{R}}_{\varphi }^{†}={\mathrm{e}}^{-\mathrm{i}k\varphi }{\hat{\rho }}_k$. Obviously, ${\hat{\rho }}_k$ contains all coherences between eigenstates of $\hat{A}$ differing by k. With ${\hat{\rho }}_k$, the spectrum of MQCs is expressed in a standard form aswhere we have utilized the fact that ${\hat{\rho }}_k={\hat{\rho }}_{-k}^{†}$. For a pure quantum state, ${\sum }_k{\mathcal{J}}_k=1$. As is well known in NMR, non-zero ${\mathcal{J}}_{k\ne 0}$ are associated to quantum coherences.

We consider applications of MQCs in a spin-1 BEC and employ the MQCs to probe the QPTs. Under the control of an external magnetic or microwave field, a spin-1 BEC can be in different quantum phases, some of which are noticeable for deep entanglement. Under the single-mode approximation, the system is governed by the following Hamiltonian (ℏ = 1) [29, 30]where N is the total number of atoms in the BEC and ${\hat{N}}_{0,\pm 1}={\hat{a}}_{0,\pm 1}^{†}{\hat{a}}_{0,\pm 1}$ is the number operator of corresponding spin component. The operator ${\hat{a}}_{0,\pm 1}$ (${\hat{a}}_{0,\pm 1}^{†}$) annihilates (creates) an m_F = 0, ± 1 atom in the mode. In the Hamiltonian (3) we have used the fact that the z-component of the total magnetic quantum number is conversed and the system is initially prepared in the subspace N₊ − N₋ = 0. Obviously, the system keeps N₊ − N₋ = 0 during its evolution so that the linear Zeeman term by external field makes no contribution. The term proportional to c₂ is from the spin-dependent collision between atoms. Without loss of generality, we choose c₂ = − 2π × 2.7 Hz in the following numerical calculations for the case of ⁸⁷Rb [19, 30]. The ratio of the second-order Zeeman $-q{\hat{N}}_0$ and the spin-dependent interaction determines the ground state phase diagram as shown in figure 1(a). By diagonalizing the Hamiltonian (3), we obtain the energy gap Δ (the energy difference between the ground state energy and the first excitation state). The energy gap develops to a minimum at two critical points, leading to three quantum phases, polar (P) phase, broken-axisymmetry (BA) phase, and TF phase. When q/∣c₂∣ ≫ 0, the second-order Zeeman term dominates and the ground state is in the polar phase, where atoms tend to occupy the m_F=0 component. Contrarily, q/∣c₂∣ ≪ 0 leads to the TF phase, where atoms occupy the m_F = ± 1 component with equal weight. Between the two, it is the BA phase, where all three components are populated. The Dicke state in the BA phase and the TF state in the TF phase are remarkable for their deep entanglement [19, 20].

We now turn to the MQCs by choosing $\hat{A}={\hat{N}}_0$. Since the Hamiltonian equation (3) annihilates or creates pairs of m_F = 0 atoms, only even order MQCs remain. The MQC spectrum in ground states is shown in figure 1(b). At both ends of the phase diagram q/∣c₂∣ ≫ 1 (q/∣c₂∣ ≪ 1), the ground state of the spinor BEC is ∣N₀ = N⟩ (∣N₀ = 0⟩), leading to ${\mathcal{J}}_0=1$ and ${\mathcal{J}}_k=0$ if k ≠ 0. In contrast, it shows a wide spectrum in the BA phase with all ${\mathcal{J}}_k$ being small. The dramatic change of spectrum only occurs at critical points, clearly indicating the QPTs. Interestingly, as shown in figure 1(c), low-order MQCs are enough to catch the QPTs since their changes around the critical points are more prominent. This observation suggests that only low-order MQCs are actually needed, which is of great benefit in experiments.

In order to show the power of the MQCs further, we consider experimental generation dynamics of the TF state [19]. Both the theory and experiment suggests that an adiabatic ramp from q ≫ 2∣c₂∣ to q ≪ − 2∣c₂∣, is robust enough to generate a deeply entangled TF state [19, 20, 22]. In the simulation, we adopt the optimized ramp curve as given in [20] (see Appendix for details). The optimized ramp shows a slower scanning speed near the critical points. Due to the limited total evolution time, a real experiment can not follow the adiabatic process exactly and excitations are inevitable. As a result, the generated quantum state may significantly deviate from the ideal TF state. To characterize the many-body dynamical process and the final fidelity of the TF state, conversion efficiency p_c = 1 − ⟨N₀⟩/N is employed in experiments. Compared to the fidelity, F_K(t) = ∣⟨ψ(t)∣N₀ = K⟩∣², the conversion efficiency is insensitive to the QPTs and poorly catches the critical points, as shown in figure 2.

We resort to the MQCs to characterize the generation dynamics as an alternate protocol. The system is initially set to the polar state, i.e., all atoms occupy the m_F = 0 component. By performing real-time evolution under the Hamiltonian (3), we can obtain the time-dependent wave function and the corresponding MQC dynamics. Figure 2 presents the MQC dynamics with N = 1024 and total evolution time T = 6 s. As shown in figure 2(a), the MQC spectrum is similar to the ideal one, i.e., the MQCs suddenly change near critical points. However, oscillations of the spectrum are clearly visible in time-dependent evolution, suggesting that the spinor BEC does not always stay in the ground state.

Detailed dynamics of ${\mathcal{J}}_{0,2}$ are shown in figures 2(b) and 2(c). We have set the initial state to the polar state so that ${\mathcal{J}}_0(t=0)=1$ and ${\mathcal{J}}_2(t=0)=0$. As the system evolves before crossing the first critical point, the limited time scanning leads to clear deviation from the adiabatic evolution. ${\mathcal{J}}_0$ decreases with rather large amplitude oscillations. As a result of the constraint ${\sum }_k{\mathcal{J}}_k=1$, high order MQCs also oscillate. Once the system crosses the first critical point, all ${\mathcal{J}}_k$ components become small with even smaller rapid oscillations. After crossing the second critical point, ${\mathcal{J}}_0$ increases with large oscillations. Similarly, ${\mathcal{J}}_2$ also increases with oscillations. Obviously, one may approximately determine the P phase, the BA phase and the TF phase by observing the change of ${\mathcal{J}}_0$ and ${\mathcal{J}}_2$.

To further understand the oscillatory behavior of ${\mathcal{J}}_{0,2}$, we present the fidelities F_N and F_N−2 in figure 2(b) and F₀ and F₂ in 2(c). Interestingly, ${\mathcal{J}}_0$ approximately follows the oscillations of F_N in the P phase and F₀ in the TF phase. Such a similarity between ${\mathcal{J}}_0$ and F_N (or F₀) manifests the fact that the density matrix of the ground state in the P (or TF) phase has dominant diagonal elements. In addition, we also find similarities in the oscillation period and the phase between ${\mathcal{J}}_2$ and F_N−2 (or F₂), indicating that excitations are limited in the low-lying excited states.

As we have shown that one may catch the TF state with ${\mathcal{J}}_0$ alone, due to the similarity of ${\mathcal{J}}_0$ to F₀. Next we simulate the generation of the TF state in different total evolution time T with different atom number N, to systematically explore the relevance between ${\mathcal{J}}_0$ and F₀. As the total evolution time T increases, it is natural to expect that the generation process of the TF state becomes more adiabatic, and a higher fidelity to the TF state is reached in principle. However, the fidelity does not grow monotonically with T due to the oscillations between the ground state and the low excited states. As the number of atoms increases, the energy gap narrows, so the evolution deviates further from the adiabatic. In such complex evolutions, it is worth carefully investigating the relation between ${\mathcal{J}}_0$ and F₀.

The results are presented in figure 3. Up to thousands of atoms, ${\mathcal{J}}_0$ follows roughly the fidelity of the TF state for various evolution time T. As T increases, when the evolution is nearly adiabatic, ${\mathcal{J}}_0$ follows closer to F₀, as shown in figures 3(a) and 3(b). When the number of atoms continues to increase to N = 11800, both ${\mathcal{J}}_0$ and F₀ are small and only share the similar increasing trend, as illustrated in figure 3(c). The result is not strange, as large N causes a small gap and the evolution deviates further from the adiabatic. As to the conversion efficiency p_c shown in figure 3, this stays near 1, regardless of the diabatic evolution.

Usually, it is challenging to obtain the full MQC spectrum in a spinor BEC, due to lacking the ability to flip the sign of the whole Hamiltonian [6, 7]. Fortunately, as we have seen, low order MQCs manifest important information of many-body dynamics. In fact, ${\mathcal{J}}_0$ can solely catch QPTs, because of its strong relevance to the fidelity in the evolution of the spinor BEC. Therefore a preferable experimental strategy is to focus only on ${\mathcal{J}}_0$. By measuring the number of atoms N₀, and its distribution f(N₀), ${\mathcal{J}}_0={\sum }_{N_0}f{\left(N_0\right)}^2$ can be easily calculated.

To obtain low but non-zero order ${\mathcal{J}}_k$, Lewis-Swan et al provide a pseudo-echo protocol [15]. By replacing the time reversal operator with a pseudo-echo sequence, where the sign of the Hamiltonian is not flipped, and instead,reversing the ramping of q. More explicitly, the time-dependent second-order Zeeman shift during the pseudo-echo is q(t) = q(2T − t) with T ≤ t ≤ 2T. By defining ${\mathcal{J}}_k^{\mathrm{P}\mathrm{E}}=|\mathrm{t}\mathrm{r}[{\hat{\rho }}_k{\hat{\rho }}_{-k}^{\mathrm{P}\mathrm{E}}]|$, where ${\hat{\rho }}^{\mathrm{P}\mathrm{E}}={\mathrm{e}}^{-{\mathrm{i}H}^{\mathrm{P}\mathrm{E}}t}\rho e^{{iH}^{\mathrm{P}\mathrm{E}}t}$, one finds ${\mathcal{J}}_k^{PE}\approx {\mathcal{J}}_k$ (as shown in figure 4). It is easy to prove that, ${\mathcal{J}}_0^{\mathrm{P}\mathrm{E}}$ obtained from such a pseudo-echo is exactly equal to ${\mathcal{J}}_0$.

We proposed a practical scheme to probe the QPT dynamics of many-body spinor BECs, based on MQCs. Numerical simulations of the generation of the TF state in a spin-1 BEC show that MQCs are capable of catching the many-body polar phase and the TF phase. Importantly, only a few low order MQCs are necessary to depict the QPTs. In addition, ${\mathcal{J}}_0$ itself can well characterize the generation of the TF state, due to the similarity between ${\mathcal{J}}_0$ and the fidelity F₀. Finally, although ${\mathcal{J}}_0$ can be calculated from the measurement of particle distribution in spin-1 BECs, non-zero order MQCs may be extracted experimentally with pseudo-echo protocol, where the Hamiltonian is partially time-reversed.

This work is supported by the NSAF under Grant No. U1930201, the National Natural Science Foundation of China (NSFC) under Grant Nos. 12274331, 91836101, ​​​​12135018 and 12204428, and the Innovation Program for Quantum Science and Technology under Grant No. 2021ZD0302100. The numerical calculations in the paper have been partially done on the supercomputing system in the Supercomputing Center of Wuhan University.

Since the energy gap is the smallest around the quantum critical point, it is advantageous to cross the point as slow as possible to reduce excitations. The ramping profile we adopt in numerics is the same as [20], which has the following form (during the $0<t⩽T^{\mathrm{\prime }}$ with $T^{\mathrm{\prime }}=T/2$),The first part (0 ≤ t < τ) gives a slow ramp around q/∣c₂∣ = 2, though the second part ($\tau <t⩽T^{\mathrm{\prime }}$) is linear. The parameter γ is fixed by the condition that the slopes at t = τ in two parts are the same and is given byAnother parameter τ is fixed by the condition $q(t=T^{\mathrm{\prime }})=0$, i.e., $\tau =(4-\alpha )/\gamma +T^{\mathrm{\prime }}$. The rest of the parameters we use are experimental ones [20], α = 2.8, β = 0.16, and ε = 0.956. During the second half of evolution, the ramp profile is centro-symmetric with (A.1), q(t) = − q(T − t) for $T^{\mathrm{\prime }}<t<T$, to guarantee a slow ramp across the second quantum critical point.