The failure of semiclassical approach in the dissipative fully-connected Ising model

Dissipative spin models are currently attracting wide interest, because they describe genuine nonequilibrium states of matter and at the same time, can be realized in various platforms from superconducting circuits to Rydberg atoms [1, 2, 3, 4, 5, 6, 7, 8, 9]. The model consists of an ensemble of spins that are subject to dissipation caused by external baths. The dynamics is described by the Lindblad equation, obtained by integrating out the baths’ degrees of freedom [10].

Great efforts have been taken in the investigation of various spin models. In the central spin model [11], the dissipative phase transition was located according to the closing of the Liouvillian gap. The transverse-field Ising model was thoroughly studied both under the mean-field approximation and beyond [12, 13, 14, 15, 16, 17, 17], in which the bistability of steady states in some region of the parameter space was found to be replaced by a first-order phase transition after the spatial correlation was correctly considered. As the couplings between spins are all-to-all and then the system can be seen as a huge spin, the semiclassical approach was adopted. The magnetization was found to display an everlasting oscillation in the thermodynamic limit, indicating that the time translational symmetry is spontaneously broken into a discrete one [18, 19, 20]. If the dissipation acts in the eigenbasis of the transverse field, then a continuous dissipative phase transition manifests itself as continuous order parameters with discontinuous derivatives [21, 22, 23, 24]. The XYZ-Heisenberg model was also studied by different approximation schemes [25, 26, 27, 28, 29] to clarify its phase diagram.

These studies focused on the time-independent Liouvillians. But much less is known as the Liouvillian changes periodically with time [30]. The topic of this work is to study the response to a time-periodic Liouvillian. On the other hand, in closed quantum systems, the response to a time-periodic Hamiltonian has been under intensive investigations. The kicked top models were studied both theoretically [31, 32, 33, 34, 35, 36, 37, 38] and experimentally[39, 40, 41], which is known to exhibit a transition between a regular dynamical phase and a chaotic one, depending on the value of the kicking strength. Similar chaotic behavior was found in the Lipkin-Meshkov-Glick model as the parameters change periodically with time [42, 43, 44, 45]. The study in this paper can be seen as an investigation of the dissipation effect on the chaotic dynamics in the periodically-driven spin models.

As a concise example, we study the fully-connected Ising model in the presence of collective dissipation. Without dissipation, chaotic behavior appears in the presence of a strongly oscillating external field [43, 44, 45]. We find that, the chaotic behavior is robust against weak dissipation, if the semiclassical approximation is taken. As the oscillating amplitude of field increases, a periodic response changes into a subharmonic oscillation, and then into a chaotic behavior. But the numerical simulation shows that, beyond the semiclassical approximation, only the periodic response can survive the quantum fluctuation, but neither the subharmonic nor the chaotic dynamics can be observed. The semiclassical approximation works only as the oscillating amplitude is small in the periodic-response regime, but it fails as the amplitude is large.

The paper is organized as follows. We introduce the model in Sec. II. Section III contributes to the discussion of the semiclassical results. The exact numerical simulations of the dynamics of magnetizations are discussed in Sec. IV. The Floquet Liouvillian spectrum is studied in Sec. V. Finally, Sec. VI summarizes our results.

We consider the transverse field Ising model with all-to-all couplings, and a sinusoidal modulation added to the external field. The Hamiltonian is written as

where $N$ denotes the total number of spins, $g$ denotes the coupling strength, and $\hat{J}_{\alpha}=\sum_{j}\hat{\sigma}_{j}^{\alpha}/N$ denote the collective spin operators with $\hat{\sigma}_{j}^{\alpha}$ being the Pauli matrices of the $j$th spin and $\alpha=x,y,z$. $\Gamma(t)$ denotes the time-dependent external field, which is supposed to be $\Gamma(t)=\Gamma_{0}+A\ sin(\omega_{0}t)$, where $A$, $\Gamma_{0}$ and $\omega_{0}$ are the oscillating amplitude, mean value and frequency, respectively. We set $g=1$ as the unit of energy throughout the paper.

In the presence of dissipation, the dynamics of the system is described by the Lindblad equation [46]. The density matrix satisfies

