Multiparticle Entanglement Dynamics of Quantum Chaos
in a Bose-Einstein condensate

For the OAT model, it is possible to obtain an analytical expression for the QFI [68, 30, 37, 86] at all times and for arbitrary number of particles. For times $ct\leq N\pi/2-2\sqrt{N}$, the QFI is given by

$$\frac{F_{Q}(t)}{N}=1+\frac{(N-1)}{4}\big{[}\alpha_{Q}(t)+\sqrt{\alpha_{Q}^{2}(t)+\beta_{Q}^{2}(t)}\big{]},$$

(11)

where $\alpha_{Q}(t)=1-\cos^{N-2}({2ct}/{N})$ and $\beta_{Q}(t)=4\sin({ct}/{N})\cos^{N-2}({ct}/{N})$. The short-time evolution of the QFI is characterized by a sudden increase starting from $F_{Q}(0)=N$, see Fig. 1. We consider $N\gg 1$, such that $N-1\simeq N-2\simeq N$. A Taylor expansion gives

$$\frac{F_{Q}(t)}{N}=1+ct+\frac{c^{2}t^{2}}{2}+\frac{c^{3}t^{3}}{8}-\frac{3c^{3}t^{3}}{4N}+O\big{(}c^{4}t^{4}\big{)}.$$

(12)

Notice that, up to the second order in $ct$, the Taylor expansion is analogous to that of $e^{ct}$. For times $ct\geq 1$ we have $\alpha_{Q}(t)\gg\beta_{Q}(t)$ and Eq. (11) becomes

$$\frac{F_{Q}(t)}{N}\simeq 1+\frac{N\alpha_{Q}(t)}{2}.$$

(13)

For times $1\leq ct\leq 2\sqrt{N}$ we have $\alpha_{Q}(t)\simeq 2c^{2}t^{2}/N$ such that

$$\frac{F_{Q}(t)}{N}\simeq c^{2}t^{2},$$

(14)

having a quadratic scaling with time. For $ct\geq 2\sqrt{N}$, we have $\cos^{N-2}({2ct}/{N})\ll 1$ and thus obtain $\alpha_{Q}(t)\simeq 1$. Equation (13) thus predicts a constant plateau $F_{Q}(t)=N^{2}/2$ as a function of time, up to $ct\leq N\pi/2-2\sqrt{N}$. Finally, for $ct\geq N\pi/2-2\sqrt{N}$ the QFI further increases and reaches the maximum value $F_{Q}(t)=N^{2}$ at time $ct=N\pi/2$. It should be mentioned that the behavior of QFI given by Eq. (11) is obviously $\pi N$-periodic.

We compare the exact QFI to that calculated within the BMF approach, which can also be solved analytically. For the initial state $s_{x}(0)=1$, the covariance matrix $\mathbf{\Lambda}_{B}$ reads

$$\mathbf{\Lambda}_{B}(t)=\frac{N^{2}}{8}\left(\begin{array}[]{ccc}0&0&0\\ 0&\frac{2+N(1-\cos\frac{2ct}{\sqrt{N}})}{N}&\frac{2\sin\frac{ct}{\sqrt{N}}}{\sqrt{N}}\\ 0&0&\frac{2}{N}\\ \end{array}\right).$$

(15)

The QFI is given by four times the maximum eigenvalue of $\mathbf{\Lambda}_{B}(t)$, namely

$$\frac{F_{B}(t)}{N}=1+\frac{N\alpha_{B}(t)+\sqrt{8N\alpha_{B}(t)+N^{2}\alpha_{B}^{2}(t)}}{4},$$

(16)

where $\alpha_{B}(t)=2\sin^{2}({ct}/{\sqrt{N}})$. Expanding Eq. (16) in a Taylor series for ${ct}/{N}\ll 1$, we obtain

$$\frac{F_{B}(t)}{N}=1+ct+\frac{c^{2}t^{2}}{2}+\frac{c^{3}t^{3}}{8}-\frac{c^{3}t^{3}}{6N}+O\big{(}c^{4}t^{4}\big{)}.$$

(17)

