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

Sheng-Chang Li1,2 Luca Pezzè2 Augusto Smerzi2 1\. MOE Key Laboratory for
Nonequilibrium Synthesis and Modulation of Condensed Matter, Shaanxi Province
Key Laboratory of Quantum Information and Quantum Optoelectronic Devices, and
School of Physics, Xi’an Jiaotong University, Xi’an 710049, China
2\. QSTAR, INO-CNR and LENS, Largo Enrico Fermi 2, Firenze 50125, Italy

###### Abstract

We study the particle-entanglement dynamics witnessed by the quantum Fisher
information (QFI) of a trapped Bose-Einstein condensate governed by the kicked
rotor Hamiltonian. The dynamics is investigated with a beyond mean-field
approach. We link the time scales of the validity of this approximation in,
both, classical regular and chaotic regions, with the maximum Lyapunov
exponents of the classical system. This establishes an effective connection
between the classical chaos and the QFI. We finally study the critical point
of a quantum phase transition using the beyond mean-field approximation by
considering a two-mode bosonic Josephson junction with attractive
interparticle interaction.

##  I Introduction

Quantum chaos aims to investigate the properties of quantum systems that have
a chaotic dynamics in the classical limit. First works on quantum chaos can be
traced back to 1950s [1] or even earlier [2]. While systematic theoretical
studies emerged in the late 1980s [3], the experimental activity – requiring
highly controllable many-body systems – has begun in the last years [4].
Atomic Bose-Einstein condensates (BECs) can be trapped in time-dependent
potentials having almost arbitrary geometries and are therefore emerging as
ideal candidates for the study of quantum chaos in many-body systems [5, 6, 7,
8, 9, 10]. Several efforts have been devoted to investigate chaotic dynamics
in both scalar and spinor BECs, including driven two-mode [11], tilted triple-
well [12], kicked ring [13], kicked optical-lattice [14], barrier-height-
modulated double-well [15], dissipative double-well [16], dissipative optical-
lattice [3], spin-1 [17], spin-2 [18], and coupled spinor atom-molecule BECs
[19]. A paradigmatic model to study quantum chaos is the kicked rotor [20, 21]
which can be realized with trapped BECs [22, 23, 24, 25, 26]. This model has
been widely studied to illustrate classical and quantum chaos and the
transition between them [27, 28].

An intriguing problem emerging in the recent literature is the study of the
interplay between quantum chaos and entanglement. It has been indeed
recognized that chaotic dynamics and entanglement generation are intimately
related in many systems [31, 30, 29]. In particular, Ref. [31] presented an
unifying picture connecting the multiparticle entanglement growth to the
quantifiers of classical and quantum chaos. In this paper, we study the
multiparticle entanglement dynamics of a BEC governed by the quantum kicked
rotor (QKR) model. Multipartite entanglement is witnessed by studying the
quantum Fisher information (QFI) [32, 33, 34, 35, 36, 37]. The QFI [38, 39]
has proved a quite valuable quantity for the analysis of entanglement in
complex systems [31, 30, 40, 47, 41, 42, 43, 44, 45, 46]. Here, we study the
chaotic dynamics of the QFI under the framework of a beyond mean-field (BMF)
approximation [48]. We show that the time scalings of the validity of the BMF
dynamics are related with the maximum Lyapunov exponents of the corresponding
classical chaotic regimes. This provides an effective relation between the
time scales of quantum evolution in a chaotic system and the QFI. We further
extend the BMF method to the analysis of the ground state of the system. We
show that the QFI evaluated within the BMF method is in excellent agreement
with exact results able to characterize the quantum phase transition in
finite-size systems. The multipartite entanglement detected by the QFI is a
resource for parameter estimation beyond the sensitivity limit achievable with
separable (or classically-correlated) states [32, 33, 34, 35, 36, 37]. In this
respect, our results show that the generation of metrologically-useful
entanglement can be accelerated in chaotic systems, when compared to
entanglement generation in integrable systems. As shown in Ref. [49], quantum
chaos can also enhance the sensitivity in the estimation of Hamiltonian
parameters.

##  II Quantum Fisher information within a beyond mean-field approach

In the following, we study the dynamics of a two-mode model governed by the
Hamiltonian ($\hbar=1$)

$\hat{H}(t)=a(t)\hat{J}_{x}+\frac{c}{N}\hat{J}_{z}^{2},$ (1)

where
$\hat{J}_{x}=(\hat{a}_{1}^{\dagger}\hat{a}_{2}+\hat{a}_{2}^{\dagger}\hat{a}_{1})/2$,
$\hat{J}_{y}=(\hat{a}_{1}^{\dagger}\hat{a}_{2}-\hat{a}_{2}^{\dagger}\hat{a}_{1})/2{i}$,
and
$\hat{J}_{z}=(\hat{a}_{1}^{\dagger}\hat{a}_{1}-\hat{a}_{2}^{\dagger}\hat{a}_{2})/2$
are the SU(2) generators satisfying angular momentum commutation relation
$[\hat{J}_{x},\hat{J}_{y}]=i\hat{J}_{z}$ with cyclic permutation. Here,
$\hat{a}_{1,2}$ ($\hat{a}_{1,2}^{\dagger}$) are bosonic annihilation
(creation) operators for the two modes (labelled as “1” and “2”). The time-
dependent parameter $a(t)$ describes the coupling between the two modes, $c$
depends on the two-body interaction, and $N$ is the total number of particles.
In atomic BEC experiments, the two-mode model (1) can be realized by using
either a double-well trapping potential [50, 51, 52, 53, 54, 56, 57, 55, 58,
59, 60, 61] or two hyperfine atomic states [62, 63, 64, 65]. As initial state,
we consider a coherent spin state [66]

$\displaystyle|\psi_{s}\rangle=$
$\displaystyle\sum_{\mu=-\frac{N}{2}}^{\frac{N}{2}}\sqrt{\bigg{(}\begin{array}[]{c}N\\\
\frac{N}{2}+\mu\end{array}\bigg{)}}\bigg{(}\cos\frac{\theta}{2}\bigg{)}^{\frac{N}{2}-\mu}\bigg{(}\sin\frac{\theta}{2}\bigg{)}^{\frac{N}{2}+\mu}$
(5) $\displaystyle\times e^{i(\frac{N}{2}+\mu)\phi}|\mu\rangle,$