where $\hat{J}_{\pm}=\left(\hat{J}_{x}\pm i\hat{J}_{y}\right)/2$ are the jump operators, and $\gamma_{c}$ is the dissipation rate. As $A=0$ and then $\Gamma(t)=\Gamma_{0}$ is a constant, the solution of Eq. (2) was studied previously [24]. The jump operator forces the spins to be aligned in the negative $z$-direction, while the interaction between spins favors an alignment in the $x$-direction. Their interplay results in a steady state, which is either ferromagnetic (with nonzero magnetization in the $x$-direction) or paramagnetic. Here we extend to the case of $A\neq 0$, in which the field oscillation prevents a steady state being reached, and then we expect a nontrivial dynamical behavior.

The semiclassical approach is frequently employed for solving the fully-connected models. We choose the order parameters to be $m_{\alpha}=\left\langle\hat{J}_{\alpha}\right\rangle=Tr[\hat{\rho}\hat{J}_{\alpha}]$. By ignoring the correlations (i.e., setting $\left\langle\hat{J}_{\alpha}\hat{J}_{\beta}\right\rangle=\left\langle\hat{J}_{\alpha}\right\rangle\left\langle\hat{J}_{\beta}\right\rangle$) in the limit $N\to\infty$, we obtain a nonlinear system of differential equations, which read

It is easy to see that $\left|{\bf{m}}(t)\right|=\sqrt{m_{x}^{2}+m_{y}^{2}+m_{z}^{2}}$ is a constant of motion, which can be set to unity without loss of generality. ${\bf{m}}$ is moving on a Bloch sphere, and the initial state can be described by the azimuthal angles $\left(\theta,\phi\right)$, which are defined by $m_{x}=\sin\theta\cos\phi$, $m_{y}=\sin\theta\sin\phi$ and $m_{z}=\cos\theta$.

Since the coefficient $\Gamma(t)$ is a periodic function of $t$ with the period $2\pi/\omega_{0}$, one may naively think that the solution ${\bf{m}}(t)$ is also a periodic function. This is the case for small $A$, but not true for large $A$. We note that Eq. (3) bears some resemblance to the Lorenz equation [47], one famous example of the deterministic chaos in the classical systems. In Eq. (3), the possibility of chaos comes from the fact that $\Gamma(t)$ is time-periodic, even the trajectory is limited on a two-dimensional sphere. Indeed, we observed a periodic ${\bf{m}}(t)$ for small $A$, but a chaotic ${\bf{m}}(t)$ for large $A$.

Next we choose $\Gamma_{0}=1$ and $\omega=1$ as an example for the demonstration. In Fig. 1(a), we display the evolution of $m_{x}$ for small $A$ ($A=0.1$) and 25 initial states that are evenly distributed on the Bloch sphere. Except for the initial state at the south pole (${\bf{m}}=(0,0,1)$), all the others eventually evolve into one of the two oscillation modes that are symmetric to each other. The two oscillation modes have the same period, which is exactly the driving period ($2\pi/\omega_{0}$). For small $A$, the long-time response is insensitive to a small deviation in the initial state. On the other hand, Fig. 1 shows the evolution of $m_{x}$ for a large $A$ ($A=1.0$). Till the largest evolution time that is accessible, no periodicity is observed. More importantly, the long-time response is extremely sensitive to the initial condition. Even for a tiny deviation in the initial state ($\theta_{1}-\theta_{0}=10^{-7}$), $m_{x}(t)$ displays a significant difference as $t$ is as large as a few hundreds (see the lines of different colors in Fig. 1).

For the values of $A$ between $0.1$ and $1.0$, ${{m}_{x}}(t)$ displays abundant dynamical behaviors. As $A$ increases up to a certain critical value, the time period is doubled. In Fig. 2, we plot ${{m}_{x}}(t)$ for $A=0.73$. It is easy to see that the period is not $2\pi/\omega_{0}$, instead, it becomes $4\pi/\omega_{0}$. As $A$ increases further, the period-doubling happens again. For example, for $A=0.7345$, the period becomes $8\pi/\omega_{0}$ (see Fig. 2). The period-doubling bifurcation is well-known in the classical nonlinear systems. Usually, a cascade of period-doubling bifurcations lead to chaos [48]. This explains why we observe a chaotic dynamics as $A$ is as large as $A=1$.

