Quantum quench dynamics of the Jaynes-Cummings-Hubbard model with weak nearest-neighbor hopping

Every term in Hamiltonian (1) is quadratic and can be usually solved by Fourier transform, especially all cavities are independent as no photon tunneling among them. Meanwhile, the spin operator can be replaced by the fermionic operator because they share the same standard fermionic anticommutation relations. But the breakdown of this simple replacement occurs as the photon tunneling strength increases owing to the indirect coupling between atoms induced by the nonlocality of the photon, namely, on the different sites spin operators obey commutation relations other than anticommutation relations of the fermion. Strictly speaking, this nonlocality can be solved using the Jordan-Wigner transform. Unfortunately, the linear terms in the Hamiltonian (1) bring indelible Jordan-Wigner factors, then the Fourier transform does not work well. Thus, we have to consider the case in which the coupling strength between cavities is weak enough to introduce fermionic approximation securely.

The Hamiltonian (1) can be rewritten as follows in the fermionic approximation [40]:

$$\begin{eqnarray} H(t)&= \sum^{M}_{j}(\omega_ca_{j}^{\dagger}a_{j}-J(a_{j}^{\dagger}a_{j+1}+a_{j+1}^{\dagger}a_{j}) +\epsilon b^{\dagger}_{j}b_{j}{}\nonumber\\ & \quad+g(t)(a_{j}b^{\dagger}_{j}+a^{\dagger}_{j}b_{j})). \end{eqnarray} \tag{ 2 }$$

Here the spin operators $\sigma^{+}(\sigma^{-})$ are replaced by fermionic operators $b^{\dagger}(b)$. If the periodic boundary condition is imposed, the Hamiltonian will be invariant under the transform $(a_{j}^{\dagger}, a_{j})\longrightarrow (a_{j+1}^{\dagger}, a_{j+1}$) and $(b^{\dagger}_{j}, b_{j})\longrightarrow (b^{\dagger}_{j+1}, b_{j+1})$, which indicates that the total quasimomentum $q=q_a+q_b$ of the system is conserved. Within this approximation, the total quasimomentum of the cavity field is $q_a\equiv\sum^{M-1}_{k=0}(\frac{2\pi k}{M})a_{k}^{\dagger}a_{k}$ (mod $2\pi$) and the total quasimomentum of the atom is $q_b\equiv\sum^{M-1}_{k=0}(\frac{2\pi k}{M})b_{k}^{\dagger}b_{k}$ (mod $2\pi$). $a_{k}^{\dagger}=\frac{1}{\sqrt M}\sum^{M}_{j=1}e^{i(j2\pi k/M)}a_{j}^{\dagger}$ and $b_{k}^{\dagger}=\frac{1}{\sqrt M}\sum^{M}_{j=1}e^{i(j2\pi k/M)}b^{+}_{j}$ are the creation operators in the k-th Bloch state.