where $|\mu\rangle$ are eigenstates of $\hat{J}_{z}$ and
$s=(\sin\theta\cos\phi;\sin\theta\sin\phi;\cos\theta)$ is the mean spin
direction. We will focus on the dynamical evolution of the initial state
$s_{x}=1$ ($s_{y}=s_{z}=0$). This corresponds to the coherent spin state given
with all $N$ spins pointing along the $+x$ direction in the Bloch sphere. For
the initial state, the collective spin mean values and variances are
$\langle\hat{J}_{x}\rangle={N}/{2}$ and $\langle\hat{J}_{\perp}\rangle=0$,
$(\Delta\hat{J}_{x})^{2}=0$ and $(\Delta\hat{J}_{\perp})^{2}={N}/{4}$,
respectively, where $\perp$ indicates an arbitrary direction orthogonal to the
mean spin direction.

We characterize the quantum dynamics by studying the QFI [38, 39, 35, 37]. The
QFI is a property of a generic state $\hat{\rho}$ with respect to a chosen
operator $\hat{J}$ and is formally defined as the variance
$F_{Q}[\hat{\rho},\hat{J}_{n}]=(\Delta\hat{L})^{2}$ of the symmetric
logarithmic derivative $\hat{L}$, satisfying
$2[\hat{J},\hat{\rho}]=\\{\hat{L},\hat{\rho}\\}$ [38, 39, 35, 37]. For a pure
quantum state $|\psi_{t}\rangle$ at time $t$, the QFI can be expressed as

$F_{Q}(t)=4(\Delta\hat{J})^{2},$ (6)

where
$(\Delta\hat{J})^{2}=\langle\psi_{t}|\hat{J}^{2}|\psi_{t}\rangle-\langle\psi_{t}|\hat{J}|\psi_{t}\rangle^{2}$.
In the following we will restrict to the calculations of the QFI for the
operators $\hat{J}=\hat{J}_{x,y,z}$ introduced above. In particular, the
optimization of the QFI over $\hat{J}_{x,y,z}$ can be obtained as the maximum
eigenvalue of the $3\times 3$ covariance matrix

$[\mathbf{\Lambda}_{Q}]_{kl}=2[(\langle\hat{J}_{k}\hat{J}_{l}\rangle+\langle\hat{J}_{l}\hat{J}_{k}\rangle)-2\langle\hat{J}_{k}\rangle\langle\hat{J}_{l}\rangle],$
(7)

with $k,l=x,y,z$, and the optimal direction is given by the corresponding
eigenvector.

The QFI is related to multipartite entanglement [32, 33, 34, 35, 36, 37]: a
QFI satisfying $F_{Q}(t)>\mathscr{s}k^{2}+\mathscr{r}^{2}$, where
$\mathscr{s}=\left\lfloor{N}/{k}\right\rfloor$ (being the maximum integer
smaller than or equal to ${N}/{k}$) and $\mathscr{r}=N-\mathscr{s}k$, is a
witness of $(k+1)$-partite entanglement among the $N$ bosonic particles in the
state $|\psi_{t}\rangle$. Multipartite entanglement generated the dynamical
evolution of an initial coherent spin state via Eq. (1) has been largely
studied in the literature. Theoretical works [32, 67] have focused on the one-
axis twisting (OAT) Hamiltonian, corresponding to $a(t)=0$ in Eq. (1) [68,
69]. The OAT dynamics has been realized experimentally by exploiting
interatomic interaction in BECs [65, 70, 71] or the atom-light interactions in
cold thermal atoms [73, 72, 74, 75, 76], see Ref. [37] for a comprehensive
review. Another model of interest is the twist-and-turn (TAT) Hamiltonian,
corresponding to $a(t)=A>0$ in Eq. (1) [51, 62, 63, 77]. The TAT dynamics is
characterized by the competition between the twisting (due to the quadratic
term in the Hamiltonian) and the linear coupling between the two modes that
turns the state along an orthogonal direction. It has been shown, both
experimentally [78, 79] and theoretically [47] that the TAT model creates
entanglement faster than the OAT model. More recently, multipartite
entanglement witnessed by the QFI has been studied for chaotic systems, where
$a(t)$ has a $\delta$-function time dependence [31].

While exact numerical calculations of the QFI for the model (1) can be
obtained for long times and large systems, it is interesting to study various
approximations to the quantum dynamics of the QFI. For instance, Ref. [47] has
considered an exact quantum phase model, Refs. [80, 81] have used the a Fock
state expansion [54, 82] for the study of the spin-squeezing parameter, which
is a saturable lower bound to the QFI [37], while Ref. [31] has proposed a
semiclassical approximation. In this paper, we study the time evolution of the
QFI within a BMF approach [48].

Here, rather than using a Schrödinger picture, we consider the dynamics of the
spin operators. According to the Heisenberg equation for the operator
$\hat{J}_{l}$, $i({d\hat{J}_{l}}/{dt})=[\hat{J}_{l},\hat{H}(t)]$, the
expectation value of $\hat{J}_{l}$ depends not only on themselves, but also on
the second-order moments $\langle\hat{J}_{l}\hat{J}_{k}\rangle$. Similarly,
the time evolution of the second-order moments depends on third-order moments,
and so on. Consequently, we can obtain the hierarchy of equations of motion
for the expectation values $\langle\hat{J}_{l}\rangle$,
$\langle\hat{J}_{l}\hat{J}_{k}\rangle$, $\dots$. The second-order truncation
of the hierarchy equations and the approximation
$\langle\hat{J}_{l}\hat{J}_{k}\hat{J}_{j}\rangle\simeq\langle\hat{J}_{l}\hat{J}_{k}\rangle\langle\hat{J}_{j}\rangle+\langle\hat{J}_{l}\rangle\langle\hat{J}_{k}\hat{J}_{j}\rangle+\langle\hat{J}_{l}\hat{J}_{j}\rangle\langle\hat{J}_{k}\rangle-2\langle\hat{J}_{l}\rangle\langle\hat{J}_{k}\rangle\langle\hat{J}_{j}\rangle$
will give the equations of motion of the expectation values for the first- and
second-order moments, which are [48]