Depending on the values of $A$, the system exhibits the periodic, subharmonic or chaotic responses. We plot the trajectory of the vector ${\bf{m}}(t)$ on the Bloch sphere in Appendix A, which provides more evidence for the existence of different dynamical behaviors.

In our model, the subharmonic response (i.e., doubled period) to a time-periodic Liouvillian must be distinguished from that in the Floquet time-crystals. As will be shown next, the subharmonic response in our model can only be observed in the semiclassical limit. It cannot survive the quantum fluctuation that is unavoidable at finite $N$.

To locate the chaotic phase in the parameter space, we calculate the Lyapunov exponent. The defining property of a chaotic system is the extreme sensitivity of trajectories to the initial condition. Two points that are initially close will drift apart exponentially over time. The Lyapunov exponent [49] provides a quantitative measure for this, which is defined as the average exponential rate of convergence or divergence between adjacent trajectories in the phase space. Especially, the largest Lyapunov exponent (LLE) is frequently employed for determining whether a nonlinear system is chaotic. If the LLE is greater than zero, two initial points will depart from each other exponentially, and then the system is chaotic [50], otherwise, it is not. Figure 3 displays the dependence of the LLE on $A$ and $\gamma_{c}$. The area in red or yellow has a positive LLE, while the area in blue or green has a negative LLE. The chaotic phase is clearly distinguishable in the parameter space. An interesting finding is that for a fixed $\gamma_{c}$ in the interval $\left(0,1\right)$, the dynamics is chaotic only for an intermediate amplitude of oscillating field, but the dynamics is regular either if $A$ is too large or too small. In the presence of strong dissipation ($\gamma_{c}>1.5$), there is no chaotic dynamics for whatever $A$.

To check the validity of semiclassical approximation, we numerically simulate the real-time dynamics at finite $N$, by exploiting the permutation symmetry of fully-connected models. The Dicke basis with maximum angular momentum is defined as [51]

where $C_{N}^{\frac{N}{2}+M}$ is the binomial coefficient, and $\sigma_{j}=\pm 1/2$ represents the spin-up and down states, respectively. And $M=-N/2,-N/2+1,\cdots,N/2$ is the magnetization in the $z$-direction. The initial state is supposed to be a pure state with all the spins aligned in the same direction. We use the azimuthal angles $(\theta,\phi)$ to indicate the initial direction, then the initial state can be expressed in the Dicke basis as

Equation (2) has the permutation symmetry, its solution can then be expressed as $\hat{\rho}(t)=\sum_{M,M^{\prime}}\rho_{M,M^{\prime}}(t)\ket{M}\bra{M^{\prime}}$, where $\rho_{M,M^{\prime}}(t)$ is the density matrix in the Dicke basis. Now Eq. (2) changes into a system of ordinary differential equations for $\rho_{M,M^{\prime}}$, which are solved numerically. The permutation symmetry reduces the dimension of Hilbert space from $2^{N}$ to $N$, and then allows us to access a large system with the number of spins $N\sim 10^{2}-10^{3}$.

In Fig. 4, we compare the dynamics of magnetizations at finite $N$s and in the semiclassical limit ($N=\infty$), as $A=0.1$ is in the semiclassical periodic regime. At finite $N$s, the magnetizations display periodic oscillations in the asymptotic long time. The curve of $m_{z}(t)$ at $N=50$ is already very close to the semiclassical one. As $N$ increases, $m_{z}(t)$ goes even closer to the semiclassical result (see Fig. 4 the inset). As $N\to\infty$, we expect that the numerical results should coincide with the semiclassical results. If $A$ is small, then the semiclassical approximation is good for large enough $N$s.

Figure 4 shows the comparison as $A=1.0$ is in the semiclassical chaotic regime. The results are significantly different from those at $A=0.1$. Up to $N=200$, we find no similarity between the numerical result and the semiclassical one. For $N$ ranging between $50$ and $200$, the initial transient behavior of $m_{z}(t)$ always quickly evolves into the asymptotic behavior - periodic oscillation. And the $m_{z}(t)$s for $N=50$ and $N=200$ have the same period, with only their amplitude being different. But in the semiclassical result, $m_{z}(t)$ is aperiodic at arbitrarily long time. These observations suggest that the exact numerical solutions do not conserve to the semiclassical one in the limit $N\to\infty$. If $A$ is large, then the semiclassical approximation is bad, even for large $N$s. This seems to be strange for a fully-connected model. We will further discuss this discrepancy in next section.