Equation (17) coincides with Eq. (12) up to the second order in $ct$. The difference in the third order is the term $7c^{3}t^{3}/12N$, which is negligible asymptotically for large $N$. Similar considerations can be obtained for higher order terms in the Taylor expansion. For time $ct\geq 1$ we have $\alpha_{B}^{2}(t)\gg\alpha_{B}(t)$ and Eq. (16) becomes

$$\frac{F_{B}(t)}{N}\simeq 1+\frac{N\alpha_{B}(t)}{2},$$

(18)

in analogy with Eq. (13). A Taylor expansion for $ct/\sqrt{N}\gg 1$, predicts the quadratic scaling

$$\frac{F_{B}(t)}{N}\simeq c^{2}t^{2},$$

(19)

for $1\leq ct\leq\sqrt{N}$, in agreement with Eq. (14). For longer times, Eq. (16) is characterized by harmonic oscillations of period $\sqrt{N}\pi/c$. At half period, namely, $t={m\sqrt{N}\pi}/{2c}$ with $m=1,3,\dots$, $\alpha_{B}=2N$, Eq. (16) (wrongly) predicts the Heisenberg limit $F_{B}\simeq N^{2}$, see Fig. 1. We thus conclude that the BMF approach fails for times $ct\geq\sqrt{N}$.

A different approximation to the exact QFI dynamics has been recently studied in Refs. [87, 29] using the Holstein-Primakoff (HP) transformation around the instantaneous mean spin polarization (lying on the $x$ axis, in the present case). Hence, the HP transformation is performed by mapping the transverse ($z$ and $y$) spin components to canonical bosonic variables

$$\hat{J}_{x}=\frac{N}{2}-\hat{n},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \hat{J}_{y}\simeq\sqrt{\frac{N}{2}}\hat{p},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \hat{J}_{z}\simeq\sqrt{\frac{N}{2}}\hat{q},$$

(20)

where the bosonic variables $\hat{q}$ and $\hat{p}$ satisfy $[\hat{q},\hat{p}]=i$, $\hat{n}=({\hat{q}^{2}+\hat{p}^{2}-1})/{2}$. The quantum number $n$ ($=0,1,2,\dots$) labels the quantized spin projection along the $x$ direction of the classical minimum. Using Eqs. (20), the Hamiltonian (1) can be approximated by neglecting terms of order ${1}/{N}$,

$$\hat{H}(t)\simeq a(t)\left(\frac{N}{2}-\hat{n}\right)+c\frac{\hat{q}^{2}}{2},$$

(21)

and thus one can easily obtain the equations of motion for the quantum fluctuations

$$\frac{d\hat{q}}{dt}=-a(t)\hat{p},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \frac{d\hat{p}}{dt}=[a(t)-c]\hat{q}.$$

(22)

This evolution of the quantum fluctuations is exact at the quadratic level and corresponds to the classical one. We start from initial coherent spin state $s_{x}(0)=1$, for which $\langle\hat{q}^{2}(0)\rangle=\langle\hat{p}^{2}(0)\rangle={1}/{2}$ and $\langle\hat{q}(0)\hat{p}(0)\rangle=\langle\hat{p}(0)\hat{q}(0)\rangle=0$. As a result, the QFI can be given by [87, 29]

$$\frac{F_{\rm HP}(t)}{N}=1+2\langle\hat{n}_{\mathrm{exc}}(t)\rangle+2\sqrt{\langle\hat{n}_{\mathrm{exc}}(t)\rangle(\langle\hat{n}_{\mathrm{exc}}(t)\rangle+1)},$$

(23)

where $\langle\hat{n}_{\mathrm{exc}}(t)\rangle=({\langle\hat{q}^{2}\rangle+\langle\hat{p}^{2}\rangle-1})/{2}$ is the number of collective excitations. For the OAT model, $\langle\hat{n}_{\mathrm{exc}}(t)\rangle=c^{2}t^{2}/4$. The Taylor expansion at short times leads to [87, 29, 88]

$$\frac{F_{\rm HP}(t)}{N}=1+ct+\frac{c^{2}t^{2}}{2}+\frac{c^{3}t^{3}}{8}+O(c^{5}t^{5}).$$