$\frac{d}{dt}\left(\begin{array}[]{c}s_{x}\\\ s_{y}\\\ s_{z}\\\ \Delta_{xz}\\\
\Delta_{yz}\\\ \Delta_{xy}\\\ \Delta_{xx}\\\ \Delta_{yy}\\\ \Delta_{zz}\\\
\end{array}\right)=\left(\begin{array}[]{ccccccccc}0&-cs_{z}&0&0&-\frac{c}{2}&0&0&0&0\\\
cs_{z}&0&-a(t)&\frac{c}{2}&0&0&0&0&0\\\ 0&a(t)&0&0&0&0&0&0&0\\\
0&0&0&0&-cs_{z}&a(t)&0&0&-cs_{y}\\\ 0&0&0&cs_{z}&0&0&0&a(t)&cs_{x}-a(t)\\\
0&0&0&cs_{x}-a(t)&-cs_{y}&0&cs_{z}&-cs_{z}&0\\\
0&0&0&-2cs_{y}&0&-2cs_{z}&0&0&0\\\ 0&0&0&0&2(cs_{x}-a(t))&2cs_{z}&0&0&0\\\
0&0&0&0&2a(t)&0&0&0&0\\\ \end{array}\right)\left(\begin{array}[]{c}s_{x}\\\
s_{y}\\\ s_{z}\\\ \Delta_{xz}\\\ \Delta_{yz}\\\ \Delta_{xy}\\\ \Delta_{xx}\\\
\Delta_{yy}\\\ \Delta_{zz}\\\ \end{array}\right),$ (8)

where
$(s_{x},s_{y},s_{z})={2}(\langle\hat{J}_{x}\rangle,\langle\hat{J}_{y}\rangle,\langle\hat{J}_{z}\rangle)/{N}$
and
$\Delta_{lk}={4}(\langle\hat{J}_{l}\hat{J}_{k}+\hat{J}_{k}\hat{J}_{l}\rangle-2\langle\hat{J}_{l}\rangle\langle\hat{J}_{k}\rangle)/{N^{2}}$.
This BMF model is also referred to as the Bogoliubov backreaction model [48,
83]. Taking ${N^{2}}\Delta_{lk}/{8}$ as the matrix elements can also form a
similar covariance matrix $\mathbf{\Lambda}_{B}$ with
$\Delta_{lk}=\Delta_{kl}$ and four times the maximum eigenvalue gives the QFI
within the BMF approximation, namely, $F_{B}$.

##  III Benchmark models

While the main focus of this paper is on the QKR model, which is discussed in
the next section, we first analyze here the QFI calculated within the BMF
approach for the OAT and TAT models.

###  III.1 One-axis-twisting dynamics

For the OAT model, it is possible to obtain an analytical expression for the
QFI [68, 32, 37, 84] at all times and for arbitrary number of particles. For
times $ct\lesssim 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{]},$
(9)

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{)}.$
(10)

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

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

For times $1\lesssim ct\lesssim 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},$ (12)

having a quadratic scaling with time. For $ct\gtrsim 2\sqrt{N}$, we have
$\cos^{N-2}({2ct}/{N})\ll 1$ and thus obtain $\alpha_{Q}(t)\simeq 1$. Equation
(11) thus predicts a constant plateau $F_{Q}(t)=N^{2}/2$ as a function of
time, up to $ct\lesssim N\pi/2-2\sqrt{N}$. Finally, for $ct\gtrsim
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$.

We compare the exact QFI to that calculated within the BMF approach, that 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).$ (13)

The QFI is given by the 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},$
(14)

where $\alpha_{B}(t)=2\sin^{2}({ct}/{\sqrt{N}})$. Expanding Eq. (14) 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{)}.$
(15)

Equation (15) coincides with Eq. (10) 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\gtrsim
1$ we have $\alpha_{B}^{2}(t)\gg\alpha_{B}(t)$ and Eq. (14) becomes

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

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

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

for $1\lesssim ct\lesssim\sqrt{N}$, in agreement with Eq. (12). For longer
times, Eq. (14) 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. (14) (wrongly) predicts the Heisenberg
limit $F_{B}\simeq N^{2}$, see Fig. 1. We thus conclude that the BMF approach
fails for times $ct\gtrsim\sqrt{N}$.

A different approximation to the exact QFI dynamics has been recently studied
in Refs. [85, 31] 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},$ (18)

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. (18), 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},$ (19)

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}.$ (20)

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 [85, 31]

$\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)},$
(21)

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 [85, 31, 86]

$\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}).$ (22)

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},$ (23)

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

In Fig. 1(a), we compare the exact $F_{Q}(t)$, Eq. (10), the BMF approximation
$F_{B}(t)$, Eq. (15), and the HP approximation, Eq. (22), as a function of
time, up to $ct\lesssim\sqrt{N}$. As discussed above, we see that all models
agree very well for short times, see also Refs. [47, 81]. 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 (8) 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.

Figure 1: QFI dynamics for (a) the OAT model [$a(t)=0$ in Eq. (1)] and (b) the
TAT model [for $A=\frac{\pi}{2}$]. The solid black line is the exact solution
(see text), the dashed line is the BMF approach $F_{B}(t)$, while the $F_{\rm
HP}$ is the semiclassical HP Eq. (21). We have set $c=\pi$, $N=400$ and an
initial coherent spin state pointing along the $x$ axis ($s_{x}=1$).

###  III.2 Twist-and-turn dynamics

The exact dynamics for $F_{Q}(t)$ is not analytical, except for weak
interaction [87]. 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. (8) 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}$
(24)
$\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 [81], the Taylor expansion up
to the second order is analogous to that of $e^{ct}$. Nevertheless, the third
order of Eq. (24) is larger than that of Eq. (15) 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{)},$
(25)

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.
[31]. The analytical solution of $F_{\rm HP}(t)$ is still given by Eq. (21)
where, in the unstable regime,
$\langle\hat{n}_{\mathrm{exc}}(t)\rangle=({A}/{\lambda}+{\lambda}/{A})^{2}\sinh^{2}(\lambda
t)/{4}$ [86]. The Taylor expansion of $F_{\rm HP}(t)$ at short times reads
[85, 31, 86]

$\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}).$
(26)

Similarly to the OAT case, the agreement between Eqs. (24) and (26) 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. (24) and (26)
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\lesssim ct\lesssim\log(N\lambda^{2}/cA)$, the HP
approach predicts [31, 86]

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

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).

##  IV Quantum kicked rotor model

In the following we consider the case