We also compare the dynamics of magnetizations for different $A$s, as $N$ is fixed. Figure 5 and 5 display $m_{z}(t)$s and $m_{x}(t)$s, respectively, for the values of $A$ ranging from $0.1$ to $1.0$. The magnetization in the $z$-direction always displays a periodic oscillation. And the oscillations for different $A$s are in phase (see Fig. 5 the inset). On the other hand, the magnetization in the $x$-direction displays two different dynamical modes, depending on the value of $A$. As seen in Fig. 5, as $A$ is small ($A=0.1,0.3$), $m_{x}$ displays an everlasting oscillation with the period $2\pi/\omega_{0}$. But as $A$ is large ($A=0.73,1.0$), $m_{x}$ rapidly decays to zero. For an intermediate $A$ ($A=0.5$), $m_{x}$ maintains an oscillation for a relatively long time, but the decay can still be clearly seen. The decay of $m_{x}$ is a signal of the discrepancy between the exact results and the semiclassical approximations, for the latter have non-decaying magnetizations in all directions (see App. A).

In the closed systems, it is widely believed that the semiclassical approximation becomes exact for the fully-connected models if the number of spins goes to infinity. In the open systems, a similar result was obtained [52], as the Liouvillian is time-independent. But in this paper, our numerical results suggest that the semiclassical approximation is bad even in the limit $N\to\infty$, if there is a big time-periodic term in the Liouvillian. Since we can only obtain the exact results at finite $N$, a scaling analysis is helpful for confirming our viewpoint. Next we give a scaling analysis of the Floquet Liouvillian spectrum.

The Lindblad equation can be expressed in a vectorized form as $d\hat{\rho}/dt=\hat{\hat{\mathcal{L}}}\left(\hat{\rho}\right)$, where $\hat{\hat{\mathcal{L}}}$ is the so-called Liouvillian superoperator (or Liouvillian in short), which is a non-Hermitian linear operator acting on the vector space of density matrices. For the dissipative systems with time-independent $\hat{\hat{\mathcal{L}}}$, it is well known that the eigenvalues and eigenvectors of $\hat{\hat{\mathcal{L}}}$ determine the dynamics. The eigenvalues of $\hat{\hat{\mathcal{L}}}$ (Liouvillian spectrum) are complex numbers. More important, if the system size $N$ is finite, all the eigenvalues must have negative real part, except one that is zero.

These notations can be generalized to the case of time-periodic $\hat{\hat{\mathcal{L}}}$. For a Lindblad equation with $\hat{\hat{\mathcal{L}}}(t)=\hat{\hat{\mathcal{L}}}(t+T)$, there exist a complete set of solutions written as $\hat{\rho}(t)=e^{\lambda t}\hat{\varrho}(t)$, where $\hat{\varrho}(t)=\hat{\varrho}(t+T)$ is the time-periodic part of density matrix, according to the Floquet theorem. And $\hat{\varrho}(t)$ satisfies

Then $\lambda$ can be seen as the eigenvalue of $\left(\hat{\hat{\mathcal{L}}}(t)-d/dt\right)$, which is an operator acting on the generalized vector space of density matrices, just as $\left(\hat{H}-id/dt\right)$ is an operator acting on the generalized Hilbert space (Sambe space) for a closed system with time-periodic Hamiltonian. While the eigenstates of $\left(\hat{H}-id/dt\right)$ are called the Floquet spectrum, the eigenstates of $\left(\hat{\hat{\mathcal{L}}}(t)-d/dt\right)$ are called the Floquet Liouvillian spectrum.

Compared to the Floquet spectrum or Liouvillian spectrum, much less is known about the Floquet Liouvillian spectrum. For the dissipative systems, we guess that the Floquet Liouvillian spectrum has the same properties as the Liouvillian spectrum, i.e., all the complex eigenvalues have negative real parts except for a unique one that is zero. The numerics support our guess. Figure 6 displays the Floquet Liouvillian spectrum at $N=30$. For both $A=0.1$ and $A=1.0$, we see that the rightmost eigenvalue in the complex plane is zero, which is non-degenerate. And the others are to the left of zero.