(24)

The HP approximation is expected to hold for time sufficiently long as far as the mean spin length is $\big{\langle}\hat{J}_{x}\big{\rangle}\simeq N/2$, namely, as far as we can locally approximate the generalize ($N$-spin) Bloch sphere with the tangent plane [37]. Using the exact result $\big{\langle}\hat{J}_{x}\big{\rangle}=(N/2)\cos^{N-1}(ct)$ [68], we obtain the HP approximation is expected to hold up to times $ct=O(\sqrt{N})$. At intermediate times, for $1\ll ct\ll\sqrt{N}$, we can approximate

$$\frac{F_{\rm HP}(t)}{N}\simeq c^{2}t^{2},$$

(25)

predicting a quadratic scaling of the QFI with time, in agreement with Eqs. (14) and (19).

In Fig. 1(a), we compare the exact $F_{Q}(t)$, Eq. (12), the BMF approximation $F_{B}(t)$, Eq. (17), and the HP approximation, Eq. (24), as a function of time, up to $ct\leq\sqrt{N}$. As discussed above, we see that all models agree very well for short times, see also Refs. [47, 83]. We also see that, as expected, the BMF approach is closer than the HP approximation to the exact solution. This is mainly because the BMF approximation (10) considers both the coupling between quantum fluctuations and the mean spin value and the coupling between fluctuations, while the semiclassical method based on HP transformation only considers the coupling between quantum fluctuations.

The exact dynamics for $F_{Q}(t)$ is not analytical, except for weak interaction [89]. The most interesting situation is that corresponding to the unstable regime [47] (i.e., $c>A$). The quantity $F_{B}(t)$ is obtained by diagonalizing the covariance matrix $\mathbf{\Lambda}_{B}$ on the basis of solving Eq. (10) numerically. For short time, we can obtain an analytical solution when taking $s_{x}(t)\simeq 1-(A^{2}+\lambda^{2})^{2}t^{2}/2NA^{2}$ with $\lambda^{2}=A(c-A)$. The Taylor expansion reads

	$\displaystyle\frac{F_{B}(t)}{N}=$		$\displaystyle 1+ct+\frac{c^{2}t^{2}}{2}+\frac{c\lambda^{2}t^{3}}{6}+\frac{c^{3}t^{3}}{8}-\frac{c^{3}t^{3}}{2N}$		(26)
			$\displaystyle+\frac{c^{2}\lambda^{2}t^{4}}{6}-\frac{c^{4}t^{4}}{2N}+O(c^{5}t^{5}),$

that agrees with Ref. [47]. Similarly to the OAT [83], the Taylor expansion up to the second order is analogous to that of $e^{ct}$. Nevertheless, the third order of Eq. (26) is larger than that of Eq. (17) by a factor ${\lambda^{2}}/{c^{2}}={A}(1-{A}/{c})/{c}$, that has a maximum at $c/A=2$ [47]. We can thus conclude that the BMF approximation breaks down for times of the order of

$$t_{c}\simeq\frac{1}{\lambda}\log\bigg{(}\frac{N\lambda^{2}}{cA}\bigg{)},$$

(27)

where $t_{c}$ has been identified in Ref. [77] as the time needed to generate a macroscopic superposition state.

We can compare this prediction with the semiclassical HP approach of Ref. [29]. The analytical solution of $F_{\rm HP}(t)$ is still given by Eq. (23) where, in the unstable regime, $\langle\hat{n}_{\mathrm{exc}}(t)\rangle=({A}/{\lambda}+{\lambda}/{A})^{2}\sinh^{2}(\lambda t)/{4}$ [88]. The Taylor expansion of $F_{\rm HP}(t)$ at short times reads [87, 29, 88]

$$\frac{F_{\rm HP}(t)}{N}=1+ct+\frac{c^{2}t^{2}}{2}+\frac{c\lambda^{2}t^{3}}{6}+\frac{c^{3}t^{3}}{8}+\frac{c^{2}\lambda^{2}t^{4}}{6}+O(c^{5}t^{5}).$$

(28)