$a(t)=\frac{AT}{\tau_{1}}\sum_{n=0}^{\infty}\delta(t-\tau_{n}),$ (28)

such that the two-mode time-dependent Hamiltonian (1) becomes

$\hat{H}(t)=\frac{AT}{\tau_{1}}\sum_{n=0}^{\infty}\delta(t-\tau_{n})\hat{J}_{x}+\frac{c}{N}\hat{J}_{z}^{2},$
(29)

where $A$ is the strength of the kicked field. The Dirac $\delta$-function
presents a rectangular-pulse-type kicked field with
$\tau_{n}\in[nT,nT+\tau_{1}]$. Here, $T=\tau_{0}+\tau_{1}$ is the period of
the kicked field with $\tau_{0}$ and $\tau_{1}$ being the interval and the
width of the kicks, respectively. Equation (29) corresponds to the QKR model
[20, 21]. The periodicity of the two-mode QKR model (29) allows us to describe
stroboscopically the evolution of quantum state with discrete time in units of
$T$, i.e.,

$|\psi(n)\rangle=U_{R}(T)|\psi(n-1)\rangle=U_{R}^{n}(T)|\psi(0)\rangle,$ (30)

with the unitary Floquet operator

$U_{R}(T)={e}^{-i[\frac{c}{N}\hat{J}_{z}^{2}]\tau_{0}}{e}^{-i[\frac{AT}{\tau_{1}}\hat{J}_{x}+\frac{c}{N}\hat{J}_{z}^{2}]\tau_{1}}.$
(31)

When $\tau_{1}\rightarrow 0$, $\tau_{0}\rightarrow T$, the kicked field
becomes a $\delta$-function-type form and the Floquet operator reduces to

$U_{\delta}(T)={e}^{-i[\frac{c}{N}\hat{J}_{z}^{2}]T}{e}^{-i[{A}\hat{J}_{x}]T}.$
(32)

In our calculations, we adopt a rectangle-pulse-type kicked field by setting
$\tau_{0}=1$ and $\tau_{1}=0.01$, which is very close to the
$\delta$-function-type kicked field.

Figure 2: (a) Maximum Lyapunov exponents for the initial state $s_{x}(0)=1$
[i.e., $(s_{z},\phi)=(0,0)$]. (b), (c), and (d) show the Poincaré sections
with $A/\pi=0.4$ for $c/\pi=0.2$, $0.8$, and $1.4$, respectively. Green solid
points indicate the initial states adopted in (a).

###  IV.1 Mean-field approximation and chaotic motion

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),$ (33)

where $a(t)$ is given by Eq. (28). 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., [88]

$\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))||},$ (34)

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$.

###  IV.2 Quantum Fisher Information and maximum Lyapunov exponent

Figure 3: Time evolution of QFI for the QKR ($A/\pi=0.4$, black solid lines),
the OAT ($A=0$, red dashed lines), and the TAT ($A/\pi=0.4$, green dash-dotted
lines) models with different interaction strengths: (a) $c/\pi=0.2$ (the
initial state being a stable point), (b) $c/\pi=0.5$ (the initial state being
a saddle point), (c) $c/\pi=0.8$ [The initial state is in chaotic regime, see
Fig. 2(c)], and (d) $c/\pi=2.0$ (The initial state is in deep chaotic regime).
We have set $N=200$ and $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).

Figure 4: Comparison of QFI dynamics for the QKR model with the initial state
$s_{x}(0)=1$ being (a) a stable fixed point (i.e., $c/\pi=0.2$), (b) a saddle
point (i.e., $c/\pi=0.5$), (c) in chaotic regime (i.e., $c/\pi=0.8$), and (d)
in deep chaotic regime (i.e., $c/\pi=1.4$). The solid black line is the exact
solution $F_{Q}(t)$, the dashed line is the BMF approach $F_{B}(t)$, while the
$F_{\rm HP}$ is the semiclassical HP Eq. (21). We have set $A/\pi=0.4$ and
$N=400$. Figure 5: Top [Bottom] row panels show the break time $t_{IB}$, as
defined in Eq. (35) [$t_{\rm HP}$, Eq. (40)], as a function of the total
particle number $N$ for the QKR model. The different panels corresponds to the
initial state being (a,d) a classical stable point, (b,e) a classical saddle
point, and (c,f) in the deep chaotic region. In all panels, symbols are
numerical results while lines are fits, see text. Values of Hamiltonian
parameters are reported in each panel, while parameters of the parameters of
the fitting lines are reoported in the text. Panel (a,d) shows the results for
the OAT model (black diamonds) and the QKR model in the stable region (red
squares). Panel (b,e) shows the results for the TAT model (black diamonds) and
the QKR model at the saddle point (red squares). Panel (c,f) shows the results
of the QKR model in the chaotic regime.  Figure 6: Break time $t_{IB}$ as a
function of the maximum Lyapunov exponent $\lambda_{L}$. Symbols are numerical
results, lines are fits according to Eq. (39). In panel (a), the parameter $A$
is fixed and the Lyapunov exponent only depends on $c$. The fitting lines are
$2.83\times{e}^{-1.46\lambda_{L}}$ (solid line),
$2.38\times{e}^{-1.22\lambda_{L}}$ (dash line), and
$1.82\times{e}^{-1.03\lambda_{L}}$ (dash-dot line), respectively. In panel
(b), the parameter $c$ is fixed and the maximum Lyapunov exponent is a double-
valued function of $A$. The display of the lines between the scattered points
is only for guiding eyes. The initial state is $s_{x}(0)=1$ and we set
$N=300$.

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},$ (35)

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}.$ (36)

In particular, we can demonstrate analytically that the $O(\sqrt{N})$
behaviour holds for the OAT dynamics. To support the claim (36) 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. (36).

  * $\bullet$

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

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

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. (37). 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. (25), 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},$ (38)

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. (38). 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. (38) 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)},$ (39)

where the coefficients $\gamma(A)$ and $\eta(A)$ depend only on the kick
strength $A$. A numerical investigation of the Eq. (39) 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
(39) 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. (39) 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$.