The Floquet Liouvillian spectrum completely determines whether the dynamics is periodic, subharmonic or chaotic. To see it, we arrange all the eigenvalues as

with the corresponding eigenvectors being $\hat{\varrho}_{0}(t)$, $\hat{\varrho}_{1}(t)$, $\hat{\varrho}_{2}(t),\cdots$. For an arbitrary initial state, the solution of Lindblad equation can be formally expressed as

where the coefficients $K_{j}$ depend on the initial state. In the asymptotic long time, the terms with $\text{Re}\lambda_{j}<0$ all decay to zero, and $1/\left|\text{Re}\lambda_{j}\right|$ is just the decay time of the $j$-th mode. At finite $N$, all the $\lambda_{j}$ with $j>0$ have negative real parts, therefore, the density matrix in the asymptotic long time becomes $\hat{\rho}(t)=K_{0}\hat{\varrho}_{0}(t)$, which is exactly periodic with the period $T$. We then expect that the asymptotic behavior is always periodic at finite $N$. The situation is more complicated in the thermodynamic limit. As found in the previous studies of time-independent Liouvillians, there exist possibilities that the real parts of some $\lambda_{j}$ decrease with increasing $N$ and vanish in the limit $N\to\infty$. As a consequence, the asymptotic density matrix becomes $\hat{\rho}(t)=\sum_{j=0}^{L}K_{j}e^{i\text{Im}\left(\lambda_{j}\right)t}\hat{\varrho}_{j}(t)$, where $L$ is the number of eigenvalues with vanishing real parts. Additionally, if these $\lambda_{j}$ have also nonvanishing imaginary parts, then $\hat{\rho}(t)$ possibly display subharmonic or chaotic behaviors, depending on the values of $\text{Im}\left(\lambda_{j}\right)$. Conversely, if the real parts of $\lambda_{j}$ with $j>0$ are all finite in the limit $N\to\infty$, then subharmonic oscillation or chaotic behavior are both impossible.

According to the above argument, we perform a scaling analysis of the Floquet Liouvillian gap, which is defined as $\Delta=\lambda_{0}-\text{Re}\lambda_{1}$. Figure 7 plots $\Delta$ as a function of $1/N$. The dots are the numerical results, while the lines are the fitted curves. The gap at $A=1.0$ is significantly smaller than the gap at $A=0.1$. But in both cases, we clearly see that the gap does not go to zero as $N\to\infty$. This indicates that no $\lambda_{j}$ with $j>0$ has vanishing real parts, therefore, the asymptotic behaviors are periodic for both $A=0.1$ and $A=1.0$, even in the limit $N\to\infty$. Such a result is consistent with our previous simulation of the real-time dynamics of magnetizations.

In summary, we study the fully-connected Ising model with a time-periodic external field and subject to a dissipation, by using both the semiclassical approach and the exact numerical simulation. If the field amplitude is small, both the semiclassical approach and the numerical simulation predict a perfect periodic oscillation of magnetizations. And the numerical results in the thermodynamic limit are consistent with the semiclassical one.

As the field amplitude increases, the semiclassical approach predicts the period-doublings or subharmonic oscillations, and a series of period-doublings finally lead to the chaotic dynamics of magnetizations, which is confirmed by the calculations of Lyapunov exponents. On the contrary, the numerical simulation show that the magnetizations are always oscillating periodically, whatever the field amplitude is. No subharmonic or chaotic dynamics are observed, even we choose the number of spins to be as large as a few hundreds in the simulation. We then analyze the Floquet Liouvillian gap, which conserves to a finite value in the thermodynamic limit, for either small or large field amplitude. We argue that a finite gap is another evidence of the periodic oscillations of observables.

For the fully-connected models, the semiclassical approximation is generally believed to be good for sufficiently large number of spins. But we find that this is not the case if the time-periodic field and the dissipation are both present. For a large field amplitude, the predictions from the semiclassical approach and the numerical method are qualitatively different. Our finding suggests that one should be more careful when using the semiclassical approach in the case of time-periodic Liouvillians.

This paper is supported by NSFC under Grant Nos. 11774315, 11835011, and 12174346, and by the Junior Associates program of the Abdus Salam International Center for Theoretical Physics.