Similarly to the OAT case, the agreement between Eqs. (26) and (28) is up to the second order of $ct$. The difference in the third order term is $c^{3}t^{3}/2N$, which is negligible asymptotically for large $N$. Similar considerations can be obtained for higher order terms in the Taylor expansion. Furthermore, neglecting subleading terms for $N\gg 1$, Eqs. (26) and (28) agree with the results of Ref. [47]. Finally, in the limit $\langle\hat{n}_{\mathrm{exc}}(t)\rangle\gg 1$, that is obtained after a transient time, $c/\lambda\leq ct\leq\log(N\lambda^{2}/cA)$, the HP approach predicts [29, 88]

$$\frac{F_{\rm HP}(t)}{N}\simeq\bigg{(}\frac{A}{\lambda}+\frac{\lambda}{A}\bigg{)}e^{2\lambda t},$$

(29)

where the largest instability parameter $\lambda/c$ is obtained for $c/A=2$. A comparison between the exact $F_{Q}(t)$ and the approximations $F_{B}(t)$ and $F_{\rm HP}(t)$ is shown in Fig. 1(b).

Here, the total number $N$ is conserved and $1/N$ can be identified with an effective Planck constant. The mean-field (MF) limit $N\rightarrow\infty$ (and thus the interaction $c/N\rightarrow 0$) corresponds to the classical limit, where $(s_{x},s_{y},s_{z})$ become classical variables that form a single-particle Bloch vector confined to the unit sphere and thus the norm is conserved [48]. $(s_{z},\phi)$ can be identified with classical phase-space variables, where $\phi$ is the azimuthal angle. The first-order truncation of the hierarchy equations together with the approximations $\langle\hat{J}_{l}\hat{J}_{k}\rangle\simeq\langle\hat{J}_{l}\rangle\langle\hat{J}_{k}\rangle$ give the equations of motion of the expectation values for the first-order moments, namely, the MF equation [48]

$$\frac{d}{dt}\left(\begin{array}[]{c}s_{x}\\ s_{y}\\ s_{z}\\ \end{array}\right)=\left(\begin{array}[]{ccc}0&-cs_{z}&0\\ cs_{z}&0&-a(t)\\ 0&a(t)&0\\ \end{array}\right)\left(\begin{array}[]{c}s_{x}\\ s_{y}\\ s_{z}\\ \end{array}\right),$$

(35)

where $a(t)$ is given by Eq. (30). For such a periodic-driven system, the classical motion can usually be divided into two classes: i) stable and periodic motion and ii) unstable and chaotic motion. In the unstable case, the magnification of the initially small deviation is exponential in time, which can be measured by the maximum Lyapunov exponent, i.e., [90]

$$\lambda_{L}(A,c)=\lim_{m\rightarrow\infty}\frac{1}{mT}\sum_{n=1}^{m}\ln\frac{||(\delta s_{x}(nT),\delta s_{y}(nT),\delta s_{z}(nT))||}{||(\delta s_{x}(0),\delta s_{y}(0),\delta s_{z}(0))||},$$

(36)

where $\delta s_{x,y,z}(0)=|s^{\prime}_{x,y,z}(0)-s_{x,y,z}(0)|$ denote the initial deviations and $\delta s_{x,y,z}(nT)=|s^{\prime}_{x,y,z}(nT)-s_{x,y,z}(nT)|$ denote the deviations at $t=nT$, where $s_{x,y,z}(t)$ and $s^{\prime}_{x,y,z}(t)$ respectively represent the two classical trajectories without perturbation and with perturbation initially. In Fig. 2(a) we plot the maximum Lyapunov exponents in the $(A,c)$-plane for the initial state $s_{x}(0)=1$. In our numerical calculations we set $\delta s_{x}(0)=10^{-5}$ and evolve the system for 500 periods. For the stable motion, $\lambda_{L}$ is zero (dark blue region), while for the chaotic regime $\lambda_{L}$ is positive. For example, in Fig. 2(a) for $A/\pi=0.4$ and $c/\pi=0.2$, the initial state is in stable region, while $A/\pi=0.4$ both for $c/\pi=0.8$ and $c/\pi=1.4$ the initial states are in chaotic region. However, the point $(A,c)/\pi=(0.4,1.4)$ is “more chaotic” than $(A,c)/\pi=(0.4,0.8)$ is, because $(A,c)/\pi=(0.4,1.4)$ ($\lambda_{L}\simeq 1.524$) has a larger Lyapunov exponent than $(A,c)/\pi=(0.4,0.8)$ ($\lambda_{L}\simeq 0.966$), meaning that $\delta s_{x,y,z}(nT)$ deviate faster from the initial conditions.