In Ref. [31], the time evolution of the HP approximation to the QFI, $F_{\rm
HP}(t)$, has been related to the local (i.e., finite-time) maximum Lyapunov
exponent $\lambda_{t}$, showing that $F_{\rm HP}(t)\sim e^{2\lambda_{t}t}$.
This relation holds in a transient time window, up to the Ehrenfest time. We
notice that the HP approach of Ref. [31] agrees well with our results. In Fig.
5 we plot the break time of the HP approximation, $t_{\rm HP}$, defined in
analogy to $t_{IB}$ as

$\frac{|F_{Q}(t_{\rm HP})-F_{B}(t_{\rm HP})|}{F_{Q}(t_{\rm HP})}=g_{\rm HP}.$
(40)

Panels (d)-(f) of Fig. 5 show $t_{\rm HP}$ as a function of $N$. Panel (d)
corresponds to the initial state in the classical stable region. Lines in Fig.
5(d) are numerical fit $t_{\rm HP}=O(\sqrt{N})$: the dashed line is
$\sqrt{N}/5.8$ fitting the OAT results, while the solid line is $\sqrt{N}/13$
fitting the QKR results for large $N$. Panel (e) corresponds to the saddle
point, where $t_{\rm HP}=O(\log N)$: the dashed line is $\ln N/2.34$ fitting
the TAT results, while the solid line is $\ln N/2.34$ fitting the QKR results
for large $N$. Panel (f) corresponds to the deep chaotic regime, where $t_{\rm
HP}=O((\log N)^{4})$: the solid line is the fit $(\ln N)^{4}/2020$, the dashed
line is $(\ln N)^{4}/3540$, the dash-dotted line is $(\ln N)^{4}/4070$, and
the dash-double-dotted line is $(\ln N)^{4}/5100$. The scaling with $N$ is
thus analogous to Eqs. (36), (37) and (38), respectively.

Figure 7: Adiabatic dynamics comparison between full quantum model (FQM) and
BMF model. To obtain the population imbalance probability $s_{z}(t)$, we
adiabatically evolve the system from the initial ground states $s_{z}(0)=\pm
1$ according to both quantum Schrödinger equation and BMF equations with
$A=vt$ ($v=10^{-3})$ and $c=1$.

##  V Application in identifying a quantum phase transition

So far, we have shown how the BMF theory can describe the entanglement
dynamics of quantum many-body kicked system up to transient times. Here, we
discuss the application of the BMF approach to identify the critical point of
a second-order quantum phase transition (QPT) that occurs in a bosonic
Josephson junction. This QPT has been recently explored with a BEC with
tuneable interparticle interaction [89, 60] (see Refs. [90, 40, 41, 44] for
studies of entanglement in this system). The system is described by the
Hamiltonian

$\hat{H}=A\hat{J}_{x}-\frac{c}{N}\hat{J}_{z}^{2},$ (41)

that is characterized by a QPT between a gapped symmetric state (for $A/c>1$)
to two degenerate asymmetric states (i.e., $A/c<1$), with critical point (in
the thermodynamic limit) at $A/c=1$. Instead of direct diagonalization, here
we use adiabatic dynamical method to describe and identify the transition
between the two regimes in a system with finite number of particles.
Initially, we assume that the coupling between the two modes is very weak, and
the eigenstates of $\hat{J}_{z}$, i.e., $|\mu=\pm{N}/{2}\rangle$ (or
$s_{z}(0)=\pm 1$) are the ground state of the system, which are double
degenerate, where all atoms occupy in one mode. Then we adiabatically evolve
the system by slowly increasing the coupling strength between the two modes as
$A=vt$. If the coupling strength changes slowly enough, i.e., $v\rightarrow
0$, the quantum adiabatic theorem ensures that the system is always in the
instantaneous ground state. When the coupling exceeds a certain pseudo-phase
transition critical point, the ground state of the system will change from
occupying one of the modes mainly to occupying the two equally. This adiabatic
dynamics shown in Fig. 7 can be used to illustrate the process of phase
transition. Obviously, increasing the number of particles, the consistence of
the results between the full quantum model and the BMF model calculations
becomes better and better. The oscillations of the order parameter $s_{z}$ for
the phase transition after passing the pseudo-critical point will be
suppressed with the decrease of adiabatic evolution velocity $v$.

Figure 8: QFI comparison between full quantum model (FQM) and BMF model. To
obtain $F_{B,Q}$, we adiabatically evolve the system from the initial ground
state $s_{z}(0)=1$ according to both quantum Schrödinger equation and BMF
equations with $A=vt$ ($v=10^{-3})$ and $c=1$.

The corresponding QFI of the instantaneous ground state during adiabatic
evolution is also calculated and the results are shown in Fig. 8. We find that
the QFI has a maximum value at the pseudo-critical point of phase transition.
This can be used as an effective signal to recognize the phase transition of
finite size systems. It is seen that, with the increase of the number of
condensed atoms, the dynamical behaviors of the QFI and its BMF approximation
become more and more similar, and the two quantities become more and more
consistent. This means that for quantum many-body system with large but finite
number of particles, the simple BMF theory can be used to accurately predict
the position of phase transition without solving the complex quantum many-body
model.

##  VI Conclusion

To conclude, we have investigated the entanglement dynamics of a quantum
kicked rotor model by using a beyond mean-field approach. We have found
numerically characteristic scaling laws for the break time $t_{IB}$ when the
quantum Fisher information calculated within the BMF approach departs from the
exact value. The scaling depends on the properties of the initial conditions:
if the initial conditions correspond to a classical regular regime,
$t_{IB}=O(\sqrt{N})$; if the initial conditions correspond to a classical
saddle point, $t_{IB}=O(\log N)$; otherwise, if the initial conditions
correspond to a classical chaotic regime, $t_{IB}=O((\log N)^{4})$. We thus
conclude that unstable saddle point dynamics is the most irregular, in the
sense that quantum effects not captured by the BMF approach emerge more
quickly than in both the regular and the chaotic evolution. Furthermore, in
the chaotic regime, we have shown a relation between $t_{IB}$ and the maximum
Lyapunov exponent, $t_{IB}\sim e^{-\lambda_{L}}$, see Eq. (39). The
application of both the QFI in accurately identifying the critical point of
quantum phase transition in finite-size systems is also discussed by taking
the bosonic Josephson junction as an example.

## Acknowledgement

We acknowledge discussions with Alessio Lerose and Silvia Pappalardi. This
work is supported by the National Natural Science Foundation of China (Grant
No. 11974273) and the Natural Science Fundamental Research Program of Shaanxi
Province of China (Grant No. 2019JM-004).