We first restrict to the Hilbert space $\mathcal{H}$ with the total number of excitations N ($N=\sum_{k}(a_{k}^{\dagger}a_{k}+b_{k}^{\dagger}b_{k}$). The dimension of this space is defined as

$$\begin{equation} D=\sum^{N}_{s=0}M\frac{(N+M-s-1)!}{(N-s)!(M-s)!s!}, \end{equation} \tag{ 3 }$$

which grows with the system size. Here, the dimension of the Hilbert space of the JCH model grows faster with the increasing system size than the Bose-Hubbard model [41]. For example, $D_B=126$ and $D_J=1002$ for $N=M=5$, $D_B=462$ and $D_J=5336$ for $N=M=6$, $D_B=1716$ and $D_J=28814$ for $N=M=7$, and $D_B=6435$ and $D_J=157184$ for $N=M=8$. D_B and D_J represent the Hilbert space dimensions of the Bose-Hubbard model and the JCH model, respectively.

We decompose the total Hilbert space $\mathcal{H}$ into M subspaces according to the total quasimomentum q of the system, i.e., $\mathcal{H}=\bigoplus^{M-1}_{q=0}\mathcal{H}^{q}$. At this point, the density matrix is always block diagonal with respect to the q subspaces, i.e., $\rho(t)=\bigoplus^{M-1}_{q=0}\rho^{q}(t)$. Because the system behaves quantitatively similarly in all the q subspaces, so one can focus on a specific q subspace. In this paper, we assume that q has a specific value of 1 and take the normalization $\textrm{Tr}[\rho^{q}(t)] = 1$ [41]. Note that we will remove the subscript q in the following discussion for simplicity.

Suppose that the system starts to be in a thermal equilibrium state with temperature T_i and the parameter g_i . The corresponding Hamiltonian is denoted as H_i , and the n-th eigenvalue and eigenstate of H_i are denoted as $E^{i}_n$ and $|\Psi^{i}_n\rangle$, respectively. The initial density matrix of the whole system is then in the representation of $|\Psi^{i}_n\rangle$ [41]

$$\begin{equation} \rho_i=\sum^{D_q}_{n=1}\frac{1}{Z_i}e^{-\beta_iE^{i}_n} |\Psi^{i}_n\rangle\langle \Psi^{i}_n|, \end{equation} \tag{ 4 }$$

where $Z_i=\sum^{D_q}_{n=1}e^{-\beta_iE^{i}_n}$ is the partition function and $D_q\approx D/M$ is the dimension of the specific q subspace.

Assume that the system is quenched by changing the coupling strength from g_i to g_f at time $t = 0$. The Hamiltonian will be denoted as H_f , and the associated eigenvalues and eigenstates will be denoted as $E^{f}_l$ and $|\Psi^{f}_l\rangle$ when $t\geq0$. The density matrix at an arbitrary time is

$$\begin{eqnarray} \rho(t)&= U(t)\rho_iU(t)^{\dagger}\\ &= \displaystyle\sum^{D_q}_{n,l=1}e^{-i(E^{f}_n-E^{f}_l)t} \langle \Psi^{f}_n|\rho_i|\Psi^{f}_l\rangle|\Psi^{f}_n\rangle\langle \Psi^{f}_l|, \end{eqnarray} \tag{ 5 }$$

where $U(t)= \mathcal{T}e^{-i\int^t_0\textrm{d}\tau H(\tau)}$, $\mathcal{T}$ means time ordering. The time-averaged density matrix is defined as

$$\begin{equation} \bar{\rho}=\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}\textrm{d}t\rho(t)=\sum^{D_q}_{n,l=1;E^{f}_n=E^{f}_l} \langle \Psi^{f}_n|\rho_i|\Psi^{f}_l\rangle|\Psi^{f}_n\rangle\langle \Psi^{f}_l|. \end{equation} \tag{ 6 }$$

It can be seen that the time-averaged density matrix keeps a memory of the initial state after the quantum quench. Generally, the spectrum is nondegenerate in many models [6–8,42,43]. However, in the JCH model, the spectrum has exact degeneracies in the final Hamiltonian H_f . Therefore, the time-averaged density matrix is

$$\begin{eqnarray} \bar{\rho}&= \sum_{l}p_l|\Psi^{f}_l\rangle\langle \Psi^{f}_l|+\sum_{d}\sum_{q^{f}_d}\sum_{q\prime^{f}_d}p_d|q^{f}_d\rangle\langle q\prime^{f}_d|, \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} p_l&= \langle \Psi^{f}_l|\rho_i|\Psi^{f}_l\rangle=\frac{1}{Z_i}\sum^{D_q}_{n=1}e^{-\beta_iE^{i}_n}{| \langle \Psi^{i}_n|\Psi^{f}_l\rangle|}^2, \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} p_d&= \langle q^{f}_d|\rho_i|q\prime^{f}_d\rangle=\frac{1}{Z_i}\sum^{D_q}_{n=1}e^{-\beta_iE^{i}_n}\langle q^{f}_d|\Psi_n^{i}\rangle\langle \Psi_n^{i}|q\prime^{f}_d\rangle,\qquad \end{eqnarray} \tag{ 9 }$$

where l labels nondegenerate eigenstates of H_f and d labels the basis of the degenerate subspaces. In this case, the observable O has off-diagonal matrix elements in the H_f eigenstate basis. So, we must compute all the overlaps $\langle q^{f}_d|\Psi^{i}_n\rangle$ and $\langle q\prime^{f}_d|O|q^{f}_d\rangle$ to get time-averaged results for an observable O, where $|q^{f}_d\rangle$ and $|q\prime^{f}_d\rangle$ are degenerate states of the degenerate subspaces after the quantum quench [44].

When the system has a fixed particle number in thermal equilibrium with a heat bath at some temperature T_c , which can be described using a canonical ensemble (subscript c),