To show the property of classical motion of the kicked system, we plot the Poincaré sections by using the stroboscopic points at $t=nT$ in Figs. 2(b)-2(d) for different interaction strengths $c/\pi=0.2$, $0.8$, and $1.4$ with same kicked strength $A/\pi=0.4$. We see that for the integrable situation [e.g., $(A,c)/\pi=(0.4,0.2)$, see Fig. 2(b)], all trajectories (periodic orbit and elliptic fixed point) are stable with the maximum Lyapunov exponents being zero. In the completely nonintegrable case [e.g., $(A,c)/\pi=(0.4,1.4)$, see Fig. 2(d)], we see that the phase space is full of chaotic trajectories with largely positive $\lambda_{L}$. In the mixed case [e.g., $(A,c)/\pi=(0.4,0.8)$, see Fig. 2(c)], the phase space illustrates one big-size island, six middle-size islands, and four small-size islands (see the blank regions). It is clear that the motions inside the islands are stable and the motions outside the islands are unstable and mainly chaotic with medium positive $\lambda_{L}$. In Figs. 2(b)-2(d), the green points mark the initial state $s_{x}(0)=1$.

To show the influence of the chaotic motions on the dynamics of entanglement, we calculate the QFI for the same initial state [$s_{x}(0)=1$] and different interacting parameters. In Fig. 3, for a fixed value of the parameter $c$, we compare the time-evolution of the QFI between the QKR model (solid black line; obtained for parameter $A$ with $\delta$-kicked time dependence), the OAT model (dashed red line; $A=0$), and the TAT model (dot-dashed green line; for constant $A$). When the initial quantum state corresponds to a stable point (i.e., $c/\pi=0.2$) or a saddle point (i.e., $c/\pi=0.5$) in the Poincaré section of the classical model, the QKR model does not show any advantage in entanglement generation because of the small oscillation [see Fig. 3(a)] or the slow growth [see Fig. 3(b)] of the QFI with time. However, when the initial quantum state is located in the classical chaotic region (i.e., $c/\pi=0.8$, and 2.0), compared with the QFI for the OAT model, the QFI for the QKR model can reach and exceed ${N^{2}}/{2}$ in a very short period of time (e.g., 2 or 3 kicks), which indicates that the quantum chaos can significantly accelerate the generation of entanglement [see Fig. 3(c) and 3(d)]. When compared with the TAT model, the QFI for the QKR model is slower in all cases. In particular, it was shown in Ref. [47] that the maximally unstable TAT dynamics (for $c/A=2$) is characterized by the fastest generation of entanglement within quadratic Hamiltonians as in Eq. (1).

In Fig. 4, we compare the exact $F_{Q}(t)$ with both the BMF approximation $F_{B}(t)$ and the HP approximation $F_{\rm HP}(t)$ for the periodic-kicked system with the initial state $s_{x}(0)=1$. When the initial state is a classical stable point, the three QFI agree well for a long time [see Fig. 4(a)]. When the initial state is an unstable saddle point [see Fig. 4(b)] or in the classical chaotic region [see Figs. 4(c) and 4(d)] the three results remain consistent only for a short time (a few kicks in the figure). While the $F_{\rm HP}(t)$ diverges after a transient time, the BMF follows more closely the exact results for longer times. It is thus interesting to introduce the criterion

$$\frac{|F_{Q}(t_{IB})-F_{B}(t_{IB})|}{F_{Q}(t_{IB})}=g_{IB},$$

(37)