## References

  * [1] E. P. Wigner, On the statistical distribution of the widths and spacings of nuclear resonance levels, Math. Proc. Cambr. Philos. Soc. 47, 790 (1951).
  * [2] A. Einstein, Zum Quantensatz von Sommerfeld und Epstein, Verh. Dtsch. Phys. Ges. 19, 82 (1917). Reprinted in The Collected Papers of Albert Einstein, A. Engel translator (Princeton University Press, Princeton, 1997).
  * [3] A. A. Bychek, P. S. Muraev, and A. R. Kolovsky, Probing quantum chaos in many-body quantum systems by the induced dissipation, Phys. Rev. A 100, 013610 (2019).
  * [4] S. Chaudhury, A. Smith, B. Anderson, S. Ghose, and P. S. Jessen. Quantum signatures of chaos in a kicked top. Nature 461, 768-771 (2009).
  * [5] D. L. Shepelyansky, Delocalization of quantum chaos by weak nonlinearity, Phys. Rev. Lett. 70, 1787 (1993).
  * [6] S. A. Gardiner, D. Jaksch, R. Dum, J. I. Cirac, and P. Zoller, Nonlinear matter wave dynamics with a chaotic potential, Phys. Rev. A 62, 023612 (2000).
  * [7] Q. Thommen, J. C. Garreau, and V. Zehnle, Classical Chaos with Bose-Einstein Condensates in Tilted Optical Lattices, Phys. Rev. Lett. 91, 210405 (2003).
  * [8] A. D. Martin, C. S. Adams, and S. A. Gardiner, Bright Matter-Wave Soliton Collisions in a Harmonic Trap: Regular and Chaotic Dynamics, Phys. Rev. Lett. 98, 020402 (2007).
  * [9] J. Reslen, C. E. Creffield, and T. S. Monteiro, Dynamical instability in kicked Bose-Einstein condensates, Phys. Rev. A 77, 043621 (2008).
  * [10] J. Martin, B. Georgeot, and D. L. Shepelyansky, Chaotic dynamics of a Bose-Einstein condensate coupled to a qubit, Phys. Rev. E 79, 066205 (2009).
  * [11] Q. Zhang, P. Hanggi, and J. B. Gong, Nonlinear Landau-Zener processes in a periodic driving field, New J. Phys. 10, 073008 (2008).
  * [12] M. Lepers, V. Zehnle, and J. C. Garreau, Tracking Quasiclassical Chaos in Ultracold Boson Gases, Phys. Rev. Lett. 101, 144103 (2008).
  * [13] J. Liu, C. Zhang, M. G. Raizen, and Q. Niu, Transition to instability in a periodically kicked Bose-Einstein condensate on a ring, Phys. Rev. A 73, 013601 (2006).
  * [14] J. Martin, B. Georgeot, and D. L. Shepelyansky, Time Reversal of Bose-Einstein Condensates, Phys. Rev. Lett. 101, 074102 (2008).
  * [15] R. A. Kidd, M. K. Olsen, and J. F. Corney, Quantum chaos in a Bose-Hubbard dimer with modulated tunneling, Phys. Rev. A 100, 013625 (2019).
  * [16] P. Coullet and N. Vandenberghe, Chaotic dynamics of a Bose-Einstein condensate in a double-well trap, J. Phys. B: At. Mol. Opt. Phys. 35, 1593 (2002).
  * [17] J. Cheng, Chaotic dynamics in a periodically driven spin-1 condensate, Phys. Rev. A 81, 023619 (2010).
  * [18] J. Kronjäger, K. Sengstock, and K. Bongs, Chaotic dynamics in spinor Bose-Einstein condensates, New J. Phys. 10, 045028 (2008).
  * [19] J. Cheng, Regular and irregular spin-mixing dynamics in coupled spin-1 atomic and molecular Bose-Einstein condensates, Phys. Rev. A 80, 023608 (2009).
  * [20] F. Haake, M. Kuś, and R. Scharf, Classical and quantum chaos for a kicked top Z. Phys. B Condens. Matter 65, 381 (1987).
  * [21] F. Haake, Quantum Signatures of Chaos (Springer, Berlin, 2010).
  * [22] F. L. Moore, J. C. Robinson, C. F. Bharucha, B. Sundaram, and M. G. Raizen, Atom Optics Realization of the Quantum $\delta$-Kicked Rotor, Phys. Rev. Lett. 75, 4598 (1995).
  * [23] L. Fallani, L. De Sarlo, J. E. Lye, M. Modugno, R. Saers, C. Fort, and M. Inguscio, Observation of Dynamical Instability for a Bose-Einstein Condensate in a Moving 1D Optical Lattice, Phys. Rev. Lett. 93, 140406 (2004).
  * [24] M. Christiani, O. Morsch, N. Malossi, M. Jona-Lasinio, M. Anderlini, E. Courtade, and E. Arimondo, Instabilities of a Bose-Einstein condensate in a periodic potential: an experimental investigation, Opt. Express 12, 4 (2004).
  * [25] G. J. Duffy, S. Parkins, T. Muller, M. Sadgrove, R. Leonhardt, and A. C. Wilson, Experimental investigation of early-time diffusion in the quantum kicked rotor using a Bose-Einstein condensate, Phys. Rev. E 70, 056206 (2004).
  * [26] C. E. Creffield, G. Hur, and T. S. Monteiro, Localization-Delocalization Transition in a System of Quantum Kicked Rotors, Phys. Rev. Lett. 96, 024103 (2006).
  * [27] B. V. Chirikov, A universal instability of many-dimensional oscillator systems, Phys. Rep. 52, 263 (1979).
  * [28] H. J. Stockmann, Quantum Chaos: An Introduction (Cambridge University Press, Cambridge, 1999).
  * [29] R. Lewis-Swan, A. Safavi-Naini, J. Bollinger, and A. Rey, Unifying scrambling, thermalization and entanglement through measurement of fidelity out-of-time-order correlators in the Dicke model, Nat. Comm. 10, 1581 (2019).
  * [30] K. Gietka, J. Chwedeńczuk, T. Wasak, and F. Piazza, Multipartite entanglement dynamics in a regular-to-ergodic transition: Quantum Fisher information approach, Phys. Rev. B 99, 064303 (2019).
  * [31] A. Lerose and S. Pappalardi, Bridging entanglement dynamics and chaos in semiclassical systems, Phys. Rev. A 102, 032404 (2020).
  * [32] L. Pezzè and A. Smerzi, Entanglement, Nonlinear Dynamics, and the Heisenberg Limit, Phys. Rev. Lett. 102, 100401 (2009).
  * [33] P. Hyllus, W. Laskowski, R. Krischek, C. Schwemmer, W. Wieczorek, H. Weinfurter, L. Pezzè, and A. Smerzi, Fisher information and multiparticle entanglement, Phys. Rev. A 85, 022321 (2012).
  * [34] G. Tóth, Multipartite entanglement and high-precision metrology, Phys. Rev. A 85, 022322 (2012).
  * [35] G. Tóth and I. Apellaniz, Quantum metrology from a quantum information science perspective, J. Phys. A 47, 424006 (2014).
  * [36] L. Pezzè, Y. Li, W. Li, and A. Smerzi, Witnessing entanglement without entanglement witness operators, Proc. Natl. Acad. Sci. U.S.A. 113, 11459 (2016).
  * [37] L. Pezzè, A. Smerzi, M. K. Oberthaler, R. Schmied, and P. Treutlein, Quantum metrology with nonclassical states of atomic ensembles, Rev. Mod. Phys. 90, 035005 (2018).
  * [38] C. W. Helstrom, Quantum Detection and Estimation Theory (Academic Press, New York, 1976).
  * [39] S. L. Braunstein, and C. M. Caves, Statistical distance and the geometry of quantum states, Phys. Rev. Lett. 72, 3439 (1994).
  * [40] J. Ma and X. Wang, Fisher information and spin squeezing in the Lipkin-Meshkov-Glick model, Phys. Rev. A 80 012318 (2009).
  * [41] P. Hauke, M. Heyl, L. Tagliacozzo, and P. Zoller, Measuring multipartite entanglement through dynamic susceptibilities. Nat. Phys. 12, 778 (2016).
  * [42] L. Pezzè, M. Gabbrielli, L. Lepori and A. Smerzi, Multipartite entanglement in topological quantum phases, Phys.Rev.Lett. 119 250401 (2017)
  * [43] S. Pappalardi, S. Russomanno, A. Silva, and R. Fazio, Multipartite entanglement after a quantum quench, J. Stat. Mech. 053104 (2017).
  * [44] M. Gabbrielli, A. Smerzi and L. Pezzè, Multipartite Entanglement at Finite Temperature, Sci. Rep. 8, 15663 (2018).
  * [45] M. Gabbrielli, L. Lepori and L. Pezzè, Multipartite-entanglement tomography of a quantum simulator New J. Phys. 21, 033039 (2019).
  * [46] I. Frérot and T. Roscilde, Reconstructing the quantum critical fan of strongly correlated systems using quantum correlations, Nat. Comm 10, 577 (2019).
  * [47] G. Sorelli, M. Gessner, A. Smerzi, and L. Pezzè, Fast and optimal generation of entanglement in bosonic Josephson junctions, Phys. Rev. A 99, 022329 (2019).
  * [48] J. R. Anglin and A. Vardi, Dynamics of a two-mode Bose-Einstein condensate beyond mean-field theory, Phys. Rev. A 64, 013605 (2001).
  * [49] L. J. Fiderer and D. Braun, Quantum metrology with quantum-chaotic sensors, Nat. Commun. 9, 1351 (2018).
  * [50] J. Javanainen, Oscillatory exchange of atoms between traps containing Bose condensates, Phys. Rev. Lett. 57, 3164 (1986).
  * [51] G. J. Milburn, J. Corney, E. M. Wright, and D. F. Walls, Quantum dynamics of an atomic Bose-Einstein condensate in a double-well potential, Phys. Rev. A 55, 4318 (1997).
  * [52] I. Zapata, F. Sols, and A. J. Leggett, Josephson effect between trapped Bose-Einstein condensates, Phys. Rev. A 57, R28 (1998).
  * [53] A. Smerzi, S. Fantoni, S. Giovanazzi, and S. R. Shenoy, Quantum coherent atomic tunneling between two trapped Bose-Einstein condensates, Phys. Rev. Lett. 79, 4950 (1997).
  * [54] J. Javanainen and M. Yu. Ivanov, Splitting a trap containing a Bose-Einstein condensate: Atom number fluctuations, Phys. Rev. A 60, 2351 (1999).
  * [55] M. Albiez, R. Gati, J. Fölling, S. Hunsmann, M. Cristiani, and M. K. Oberthaler, Direct observation of tunneling and nonlinear self-trapping in a single bosonic Josephson junction, Phys. Rev. Lett. 95, 010402 (2005).
  * [56] T. Schumm, S. Hofferberth, L. M. Andersson, S. Wildermuth, S. Groth, I. Bar-Joseph, J. Schmiedmayer, and P. Krüger, Matter-wave interferometry in a double well on an atom chip, Nat. Phys. 1, 57 (2005).
  * [57] L. Pezzè, L. A. Collins, A. Smerzi, G. P. Berman, and A. R. Bishop, Sub-shot-noise phase sensitivity with a Bose-Einstein condensate Mach-Zehnder interferometer, Phys. Rev. A 72, 043612 (2005).
  * [58] D. Ananikian and T. Bergeman, Gross-Pitaevskii equation for Bose particles in a double-well potential: Two-mode models and beyond, Phys. Rev. A 73, 013604 (2006).
  * [59] R. Gati and M. K. Oberthaler, A bosonic Josephson junction, J. Phys. B: At., Mol. Opt. Phys. 40, R61 (2007).
  * [60] A. Trenkwalder, et al., Quantum phase transitions with parity-symmetry breaking and hysteresis, Nat. Phys. 12, 826 (2016).
  * [61] G. Spagnolli, et al., Crossing over from attractive to repulsive interactions in a tunneling bosonic josephson junction, Phys. Rev. Lett. 118, 230403 (2017).
  * [62] J. I. Cirac, M. Lewenstein, K. Mølmer, and P. Zoller, Quantum superposition states of Bose-Einstein condensates, Phys. Rev. A 57, 1208 (1998).
  * [63] M. J. Steel and M. J. Collett, Quantum state of two trapped Bose-Einstein condensates with a Josephson coupling, Phys. Rev. A 57, 2920 (1998).
  * [64] T. Zibold, E. Nicklas, C. Gross, and M. K. Oberthaler, Classical bifurcation at the transition from Rabi to Josephson dynamics. Phys. Rev. Lett. 105, 204101 (2010).
  * [65] C. Gross, T. Zibold, E. Nicklas, J. Estève, and M. K. Oberthaler, Nonlinear atom interferometer surpasses classical precision limit, Nature 464, 1165 (2010).
  * [66] F. T. Arecchi, E. Courtens, R. Gilmore and H. Thomas, Atomic Coherent States in Quantum Optics, Phys. Rep. 6, 2211 (1972).
  * [67] G. Ferrini, D. Spehner, A. Minguzzi, and F. W. J. Hekking, Effect of phase noise on quantum correlations in Bose-Josephson junctions, Phys. Rev. A 84, 043628 (2011).
  * [68] M. Kitagawa and M. Ueda, Squeezed spin states, Phys. Rev. A 47, 5138 (1993).
  * [69] A. S. Sørensen, L.-M. Duan, J. I. Cirac, and P. Zoller, Many-particle entanglement with Bose-Einstein condensates, Nature 409, 63 (2001).
  * [70] M. F. Riedel, P. Böhi, Y. Li, T. W. Hänsch, A. Sinatra, and P. Treutlein, Atom-chip-based generation of entanglement for quantum metrology, Nature 464, 1170 (2010).
  * [71] R. Schmied, J.-D. Bancal, B. Allard, M. Fadel, V. Scarani, P. Treutlein, and N. Sangouard, Bell correlations in a Bose-Einstein condensate, Science 352, 441 (2016).
  * [72] M. Takeuchi, S. Ichihara, T. Takano, M. Kumakura, T. Yabuzaki, and Y. Takahashi, Spin Squeezing Via One-Axis Twisting with Coherent Light, Phys. Rev. Lett. 94, 023003 (2005).
  * [73] I. D. Leroux, M. H. Schleier-Smith, and V. Vuletić, Implementation of Cavity Squeezing of a Collective Atomic Spin, Phys. Rev. Lett. 104, 073602 (2010).
  * [74] M. H. Schleier-Smith, I. D. Leroux, and V. Vuletić, Squeezing the collective spin of a dilute atomic ensemble by cavity feedback, Phys. Rev. A 81, 021804 (2010).
  * [75] J. G. Bohnet, K. C. Cox, M. A. Norcia, J. M. Weiner, Z. Chen, and J. K. Thompson, Reduced spin measurement back-action for a phase sensitivity ten times beyond the standard quantum limit, Nat. Photonics 8, 731 (2014).
  * [76] O. Hosten, N. J. Engelsen, R. Krishnakumar, and M. A. Kasevich, Measurement noise 100 times lower than the quantum projection limit using entangled atoms, Nature 529, 505 (2016).
  * [77] A. Micheli, D. Jaksch, J. I. Cirac, and P. Zoller, Many-particle entanglement in two-component Bose-Einstein condensates, Phys. Rev. A 67, 013607 (2003).
  * [78] H. Strobel, W. Muessel, D. Linnemann, T. Zibold, D. B. Hume, L. Pezzè, A. Smerzi, and M. K. Oberthaler, Fisher information and entanglement of non-Gaussian spin states, Science 345, 424 (2014).
  * [79] W. Muessel, H. Strobel, D. Linnemann, T. Zibold, B. Juliá-Diaz, and M. K. Oberthaler, Twist-and-turn spin squeezing in Bose-Einstein condensates, Phys. Rev. A 92, 023603 (2015).
  * [80] B. Julià-Díaz, E. Torrontegui, J. Martorell, J. G. Muga, and A. Polls, Fast generation of spin-squeezed states in bosonic Josephson junctions, Phys. Rev. A 86, 063623 (2012).
  * [81] B. Julià-Díaz, T. Zibold, M.K. Oberthaler, M. Mele-Messeguer, J. Martorell, and A. Polls, Dynamic generation of spin- squeezed states in bosonic Josephson junctions, Phys. Rev. A 86, 023615 (2012).
  * [82] V. S. Shchesnovich and M. Trippenbach, Fock-space WKB method for the boson Josephson model describing a Bose-Einstein condensate trapped in a double-well potential, Phys. Rev. A 78, 023611 (2008).
  * [83] A. Lerose, J. Marino, B. Zunkovic, A. Gambassi, and A. Silva, Chaotic Dynamical Ferromagnetic Phase Induced by Nonequilibrium Quantum Fluctuations, Phys. Rev. Lett. 120, 130603 (2018); A. Lerose, B. Zunkovic, J. Marino, A. Gambassi, and A. Silva, Impact of nonequilibrium fluctuations on prethermal dynamical phase transitions in long-range interacting spin chains, Phys. Rev. B 99, 045128 (2019).
  * [84] The QFI is given by the maximum between Eq. (9) and  $\frac{F_{Q}(t)}{N}=N\left(1-\cos^{2N-2}\frac{ct}{N}\right)-(N-1)\frac{\alpha_{Q}(t)}{2}.$ (42) Equation (42) is larger than Eq. (9) for $ct\simeq N\pi/2-2\sqrt{N}$ and even values of $N$.
  * [85] A. Lerose and S. Pappalardi, Origin of the slow growth of entanglement entropy in long-range interacting spin systems, Phys. Rev. Research 2, 012041(R) (2020).
  * [86] A. Lerose and S. Pappalardi, private communication. S. Pappalardi, Entanglement dynamics and chaos in long-range quantum systems PhD thesis (2020). Available at https://iris.sissa.it.
  * [87] C. K. Law, H. T. Ng, and P. T. Leung, Coherent control of spin squeezing, Phys. Rev. A 63, 055601 (2001).
  * [88] J. Liu, W. Wang, C. Zhang, Q. Niu, and B. Li, Fidelity for the quantum evolution of a Bose-Einstein condensate, Phys. Rev. A 72, 063623 (2005).
  * [89] V. V. Ulyanov and O. B. Zaslavskii, New methods in the theory of quantum spin systems, Phys. Rep. 216, 179 (1992).
  * [90] J. Vidal, G. Palacios and R. Mosseri, Entanglement in a second-order quantum phase transition, Phys. Rev. A 69, 022107 (2004).