$$\begin{equation} \rho_c=\frac{e^{-\beta_fH_f}}{tr(e^{-\beta_fH_f})}, \end{equation} \tag{ 10 }$$

where $\beta_f=1/T_c$ is the final inverse temperature after the quantum quench.

In fig. 1, we focus on the comparison between the time-averaged density matrix $\bar{\rho}$ and the canonical ensemble $\rho_c$ as a function of the eigenvalues $E_l^f$ for a wide range of coupling strengths between the cavity field and the atom. The red dots represent the occupation vs. the eigenvalues $E_l^f$ by $\bar{\rho}$, while the green dots are obtained by $\rho_c$. As can be seen from fig. 1(a), when the coupling strength $g_f/J$ changes abruptly from 10 (the initial coupling strength) to 11, there is a tiny difference between $\bar{\rho}$ and $\rho_c$. In the cases of $g_f/J=30$ and $g_f/J=50$, the distribution of the time-averaged density matrix $\bar{\rho}$ and canonical ensemble $\rho_c$ are similar, i.e., a linear variation with $E^f_l$, but they are in a bad agreement. On the other hand, further numerical results show that the discrepancy of the distribution between the time-averaged density matrix and the canonical ensemble ones can be diminished remarkably as $g_f/J$ is close to the initial value $g_i/J$, and conversely, the larger the deviation. These results are also related to whether the inverse temperature $\beta_fJ$ approaches the initial inverse temperature $\beta_iJ$. It can also be seen from fig. 1 that the density of states exhibits surprising symmetry no matter how large the value of $g_f/J$ takes after the quantum quench, which also shows a good correspondence with the results obtained from the time-averaged density matrix. It also can be found, in figs. 1(b) and (c), that the density of states located at the edges of the spectrum is small, and the red dots are less. The density of states approaches to the peak at the center of the spectrum, and the red dots are more. In addition, at some fixed locations of the spectrum, the density of states and the time-averaged density matrix will appear or disappear simultaneously as shown in fig. 1(d). In the entire spectrum, there are many local peaks in the density of states, which indicates that the localization properties of photons and atoms are very strong, and the system shows the characteristics of bound states. Based on the analysis mentioned above, one knows that it is a failure to describe the JCH system using a canonical ensemble with a small nearest-neighbor hopping strength, which suggests that the JCH system within the parameters under consideration may be an integrable one.

Zoom In Zoom Out Reset image size

In order to identify whether the system is integrable, we analyze the level spacing distribution P(s) with the spacing s of the neighboring energy levels. The form of the distribution depends on the symmetry properties of the Hamiltonian. The quantum levels of integrable systems are not correlated and levels may cross, so the distribution is Poissonian. In nonintegrable systems, crossings are avoided and the level spacing distribution is given by the Gaussian orthogonal ensembles or the Gaussian unitary ensembles. The results of the level spacing distribution with weak nearest-neighbor hopping strength are shown in figs. 2(a)–(d) for a wide range of values of g/J (condition $J\ll g\ll \omega$ is satisfied), which shows that the level spacing distribution has similar features and approximates to the characteristics of Poisson distribution due to a large number of degenerates and small interval of the system eigenvalues. Thus, this also explains the reasons that the time-averaged density matrix does not fall close to the canonical ensemble after the quantum quench as shown in fig. 1. Actually, the JCH model can be solved analytically for weak nearest-neighbor hopping strength [40], therefore the results obtained here also agree with this conclusion. In order to further confirm whether the system can be integrable, we remove the degenerate energy level and display the Brody distribution in figs. 2(e)–(h).

Zoom In Zoom Out Reset image size

Generally, the expression of the Brody distribution is

$$\begin{equation} P(s)=b(1+\alpha)s^{\alpha}e^{-bs^{\alpha+1}}, \qquad b=\biggl[\Gamma\biggl(\frac{2+\alpha}{1+\alpha}\biggr)\biggr]^{\alpha+1}. \end{equation} \tag{ 11 }$$

In the strongly integrable regime $\alpha=0$ we observe Poissonian statistics while in the completely chaotic one $\alpha=1$. Thus, the system is close to an integrable system with the parameters considered as shown in figs. 2(e)–(h).