to quantitatively define a break time $t_{IB}$ where the BMF approximation departs from the exact result (by a relative amount $g_{IB}$). This time is closely related to the instability degree of the initial state, which we will relate below to the maximum Lyapunov exponent. In practice, we set $g_{IB}=0.01$ in the following. The main result of this paper is to study numerically the break time $t_{IB}$ in different chaotic regimes and determine characteristic behaviours. We consider different regimes:

• $\bullet$

  If the initial condition is in the classical stable region, we find

  $$t_{IB}\simeq\alpha_{\rm st}(A,c)\times\sqrt{N}.$$

  (38)

  In particular, we can demonstrate analytically that the $O(\sqrt{N})$ behaviour holds for the OAT dynamics. To support the claim (38) we show numerical results in Fig. 5(a), obtained for the QKR model for $A/\pi=0.4$ and $c/\pi=0.2$. There, we plot the break time $t_{IB}$ as a function of the number of particles $N$. The dashed line in Fig. 5(a) is $\sqrt{N}/4.7$ fitting the OAT results, while the solid line is $\sqrt{N}/15$ fitting the QKR results for large $N$, in agreement with Eq. (38).

• $\bullet$

  If the initial condition is at the saddle point, we find

  $$t_{IB}\simeq\alpha_{\rm sp}(A,c)\times\ln{N}.$$

  (39)

  In particular, this behaviour holds for the TAT dynamics as well as in the QKR model at the saddle point [see Fig. 5(b)]. The dashed line in Fig. 5(b) is $\ln N/2.02$ fitting the TAT results, while the solid line is $\ln N/2.2$ fitting the QKR results for large $N$, in agreement with Eq. (39). Notice that, for the TAT model with initial conditions at the saddle point, we expect $t_{IB}\propto t_{c}$, where $t_{c}\sim(1/\lambda)\log N$ (for $N\gg 1$), as given in Eq. (27), with $\lambda=\sqrt{A(c-A)}$.

• $\bullet$

  If the initial conditions are in the deep chaotic regime we find

  $$t_{IB}\simeq\alpha_{\rm ch}(A,c)\times(\ln{N})^{4},$$

  (40)

  see Fig. 5(c). The solid line in Fig. 5(c) is the fit $(\ln N)^{4}/1590$, the dashed line is $(\ln N)^{4}/2810$, the dash-dotted line is $(\ln N)^{4}/3200$, and the dash-double-dotted line is $(\ln N)^{4}/4000$. In the chaotic regime, all results agree well with Eq. (40). Notice that we numerically do not see a smooth transition between the $O(\log N)$ behaviour at the saddle point and the $O((\log N)^{4})$ behaviour in the deep chaotic regime, but rather a sharp transition. Furthermore, we can connect the prefactor $\alpha_{\rm ch}(A,c)$ in Eq. (40) and the asymptotically maximum Lyapunov exponent $\lambda_{L}(A,c)$ (depending on both parameters $c$ and $A$). For strong chaos and large $N$, we find numerically

  $$\alpha(A,c)_{\rm ch}\simeq\eta(A){e}^{-\gamma(A)\lambda_{L}(A,c)},$$

  (41)

  where the coefficients $\gamma(A)$ and $\eta(A)$ depend only on the kick strength $A$. A numerical investigation of the Eq. (41) is shown in Fig. 6. In Fig. 6(a) we show the break time as a function of $\lambda_{L}(A,c)$ and for different values of kicked strength, $A/\pi=0.2$, $0.3$, and $0.4$. Equation (41) predicts that, while fixing $A$, the break time $t_{IB}$ depends on $c$ through an exponential function of the maximum Lyapunov exponent. This behaviour (lines) is well reproduced by the numerical calculations (symbols) for sufficiently large values of $\lambda_{L}(A,c)$, as shown in Fig. 6(a). Furthermore, when fixing the parameter $c$, Fig. 2(a) shows that the maximum Lyapunov exponent in the chaotic region is a double-valued function of $A$. Therefore, Eq. (41) also predicts that $t_{IB}$ is double-valued function of $A$, for fixed $c$. This prediction is confirmed in Fig. 6(b), where we plot $t_{IB}$ as a function of $\lambda_{L}(A,c)$ for different values of $c$.