In the k = 0 quasimomentum sector for three different system sizes M = 6, 7, 8, we plot the results for observables $\langle a^{\dagger}_{0}a_{0}\rangle$ vs. the eigenstate energy density $E^{f}_{l}/E_{\textit{tot}}$ in figs. 5(a)–(d). The results show the widths of the matrix-element distribution of $\langle a^{\dagger}_{0}a_{0}\rangle$ in the middle of the spectrum become larger as the system size M increases. Therefore, these results violate the validity of ETH (the eigenstate-expectation values will become a smooth and sharp function of the energy density when M tends to infinity).

In addition, the key step to verify the validity of ETH is to analyze the fluctuations of the diagonal elements of observables with the finite-size scaling. The validity of ETH means that the eigenstate-to-eigenstate fluctuations will disappear when M tends to infinity. Therefore, we analyse eigenstate-to-eigenstate fluctuations vs. the size of the system in the bulk of the spectrum. A measure of eigenstate-to-eigenstate fluctuations of diagonal expectation values [45] is

$$\begin{eqnarray} &&Z_l(O)=\langle\Psi^{f}_{l+1}|\widehat{O}|\Psi^{f}_{l+1}\rangle-\langle\Psi^{f}_l|\widehat{O}|\Psi^{f}_l\rangle. \end{eqnarray} \tag{ 15 }$$

In a narrow energy window around $E=E_{\textit{tot}}$ (with a control parameter $\eta\in[0, 1]$), the mean is defined as

$$\begin{eqnarray} &&\langle Z\rangle_{\eta}(O)=||Z_{\eta}||^{-1}\sum_{|\Psi^{f}_l\rangle\in Z_{\eta}}|Z_l(O)|, \end{eqnarray} \tag{ 16 }$$

where $||Z_{\eta}||$ is the number of eigenstates in $Z_{\eta}$, $E_{\textit{tot}}$ and η are defined as follows:

$$\begin{eqnarray} & E_{\textit{tot}}=\langle \Phi_{ini}|\hat{H}_{\textit{fin}}|\Phi_{\textit{ini}}\rangle=\langle \Phi_{\textit{fin}}|\hat{H}_{\textit{fin}}|\Phi_{\textit{fin}}\rangle, \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray} & \frac{E_{\textit{tot}}-E^{f}_{l}}{E_{\textit{tot}}-E_{\textit{min}}}<\eta \quad if \quad E^{f}_{l}<E_{\textit{tot}}, \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray} & \frac{E^{f}_{l}-E_{\textit{tot}}}{E_{\textit{max}}-E_{\textit{tot}}}<\eta \quad if \quad E^{f}_{l}>E_{\textit{tot}}, \end{eqnarray} \tag{ 19 }$$

we fix $\eta=1/3$. Therefore, the corresponding set of eigenstates are limited to the above target energy-density window.

Figure 5(d) shows the finite-size scaling of $\langle Z\rangle_{\eta}(a^{\dagger}_{0}a_{0})$. We observe that the fluctuations hardly decrease as the system size increases. This result violates the validity of ETH, where the fluctuations will decrease exponentially with increasing size. The behavior shown in fig. 5(d) is in accord with one of the diagonal elements of the integrable system [46] and violates the diagonal part of the ETH ansatz in our model.

Zoom In Zoom Out Reset image size

Figures 6(a)–(c) shows the off-diagonal matrix elements of observable $\langle\Psi^{f}_l|a^{\dagger}_{0}a_{0})|\Psi^{f}_n\rangle$ as a function of the eigenenergy difference $\omega=E^{f}_n-E^{f}_l$ (we only plot results for $\omega>0$) in a narrow energy window around $E_{\textit{tot}}$. We observe that the off-diagonal matrix elements of the observable have strong fluctuations, and the fluctuations increase rapidly with the system size. In addition, the running average (black line) of those elements oscillates, but it does not show a downward trend. This is different from the running average of the off-diagonal matrix element that satisfies ETH (the running average is almost flat at small ω and then rapidly decays at larger ω in the nonintegrable case [12]). In fig. 6(d) we present such a plot for the average between the off-diagonal matrix elements of eigenstates with $0.95<(E^f_l+E^f_n)/(2E_{\textit{tot}})<1.05$. The numerical results do not exhibit an excellent agreement with $(N_\varepsilon)^{-1/2}$ behavior, which does not satisfy the nondiagonal part of ETH ansatz.

Zoom In Zoom Out Reset image size