Signatures of Integrability and Exactly Solvable Dynamics in an Infinite-Range Many-Body Floquet Spin System

Abstract

In a recent work Sharma and Bhosale [Phys. Rev. B 109, 014412 (2024)], $N$-spin Floquet model having infinite range Ising interaction was introduced. In this paper, we generalized the strength of interaction to $J$, such that $J=1$ case reduces to the aforementioned work. We show that for $J=1/2$ the model still exhibits integrability for an even number of qubits only. We analytically solve the cases of $6$, $8$, $10$, and $12$ qubits, finding its eigensystem, dynamics of the entanglement for various initial states, and the unitary evolution operator. These quantities exhibit the signature of quantum integrability (QI). For the general case of even-$N>12$ qubits, we conjuncture the presence of QI using the numerical evidences such as spectrum degeneracy, and the exact periodic nature of both the entanglement dynamics and the time-evolved unitary operator. We numerically show the absence of QI for odd $N$ by observing a violation of the signatures of QI. We analytically and numerically find that the maximum value of time-evolved concurrence ($C_{\mbox{max}}$) decreases with $N$, indicating the multipartite nature of entanglement. Possible experiments to verify our results are discussed.

I Introduction

Long-range interactions are encountered in various scientific domains, including statistical physics [1, 2, 3], quantum mechanics [4, 5], cosmology [6, 7, 8, 9], atomic and nuclear physics [10], plasma physics [11], hydrodynamics [12, 13, 14] and condensed matter physics [15]. This interaction can now be replicated in quantum simulator [16, 17, 18] with artificial ion crystals [19, 20, 21], cold atoms in cavities [22], polar molecules [23], dipolar quantum gases [24, 25, 26], Rydberg atoms [27], magnetic atoms [28, 29, 30], nonlinear optical media [31], and solid- state defects [32]. They are also present between genomic elements which can be detected by using the 3C method (chromosome conformation capture) [33] and in adatoms of graphene [34]. Long-range interactions are pervasive and give rise to qualitatively new physics, e.g. emergence of novel quantum phase and dynamical behaviors [35, 36]. Additionally, they play a crucial role in enabling speed up in quantum information processing [37, 38]. Long-range interactions are proven to be very useful in quantum technology applications like quantum heat engine [39], quantum computing [40, 41], quantum metrology [42] and ion trap [43].

Entanglement stands out as perhaps the most distinctive feature of quantum mechanics, giving rise to unique nonlocal correlations that find no counterpart in the classical domain. Long-range interactions have the potential to profoundly influence the dynamics of correlated systems [1]. Due to the breakdown of quasi-locality in long-range interaction various phenomena arise such as it is capable of faster entanglement generation [44, 45, 46] than the Lieb-Robinson bound [47], and dynamical phase transitions [48, 49]. Numerous studies have delved into the examination of entanglement in such system [50, 51, 52, 53, 54, 55, 56, 57, 58]. The time evolution of multipartite entanglement measured using quantum Fisher information (QFI) and scrambling in a spin chain has been studied having these interactions [59]. They have been measured in experiments too [60].

The long-range interaction decays with distance $r$ according to a power law ($1/r^{\alpha}$). Corresponding to different values of $\alpha$, these interactions in the natural system fall into categories such as Rydberg atoms ($\alpha=6$, van der Waals interactions), magnetic atoms ($\alpha=3$, dipole- dipole interaction), Coulomb interactions ($\alpha=1$), atoms coupled to cavities ($\alpha=0$), etc. The case $\alpha=0$ corresponds to a category of infinite-range or all-to-all interactions [61, 62, 63, 23, 64, 65, 66, 67, 68, 69, 70, 71, 72]. Models with $\alpha<d$ are categorized as long-range interaction, where $d$ is the physical dimension of the system [1, 15]. In such interactions, the energy is not extensive [3, 1, 15]. There are models corresponding to infinite range interaction that are integrable, such as the LMG model [73] and the model with Ising interaction in a transverse field [74]. The main theme of this paper revolves around integrability in one such recently introduced model [74].

From the perspective of classical mechanics, integrability can be understood by the connection between the degree of freedom and a sufficient number of constants of motion [75, 76], whereas the quantum integrability [77, 78, 79, 80, 81, 82] is typically associated with an exact solution of models, for example, based on the solution of the Yang-Baxter equation by obtaining transfer matrix [83, 78, 84, 85, 86, 87], and techniques like the Bethe ansatz [88, 89, 90, 91]. Alternatively, it can be identified through other features, for instance, the existence of a set of infinite number of conserved quantities and/or Poissonian-level statistics after taking conserved quantities into account [92, 93]. In various studies, the indication of integrability in systems is revealed through signatures like exact periodicity of entanglement dynamics and time evolution of Floquet operator and degenerated spectra or level crossing [94, 80, 95, 96, 78, 97, 98].

In a very recent work [74], a many-body Floquet spin model having an infinite range Ising interaction was introduced. They have shown that it exhibits quantum integrability. In this Ref. [74], they have analytically calculated the eigensystem, reduced density matrix, time evolution of the unitary operator, and entanglement dynamics for $5$ to $11$ qubits. They have measured the entanglement dynamics using linear entropy and concurrence. For the general case $N>11$ qubits, they resort to numerical methods due to the complexity and fairly large calculation, as analytical solutions pose mathematical challenges. The Hamiltonian was defined as follows:

$$H(t)=\sum_{l<l^{\prime}}\sigma^{z}_{l}\sigma^{z}_{l^{\prime}}+\sum_{n=-\infty}^{\infty}\delta(n-t/\tau)\sum_{l=1}^{N}\sigma^{y}_{l},$$

where the first term represents the Ising interaction with unit strength, while in the second term, $\tau$ is the period at which the magnetic field is applied periodically along the $y$ axis. In various studies, the integrability in system was identified through the signatures, such as the periodicity of entanglement dynamics [99, 98, 96, 97] and that of the Floquet operator dynamics [97], and a highly degenerated spectra [80, 98]. In Ref. [74] they have used the same signatures to show integrability in their model for any value of $N$. In their study, the authors reported that the integrability in the system is found for a specific value of the parameter $\tau=\pi/4$, and for $\tau\neq\pi/4$ it wasn’t. Their findings revealed that the pairwise entanglement using concurrence remains zero for various initial states, indicating the multipartite nature of the entanglement. Furthermore, the signatures of integrability show the same behavior (depending on the parity of $N$) for any value of $N$. The authors have shown that the special case of their model is connected to the one with nearest-neighbor Ising interaction model [99, 96, 100, 101, 98, 102] as well as to the well-known quantum chaotic kicked top (QKT) model for the specific value of parameters [103]. Once it was mapped with QKT, the integrability in the system was limited only up to four qubits [104], whereas, in their work, they generalized it to any $N$ using these signatures [74].

In this work, we have generalized this model to include any value of the strength of Ising interaction to $J$. In the context of integrability, we utilize the same signatures, and for the same value of the $\tau$, as reported in the previous work [74]. We analytically obtained the eigenvalues and eigenvectors for an even number of qubits ranging from $6$ to $12$. We explicitly derive the expression for entanglement measures, such as linear entropy and entanglement entropy for various initial unentangled states. Here, we also find that these quantities exhibit a periodic nature for specific values of $J$, other than one [74]. We also observe the periodic behavior of the time evolution of the unitary operator for the aforementioned cases. However, for cases involving $5$, $7$, $9$, and $11$ qubits, we have numerically demonstrated the absence of these signatures. For the general case, $N>12$, we provide sufficient numerical evidences of Integrability for any even number of qubits and the absence of quantum integrability for any odd qubits by using the absence of the same signatures [74]. We also observe the decay of the maximum value of concurrence with $N$, indicating the increasing multipartite nature of entanglement.

The rest of this paper is organized as follows. In Sec. II, we give a brief introduction to the model investigated. In Sec. III, we have given an exact analytical solution for entanglement measures such as linear entropy, entanglement entropy, and concurrence for six qubits for various initial unentangled states. In Sec. IV, the results on the analytical expression of entanglement and eigenvalues of time evolved unitary operator for $N=8,10$ and $12$ qubits are presented. In Sec. V, we provide sufficient numerical evidence of the signature of quantum integrability for the general case of even-$N>12$ qubits and the absence of the QI for odd-$N>12$ qubits. In Sec. VI, a summary of the results and conclusions is given.

II Model

The generalization of the Hamiltonian model from the Ref. [74] is given as follows:

$$H(t)=H_{I}(h)+\sum_{n=-\infty}^{\infty}\delta(n-t/\tau)~{}H_{k},$$ (1)

where $\delta(t)$ is Dirac delta function and we define,

$$H_{I}={J}\sum_{l<l^{\prime}}\sigma^{z}_{l}\sigma^{z}_{l^{\prime}}~{}~{}\mbox{and}~{}~{}H_{k}=\sum_{l=1}^{N}\sigma^{y}_{l}.$$ (2)

Here, the first term represents the ising interaction with a field strength ${J}$, while the second term corresponds to the periodically applied magnetic field along the $y$-axis, with a period of $\tau$. The case $J=1$ corresponds to the Hamiltonian from the Ref. [74]. The corresponding Floquet operator is given as follows:

$$\mathcal{U} = \exp\left[-i~{}\tau H_{I}(h)\right]\exp\left[-i~{}\tau H_{k}\right]$$ (3)

$$= \exp\left(-i~{}J\tau\sum_{l<l^{\prime}}\sigma^{z}_{l}\sigma^{z}_{l^{\prime}}\right)\exp\left(-i~{}\tau\sum_{l=1}^{N}\sigma^{y}_{l}\right).$$

The Ising interaction in this model is uniform and all-to-all. The nearest-neighbor (NN) interaction model, a special case of our model, has been extensively studied [99, 96, 100, 97, 101, 98, 102]. It exhibits permutation symmetry under the exchange of spins. Due to the presence of symmetries in the model, its effective Hilbert dimension reduces from $2^{N}$ to $N+1$. Scrambling in models similar to ours has been extensively explored in the Refs. [105, 106, 107, 108, 109, 110]. In Ref. [74], they have shown that the model exhibits signatures of quantum integrability for the parameter $\tau=\pi/4$ and ${J}=1$. The main objective of this work is to find the values of parameters in which the model shows quantum integrability. From that perspective, we restricted ourselves to the same $\tau$, while varying the field strength ${J}$. We find that the model shows identical signatures for the parameter ${J}=1/2$. The Floquet operator corresponding to the Eq. (1) with $\tau=\pi/4$ and ${J}=1/2$ is given as follows:

$${\mathcal{U}}=\exp\left(-i\frac{\pi}{8}\sum_{l<l^{\prime}}\sigma^{z}_{l}\sigma^{z}_{l^{\prime}}\right)\exp\left(-i\frac{\pi}{4}\sum_{l=1}^{N}\sigma^{y}_{l}\right).$$ (4)

Throughout this paper, we consistently employ the parameter values $J=1/2$ and $\tau=\pi/4$ unless otherwise stated. In their study [74], the authors reported that this model has a close connection to the very well-known QKT model [103, 111, 112] and shows quantum integrability up to $4$ qubits. In the Refs. [112, 104], this model is extensively studied for a smaller number of qubits ($N=2$, $3$ and $4$) for these parameters. Its Hamiltonian is given by,

$$H_{QKT}(t)=\frac{p}{\tau^{\prime}}\,{\tilde{J}_{y}}+\frac{k}{2j}{\tilde{J}_{z}}^{2}\sum_{n=-\infty}^{\infty}\delta(t-n\tau^{\prime}).$$ (5)

The Floquet operator in the Eq. (4) also has connection with this model, for the parameters $p=\pi/2$, $k=j\pi/2~{}(=N\pi/4)$, $\tau^{\prime}=1$, and using many-qubit transformation $\tilde{J}_{x,y,z}=\sum_{l=1}^{2j}\sigma_{l}^{x,y,z}/2$, where $\sigma_{l}^{x,y,z}$ are the standard Pauli matrices. After this transformation Eqs. (4) and (5) are related by $\exp\left(-i~{}H_{QKT}\right)=\exp\left(-i~{}N{J}\tau\right)\times\mathcal{U}$. The overall phase (global phase) $\exp\left(-i~{}N{J}\tau\right)=\exp\left(-i~{}N\pi/8\right)$ for $J=1/2$ and $\tau=\pi/4$, does not alter the entanglement dynamics. However, the effect of this phase can be observed on the time period of operator dynamics, when it is periodic.

In this work, we will be studying the time evolution from the initial states, that are, localized in spherical phase space. They lie on the unit sphere with spherical coordinates, $(\theta_{0},\phi_{0})$. These states are the standard SU(2) coherent states and given as follows [113, 114]:

$$|\theta_{0},\phi_{0}\rangle=\otimes^{N}\left(\cos(\theta_{0}/2)|0\rangle+e^{-i\phi_{0}}\sin(\theta_{0}/2)|1\rangle\right).$$ (6)

The entanglement dynamics as a function of time can be studied, by evolving these initial states using the Floquet operator $\mathcal{U}$. Our study utilizes various measures such as linear entropy [115], von Neumann entropy [116, 117] and concurrence [118, 119] to quantify the entanglement in the system. The commutation relation $[\mathcal{U},\otimes_{l=1}^{N}\sigma_{l}^{y}]=0$, implies that the system possesses a symmetry, which is an up-down or parity symmetry (for proof see the Supplemental Material of the Ref. [74]). In this work, we utilize the same general basis from Ref. [74] to solve the system analytically for any $N$ number of qubits. The basis when $N$ is odd is given as follows:

$$\ket{\phi_{q}^{\pm}}=\frac{1}{\sqrt{2}}\left(\ket{w_{q}}\pm{i^{\left(N-2q\right)}}\ket{\overline{w_{q}}}\right),0\leq q\leq\dfrac{2j-1}{2};$$ (7)

whereas for even $N$ it is :

$$\begin{split}\ket{\phi_{r}^{\pm}}&=\frac{1}{\sqrt{2}}\left(\ket{w_{r}}\pm{(-1)^{\left(j-r\right)}}\ket{\overline{w_{r}}}\right),0\leq r\leq j-1\\ \mbox{and}&\;\;\ket{\phi_{j}^{+}}=\left({1}/{\sqrt{\binom{N}{j}}}\right)\sum_{\mathcal{P}}\left(\otimes^{j}\ket{0}\otimes^{j}\ket{1}\right)_{\mathcal{P}},\end{split}$$ (8)

where $\ket{w_{q}}=\left({1}/{\sqrt{\binom{N}{q}}}\right)\sum_{\mathcal{P}}\left(\otimes^{q}\ket{1}\otimes^{(N-q)}\ket{0}\right)_{\mathcal{P}}$ and $\ket{\overline{w_{q}}}=\left({1}/{\sqrt{\binom{N}{q}}}\right)\sum_{\mathcal{P}}\left(\otimes^{q}\ket{0}\otimes^{(N-q)}\ket{1}\right)_{\mathcal{P}}$, both being definite particle states [120]. The $\sum_{\mathcal{P}}$ denotes the sum over all possible permutations. These basis states $\ket{\phi_{j}^{\pm}}$ are the eigenstate of parity operator having eigenvalues $\pm 1$ i.e. $\otimes_{l=1}^{N}\sigma_{l}^{y}\ket{\phi_{j}^{\pm}}=\pm\ket{\phi_{j}^{\pm}}$. In this permutation symmetric basis, $\mathcal{U}$ is block-diagonalized in two blocks $\mathcal{U}^{+}$ and $\mathcal{U}^{-}$. Due to block diagonalization, it becomes easy to compute its $n^{th}$ power, which can simplify further analysis. We analytically calculate the entanglement dynamics for the coherent states $|\theta_{0}=0,~{}\phi_{0}=0\rangle$ and $|\theta_{0}=\pi/2,~{}\phi_{0}=-\pi/2\rangle$. These states have special importance when the classical phase space of the QKT is considered [104].

III Exact solution for six qubits

By using the Eq. (8) for $N=6$, the permutation symmetric basis in which $\mathcal{U}$ is block diagonal are given as follows:

$$\ket{\phi_{0}^{\pm}} = \frac{1}{\sqrt{2}}(\ket{w_{0}}\mp\ket{\overline{w_{0}}}),$$ (9)

$$\ket{\phi_{1}^{\pm}} = \frac{1}{\sqrt{2}}(\ket{w_{1}}\pm\ket{\overline{w_{1}}}),$$ (10)

$$\ket{\phi_{2}^{\pm}} = \frac{1}{\sqrt{2}}(\ket{w_{2}}\mp\ket{\overline{w_{2}}}),$$ (11)

$$\ket{\phi_{3}^{+}} = \frac{1}{\sqrt{20}}\sum_{\mathcal{P}}\ket{000111},$$ (12)

where $\ket{w_{0}}=\ket{000000}$, $\ket{\overline{w_{0}}}=\ket{111111}$, $\ket{w_{1}}=\frac{1}{\sqrt{6}}\sum_{\mathcal{P}}\ket{000001}_{\mathcal{P}}$, $\ket{\overline{w_{1}}}=\frac{1}{\sqrt{6}}\sum_{\mathcal{P}}\ket{011111}_{\mathcal{P}}$, $\ket{w_{2}}$=$\frac{1}{\sqrt{15}}\sum_{\mathcal{P}}\ket{000011}_{\mathcal{P}}$ and $\ket{\overline{w_{2}}}$=$\frac{1}{\sqrt{15}}\sum_{\mathcal{P}}\ket{001111}_{\mathcal{P}}$. The seven dimensional ($N+1$) space splits into $4\oplus 3$ subspace on which operators are $\mathcal{U}_{\pm}$. The unitary operator is given by,

$$\mathcal{U}=\begin{pmatrix}\mathcal{U}_{+}&0_{A}\\ 0_{B}&\mathcal{U}_{-}\end{pmatrix},$$ (13)

where  $\mathcal{U}_{+}\left(\mathcal{U}_{-}\right)$ are $4\times 4\left(3\times 3\right)$ dimensional matrices and $0_{A}$ ($0_{B}$) is null matrices of the dimension $4\times 3\left(3\times 4\right)$. The $\mathcal{U}_{+}$ ($\mathcal{U}_{-}$) are written in the positive (negative) parity subspaces $\left\{\phi_{0}^{+},\phi_{1}^{+},\phi_{2}^{+},\phi_{3}^{+}\right\}$ $\left(\left\{\phi_{0}^{-},\phi_{1}^{-},\phi_{2}^{-}\right\}\right)$ respectively and are obtained as follows:

$$\mathcal{U_{+}}=\frac{-e^{\frac{i\pi}{8}}}{2\sqrt{2}}\left(\begin{array}[]{cccc}0&\sqrt{3}&0&\sqrt{5}\\ \sqrt{3}~{}e^{\frac{i\pi}{4}}&0&\sqrt{5}~{}e^{\frac{i\pi}{4}}&0\\ 0&\sqrt{5}&0&-\sqrt{3}\\ -\sqrt{5}~{}e^{\frac{i\pi}{4}}&0&\sqrt{3}~{}e^{\frac{i\pi}{4}}&0\\ \end{array}\right)\mbox{and}$$ (14)

$$\mathcal{U_{-}}=\frac{e^{\frac{i\pi}{8}}}{4}\left(\begin{array}[]{ccc}1&0&\sqrt{15}\\ 0&4~{}e^{\frac{i\pi}{4}}&0\\ \sqrt{15}&0&-1\\ \end{array}\right).$$ (15)

The eigenvalues of the $\mathcal{U_{+}}\left(\mathcal{U_{-}}\right)$ are ${\left\{-(-1)^{1/4},-(-1)^{3/4},(-1)^{1/4}\right\}}\left(\left\{-1,-1,1,1\right\}\right)$ and
the eigenvectors are $\left\{\left[\sqrt{\frac{3}{5}},\sqrt{\frac{3}{5}},\frac{(-1)^{3/8}}{2}\sqrt{\frac{5}{2}},-\frac{(-1)^{3/8}}{2}\sqrt{\frac{5}{2}}\right]^{T},\right.$
$\left.\left[-2(-1)^{1/8}\sqrt{\frac{2}{5}},2(-1)^{1/8}\sqrt{\frac{2}{5}},0,0\right]^{T},\left[1,1,-\frac{(-1)^{3/8}}{2}\sqrt{\frac{3}{2}},\right.\right.$
$\left.\left.\frac{(-1)^{3/8}}{2}\sqrt{\frac{3}{2}}\right]^{T},\left[0,0,1,1\right]^{T}\right\}\left(\left\{\left[-\sqrt{\frac{3}{5}},\sqrt{\frac{5}{3}},0\right]^{T},\left[0,0,1\right]^{T},\right.\right.$
$\left.\left.\left[1,1,0\right]^{T}\right\}\right)$. Thus, the $n$th time evolution of the blocks $\mathcal{U_{\pm}}$ is given as follows:

$$\mathcal{U}_{+}^{n}={\frac{e^{\frac{in\pi}{4}}}{16}\left[\begin{array}[]{cccc}{a_{n}}\left(3+5e^{\frac{in\pi}{2}}\right)&2\sqrt{{6}}~{}b_{n}{e^{\frac{7i\pi}{8}}}&4{i\sqrt{15}~{}e^{\frac{3in\pi}{4}}}\sin\left(\frac{n\pi}{4}\right)\cos\left(\frac{n\pi}{2}\right)&2\sqrt{{10}}~{}b_{n}~{}{e^{(\frac{in\pi}{2}+\frac{3i\pi}{8})}}\\ -2\sqrt{{6}}~{}b_{n}~{}e^{\frac{i\pi}{8}}&8a_{n}&-2\sqrt{{10}}~{}b_{n}~{}{e^{\frac{i\pi}{8}}}&0\\ {4i\sqrt{15}~{}e^{\frac{3in\pi}{4}}}\sin\left(\frac{n\pi}{4}\right)\cos\left(\frac{n\pi}{2}\right)&2\sqrt{{10}}~{}b_{n}{e^{\frac{7i\pi}{8}}}&a_{n}\left(5+3e^{\frac{in\pi}{2}}\right)&-2\sqrt{{6}}~{}b_{n}~{}{e^{(\frac{in\pi}{2}+\frac{3i\pi}{8})}}\\ -2\sqrt{{10}}~{}b_{n}~{}{e^{(\frac{in\pi}{2}+\frac{5i\pi}{8})}}&0&2\sqrt{{6}}~{}b_{n}~{}{e^{(\frac{in\pi}{2}+\frac{i\pi}{8})}}&{8a_{n}~{}e^{\frac{in\pi}{2}}}\\ \end{array}\right]}$$ (16)

$$\mbox{and}~{}~{}~{}~{}\mathcal{U}_{-}^{n}=\frac{e^{\frac{in\pi}{8}}}{8}{\left[\begin{array}[]{ccc}5+3~{}e^{in\pi}&0&\sqrt{15}\left(1-e^{{in\pi}}\right)\\ 0&8~{}e^{\frac{in\pi}{4}}&0\\ \sqrt{15}\left(1-e^{{in\pi}}\right)&0&3+5~{}e^{in\pi}\\ \end{array}\right]},$$ (17)

where $a_{n}$=$1+e^{in\pi}$ and $b_{n}$=$1-e^{in\pi}$. From the Eqs. (16) and (17), we can observe the periodic nature of $\mathcal{U}^{n}_{+}(\mathcal{U}^{n}_{-})$ with period $8(16)$. Hence, the unitary operator shows periodicity with period $16$ i.e. $\mathcal{U}^{n}=\mathcal{U}^{n+16}$. In Refs. [97, 74] it was found that when the Hamiltonian is integrable then the corresponding unitary operator shows a periodic nature. Now, it is straightforward to perform the time evolution of any initial state and consequently study various quantum correlations. We further conducted a detailed analysis on two states, $\ket{0,0}$ and $\ket{\pi/2,-\pi/2}$, deriving exact expressions for the linear entropy, von Neumann entropy of single-qubit reduced density matrix (RDM), and the concurrence between any two qubits.

III.1 Initial state $\ket{000000}=\ket{\theta_{0}=0,\phi_{0}=0}$

Let us first start with the state $|000000\rangle$. The $n$th time evolution of the unitary operator on this state is given by:

$$\ket{\psi_{n}} = \mathcal{U}^{n}|00000\rangle=\mathcal{U}^{n}|w_{0}\rangle=\mathcal{U}^{n}\left(|\phi_{1}^{+}\rangle+|\phi_{1}^{-}\rangle\right)/\sqrt{2}$$

$$= \left(\mathcal{U}_{+}^{n}|\phi_{1}^{+}\rangle+\mathcal{U}_{-}^{n}|\phi_{1}^{-}\rangle\right)/\sqrt{2}$$

$$= \frac{1}{\sqrt{2}}\left(a_{1}\ket{\phi_{0}^{+}}+a_{2}\ket{\phi_{1}^{+}}+a_{3}\ket{\phi_{2}^{+}}+a_{4}\ket{\phi_{3}^{+}}+a_{5}\ket{\phi_{0}^{-}}+a_{6}\ket{\phi_{2}^{-}}\right),$$

where ${a}_{1}$=$\left[\frac{e^{\frac{in\pi}{4}}}{16}\left(3+5~{}e^{\frac{in\pi}{2}}\right)\left(1+e^{in\pi}\right)\right]$, ${a}_{2}$ = $\frac{e^{\frac{i\pi}{8}+\frac{in\pi}{4}}}{4}\sqrt{\frac{3}{2}}\left(-1+e^{in\pi}\right)$, ${a}_{3}=\frac{i\sqrt{15}~{}e^{in\pi}}{8}\left(\sin\left(\frac{n\pi}{4}\right)-\sin\left(\frac{3n\pi}{4}\right)\right)$, ${a}_{4}=\frac{e^{\left(\frac{5i\pi}{8}+\frac{3in\pi}{4}\right)}}{4}\sqrt{\frac{5}{2}}\left(-1+e^{in\pi}\right)$,  ${a}_{5}=\frac{e^{\frac{in\pi}{8}}}{8}\left(5+3e^{{in\pi}}\right)$ and ${a}_{6}=\sqrt{15}~{}e^{\frac{in\pi}{8}}\left(1-e^{{in\pi}}\right)/8$. From $\ket{\psi_{n}}$, we will obtain the RDM of a single qubit $\left(\rho_{1}(n)=\mbox{Tr}_{\neq 1}\left(\ket{\psi_{n}}\bra{\psi_{n}}\right)\right)$ and the two qubits $\left(\rho_{12}(n)=\mbox{Tr}_{\neq 1,2}\left(\ket{\psi_{n}}\bra{\psi_{n}}\right)\right)$. Thus we get $\rho_{1}(n)$, which is given as follows:

$$\rho_{1}(n)=\frac{1}{4}\left(\begin{array}[]{cc}2+a_{n}&w_{n}\\ w_{n}^{*}&2-a_{n}\\ \end{array}\right),$$ (18)

$$\mbox{where}~{}a_{n} = \frac{1}{4}\cos\left(\frac{n\pi}{8}\right)\cos\left(\frac{n\pi}{2}\right)\left[2-4\cos\left(\frac{n\pi}{4}\right)+5\cos\left(\frac{n\pi}{2}\right)+5\cos(n\pi)\right]\mbox{ and}$$

$$w_{n} = {(-1)^{1/8}~{}e^{-\frac{13}{8}in\pi}}\left(-1+e^{in\pi}\right)\left[(5+5i)\sqrt{2}+5\left(-1+(-1)^{1/4}\right)e^{\frac{in\pi}{8}}+5~{}i\left((-1+i)+\sqrt{2}\right)e^{\frac{9in\pi}{8}}\right.$$

$$\left.-2~{}e^{\frac{3in\pi}{4}}-10(-1)^{1/4}e^{in\pi}-\right.\left.10~{}ie^{\frac{5in\pi}{4}}-2~{}(-1)^{3/4}e^{\frac{3in\pi}{2}}+\left((1+5i)-(2+3i)\sqrt{2}\right)\left(e^{\frac{5in\pi}{8}}+e^{\frac{13in\pi}{8}}\right)\right.$$

$$\left.+10~{}e^{\frac{7in\pi}{4}}+5~{}\left(-i+(-1)^{3/4}\right)e^{\frac{17in\pi}{8}}+10~{}(-1)^{3/4}e^{\frac{in\pi}{2}}+10~{}i~{}e^{\frac{9in\pi}{4}}\right]\Big{/}{64}.$$

The eigenvalues of $\rho_{1}(n)$ are $\lambda_{n}$ and $1-\lambda_{n}$ where $\lambda_{n}$=
$\frac{1}{2}-\frac{1}{16}\left\{17+8\left(2+\sqrt{2}\right)\cos\left(\frac{n\pi}{4}\right)+8\left(2-\sqrt{2}\right)\cos\left(\frac{3n\pi}{4}\right)+\right.$
$\left.15\left[\cos\left(\frac{n\pi}{2}\right)-\sin\left(\frac{n\pi}{2}\right)\right)-8\sqrt{2}\left(\sin\left(\frac{n\pi}{4}\right)+\sin\left(\frac{3n\pi}{4}\right)\right]\right\}^{\frac{1}{2}}$. Using these eigenvalues we can obtain the linear entropy [115] given by $2\lambda_{n}(1-\lambda_{n})$ and is plotted in the Fig. 1. The entanglement entropy [116, 117] can also be calculated using $-(\lambda_{n}\log\lambda_{n}+(1-\lambda_{n})\log(1-\lambda_{n}))$ and is plotted in the same figure. It can be shown from the expressions and the figure that both, the linear and entanglement entropy have a periodic nature with period eight i.e. $S(n+8)=S(n)$. This periodic nature has been previously observed in integrable systems, particularly those involving periodically kicked spin chains Refs. [99, 96, 98, 74]. For this particular state, we provided analytical proof that the entanglement content remains unchanged for consecutive odd and even values of $n$ [121] i.e. $S_{(0,0)}^{(6)}(2n-1)=S_{(0,0)}^{(6)}(2n)$, which is shown in the Fig. 1. The linear and entanglement entropy attains its maximum upper bound value $0.5$ and $\ln 2\approx 0.6932$ for the parameter $\tau=\pi/4$ and ${J}=1/2$, which is also observed in the previous study for this model [74] with same $\tau$ but ${J}=1$. We can observe from the figure that the entanglement content takes more time ($n=3$) to reach the maximum upper bound value compared to the previous work ($n=1$) in Ref. [74].

Figure 1: Correlations (analytical values) using linear entropy (circles), entanglement entropy (squares) are plotted for (a) $6$ qubit (b) $8$ qubit (c)$10$ qubit and (d) $12$ qubit for the initial state $\otimes^{N}{\ket{0}}_{y}$.

III.1.1 Concurrence

The linear entropy measures the entanglement between a single qubit with the rest of the qubits in a pure state, while the concurrence measures the entanglement between any pair of qubits within the system (pure or mixed). Regardless of which particular qubits are chosen, the presence of permutation symmetry in the state results in only one concurrence value [118, 119]. The concurrence is given by,

$$\mathcal{C}(\rho_{12})=\text{max}\left(0,\sqrt{\lambda_{1}}-\sqrt{\lambda_{2}}-\sqrt{\lambda_{3}}-\sqrt{\lambda_{4}}\right),$$ (19)

where $\lambda_{l}$ are eigenvalues in decreasing order of $(\sigma_{y}\otimes\sigma_{y})\rho_{12}(\sigma_{y}\otimes\sigma_{y})\rho_{12}^{*}$, where $\rho_{12}^{*}$ is complex conjugation in the standard ($\sigma_{z}$) basis. The two-qubit RDM, $\rho_{12}$, is given as follows:

$$\rho_{12}(n)={\frac{1}{4}\left(\begin{array}[]{cccc}b_{m}&a_{m}&a_{m}&d_{m}^{*}\\ a_{m}^{*}&e_{m}&e_{m}&a_{m}^{*}\\ a_{m}^{*}&e_{m}&e_{m}&a_{m}^{*}\\ d_{m}&a_{m}&a_{m}&f_{m}\\ \end{array}\right)},$$ (20)

where the coefficients are,

$$b_{m} = \frac{1}{16}\left[22+6\cos\left(\frac{n\pi}{8}\right)+10\cos\left(\frac{3n\pi}{8}\right)+4\cos\left(\frac{n\pi}{2}\right)\right.$$

$$\left.+10\cos\left(\frac{5n\pi}{8}\right)+6\cos\left(\frac{7n\pi}{8}\right)+6\cos(n\pi)\right],$$

$$f_{m} = \frac{1}{4}\left[58+94\cos\left(\frac{n\pi}{8}\right)+78\cos\left(\frac{n\pi}{4}\right)+62\cos\left(\frac{3n\pi}{8}\right)\right.$$

$$\left.+56\cos\left(\frac{n\pi}{2}\right)+46\cos\left(\frac{5n\pi}{8}\right)+46\cos\left(\frac{3n\pi}{4}\right)\right.$$

$$\left.+26\cos(n\pi)+46\cos\left(\frac{7n\pi}{8}\right)\right]\sin^{2}\left(\frac{n\pi}{16}\right),$$

$$d_{m} = \frac{1}{16}\left[4\cos\left(\frac{n\pi}{2}\right)-6\cos(n\pi)-{2~{}i}\left(i-4~{}i\sin\left(\frac{n\pi}{2}\right)+\right.\right.$$

$$\left.\left.\sin\left(\frac{n\pi}{8}\right)+\sin\left(\frac{3n\pi}{8}\right)-\sin\left(\frac{5n\pi}{8}\right)-\sin\left(\frac{7n\pi}{8}\right)\right)\right],$$

$$e_{m} = \frac{1}{8}\left[3\sin\left(\frac{n\pi}{4}\right)+\sin\left(\frac{3n\pi}{4}\right)\right]^{2}~{}~{}\mbox{and}$$

$$a_{m} = (-1)^{3/8}e^{-\frac{13}{8}in\pi}\left(-1+e^{in\pi}\right)\left[6\sqrt{2}+e^{\frac{in\pi}{8}}\left((2-2i)\right.\right.$$

$$\left.\left.-3\sqrt{2}+6~{}i\sqrt{2}e^{\frac{3in\pi}{8}}+4\left(-2+\sqrt{4-3~{}i}\right)\cos\left(\frac{n\pi}{2}\right)\right.\right.$$

$$\left.\left.-(4-4~{}i)~{}\left(e^{\frac{5in\pi}{8}}-e^{\frac{-3in\pi}{8}}\right)+\left(4-6e^{\frac{i\pi}{4}}\right)~{}e^{in\pi}\right.\right.$$

$$\left.\left.-6\sqrt{2}~{}e^{\frac{7in\pi}{8}}-(4+4~{}i)~{}\left(e^{\frac{9in\pi}{8}}-e^{\frac{in\pi}{8}}\right)+2~{}i\sqrt{2}~{}\right.\right.$$

$$\left.\left.e^{\frac{-5in\pi}{8}}+i(-2+2~{}i)-3~{}i\sqrt{2}\right)\right]\Big{/}{64\sqrt{2}}.$$

Figure 2: The concurrence (circles) are plotted for (a) $6$ qubit (b) $8$ qubit (c) $10$ qubit and (d) $12$ qubit for the initial state $\otimes^{N}{\ket{0}}_{y}$.

By evaluating the eigenvalues of $(\sigma_{y}\otimes\sigma_{y})\rho_{12}(\sigma_{y}\otimes\sigma_{y})\rho_{12}^{*}$, we obtain fairly long expressions (please refer supplemental material [121] for their numerical values). Using these values in the Eq. (19), the computed concurrence values are plotted in Fig. 2. Similarly, concurrence is also periodic in nature and it remains the same for consecutive odd and even values of $n$ [121], which can be seen in Fig. 2, whereas in Ref. [74] the pairwise concurrence remains zero for all the values of $n$. Thus all the entanglement measure quantities are periodic in nature which is the signature of quantum integrability [99, 98, 96, 97, 74].

III.2 Initial state $\ket{++++++}$= $\ket{\theta_{0}=\pi/2,\phi_{0}=-\pi/2}$.

Let us now focus on another state $\ket{++++++}_{y}$, where $|+\rangle_{y}=\frac{1}{\sqrt{2}}(|0\rangle+i|1\rangle)$ is an eigenstate of $\sigma_{y}$ with eigenvalue $+1$. The evolution of this state is entirely confined to the positive parity subspace of the seven-dimensional permutation symmetric space of six-qubits hilbert space. It can be expressed as:

$$\otimes^{6}{\ket{+}}_{y}=\frac{1}{4\sqrt{2}}\ket{\phi_{0}^{+}}+\frac{i\sqrt{3}}{4}\ket{\phi_{1}^{+}}-\frac{\sqrt{15}}{4\sqrt{2}}\ket{\phi_{2}^{+}}-\frac{i\sqrt{5}}{4}\ket{\phi_{3}^{+}}.$$ (21)

The state $\ket{\psi_{n}}$ can than be obtained by applying the $n$ iterations of unitary operator $\mathcal{U}$ on it as follows:

$$\ket{\psi_{n}} = \mathcal{U}_{+}^{n}\ket{++++++}$$

$$= e^{{\frac{in\pi}{4}}}\left({\alpha}_{n}^{\prime}\ket{\phi_{0}^{+}}+{\beta}_{n}^{\prime}\ket{\phi_{1}^{+}}+{\gamma}_{n}^{\prime}\ket{\phi_{2}^{+}}+\zeta_{n}\ket{\phi_{3}^{+}}\right),$$

where the coefficients can be expressed as,

$${\alpha}_{n}^{\prime} = \frac{e^{-\frac{3}{4}in\pi}}{{16\sqrt{2}}}\left[3\left(-1+e^{\frac{3i\pi}{8}}\right)+5\left(1+e^{\frac{7i\pi}{8}}\right)e^{\frac{in\pi}{2}}\right.$$

$$\left.-3\left(1+e^{\frac{3i\pi}{8}}\right)e^{in\pi}-5\left(-1+e^{\frac{7i\pi}{8}}\right)e^{\frac{3in\pi}{2}}\right],$$

$$~{}{\beta}_{n}^{\prime} = \frac{\sqrt{3}~{}e^{-\frac{3}{4}in\pi}}{8}\left[i-e^{\frac{i\text{ }\pi}{8}}+\left(i+e^{\frac{i\pi}{8}}\right)e^{in\pi}\right],$$

$${\gamma}_{n}^{\prime} = \frac{i~{}e^{-\frac{3}{4}in\pi}}{16}\sqrt{\frac{15}{2}}\left[1-e^{\frac{3i\text{ }\pi}{8}}+\left(i+e^{\frac{3i\text{ }\pi}{8}}\right)e^{\frac{in\pi}{2}}\right]$$

$$\left[i+(1+i)e^{\frac{in\pi}{2}}+e^{in\pi}\right]~{}~{}\mbox{and}$$

$$~{}\zeta_{n} = \frac{\sqrt{5}~{}e^{-\frac{1}{4}in\pi}}{8}\left[-i+e^{\frac{5i\text{ }\pi}{8}}-\left(i+e^{\frac{5i\text{ }\pi}{8}}\right)e^{in\pi}\right].$$

The $\rho_{1}(n)$, is given as follows:

$$\rho_{1}(n)={\frac{1}{2}\left[\begin{array}[]{cc}1&Z_{n}\\ Z_{n}^{*}&1\\ \end{array}\right]},$$ (23)

where  $Z_{n}={-i}\left[3(1+\cos(n\pi))+\sqrt{2}\left(1-\cos(n\pi)\right)+10\right.$
$\left.\cos\left(\frac{n\pi}{2}\right)\right]\Big{/}{16}$. The eigenvalues of $\rho_{1}(n)$ are $\lambda_{n}$ and $1-\lambda_{n}$
where $\lambda_{n}=\left[8-\sqrt{18+30\cos\left(\frac{n\pi}{2}\right)+16\cos(n\pi)}\right]\Big{/}16$.

Figure 3: Correlations using linear entropy (circles), and entanglement entropy (squares) are plotted for (a) $6$ qubit (b) $8$ qubit (c) $10$ qubit and (d) $12$ qubit for the initial state $\otimes^{N}{\ket{+}}_{y}$.

The linear entropy and entanglement entropy can be calculated by using eigenvalues ($\lambda_{n}$) and are plotted in Fig. 3. The entanglement content is periodic in nature having a period of $4$.

III.2.1 Concurrence

The $\rho_{12}(n)$, for this state is given as follows:

$$\rho_{12}(n)={\frac{1}{2}\left(\begin{array}[]{cccc}h_{1}&h_{4}&h_{4}&h_{2}\\ h_{4}^{*}&h_{3}&h_{3}&h_{6}\\ h_{4}^{*}&h_{3}&h_{3}&h_{6}\\ h_{2}^{*}&h_{6}^{*}&h_{6}^{*}&h_{1}\\ \end{array}\right)},$$ (24)

where the coefficients are,

$$h_{1} = \frac{1}{8}\left[5-\left(\cos\left(\frac{n\pi}{2}\right)+\sin\left(\frac{n\pi}{2}\right)\right)\right],$$

$$h_{2} = \frac{-1}{8}\left[1+3\cos\left(\frac{n\pi}{2}\right)-\sin\left(\frac{n\pi}{2}\right)\right],$$

$$h_{3} = \frac{1}{8}\left[3+\cos\left(\frac{n\pi}{2}\right)+\sin\left(\frac{n\pi}{2}\right)\right],$$

$$h_{6} = \frac{(1-i)}{64}\left[3(1+i)(1+\cos(n\pi))+4\sqrt{2}(1-\cos(n\pi))\right.$$

$$\left.+20\left(1+i\right)\cos\left(\frac{n\pi}{2}\right)\right]\mbox{and}$$

$$h_{4} = \frac{-(1+i)}{64}\left[3(1+i)(1+\cos(n\pi))+4\sqrt{2}(1-\cos(n\pi))\right.$$

$$\left.+20\left(1+i\right)\cos\left(\frac{n\pi}{2}\right)\right].$$

Figure 4: The concurrence (circles) are plotted for (a) $6$ qubit and (b) $8$ qubit (c) $10$ qubit and (d) $12$ qubit for the initial state $\otimes^{N}{\ket{+}}_{y}$.

The eigenvalues of $(\sigma_{y}\otimes\sigma_{y})\rho_{12}(\sigma_{y}\otimes\sigma_{y})\rho_{12}^{*}$ are $\left\{\frac{1}{4}\sin^{4}\left(\frac{n\pi}{4}\right),\frac{1}{256}\left[23-6\cos\left(\frac{n\pi}{2}\right)-17\cos\left({n\pi}\right)\pm 4\right.\right.$
$\left.\left.\left[\cos^{2}\left(\frac{n\pi}{2}\right)(-5+\cos\left({n\pi}\right))\left((-41+16\sqrt{2})\cos\left({n\pi}\right)-109\right.\right.\right.\right.$
$\left.\left.\left.\left.+16\sqrt{2}+2(-69+16\sqrt{2})\cos\left(\frac{n\pi}{2}\right)\right)\sin^{4}\left(\frac{n\pi}{4}\right)\right]^{\frac{1}{2}}\right],0\right\}$. The computed concurrence values are plotted in the Fig. 4. It can be observed that it is periodic in nature with a period of $4$, whereas in Ref. [74] the pairwise concurrence remains zero. The periodic nature of entanglement dynamics and operator dynamics are the signature of integrability [74], hence for $6$ qubit the model is integrable for these parameters as well.

IV Exact solutions for eight, ten and twelve qubits

We followed a similar process from section III to solve the analytically cases involving $N$= $8$, $10$, and $12$ qubits at the particular value of $J$=$1/2$ and $\tau$=$\pi/4$. We created distinct tables for the states $\otimes^{N}|0\rangle$ (refer Tables 1 and 2) and $\otimes^{N}|+\rangle$ (refer Table 3) to show the analytical expression for the eigenvalues and linear entropy for the qubits $8$, $10$ and $12$. Using these expressions, it can be shown that they exhibit a periodic behavior. For the initial state $\otimes^{N}\ket{0}$ and the cases of $6$ and $10$ qubits period is $8$, whereas for $8$ and $12$ qubits it is $24$. This can be observed in the Fig. 1. We observe that each correlation is same for consecutive odd and even values of $n$ [104]. We have plotted the analytically obtained concurrence in Fig. 2. From the values and the plot, we observe that it shows a periodic nature. For $N=6$ and $10$ ($8$ and $12$), we find $T=4(12)$. The maximum value of concurrence ($C_{max}$) for $N$= $6$, $8$, $10$ and $12$ is $0.05618622$, $0.029509$, $0.0038762$ and $0.0019455$ respectively (All values are close upto $12$ decimal places). Reduction of $C_{\mbox{max}}$ with $N$ shows an increased multipartite nature of entanglement. The signatures of QI are obtained by limiting $n$ to positive integers only. The data points in all the figures of this paper correspond to $n\in\mathcal{N}^{+}$.

The Table 3 presents the analytical expressions for the eigenvalues and linear entropy for the initial state $\otimes^{N}|+\rangle$. Using these eigenvalues expression for entanglement entropy can be obtained easily (not shown here). These correlations are plotted in Fig. 3. It can be seen that they show periodic behavior. For the cases $6$ and $10$ ($8$ and $12$) qubits the period is $4$ ($12$). The values of $C_{max}$ for $N$=$6$, $8$, $10$ and $12$ are $0.0669873$, $0.0295085$, $0.00392163$ and $0.0019455$ respectively. Thus, for this state also $C_{\mbox{max}}$ decreases with $N$, displaying the multipartite nature of entanglement.

For this case also we study the dynamics of the floquet operator in which we observe that it shows periodic nature with time i.e. ${\mathcal{U}}^{n+T}={\mathcal{U}}^{n}$, where $n\geq 1$ and $T$ is the period. For qubits $N$= $8$, $10$ and $12$, the period is $48$. Through analytical calculations, we determine the eigenvalues of $\mathcal{U}$ and identify degeneracy among them (see Supplemental Material [121]).

Figure 5: (Color online) Correlations using linear entropy (black), and entanglement entropy (red) are plotted for (a) $5$ qubit (b) $7$ qubit (c) $9$ qubit and (d) $11$ qubit for initial state $\otimes^{N}{\ket{0}}_{y}$.

Figure 6: (Color online) Correlations using linear entropy (black), and entanglement entropy (red) are plotted for (a) $5$ qubit (b) $7$ qubit (c) $9$ qubit and (d) $11$ qubit for the initial state $\otimes^{N}{\ket{+}}_{y}$.

Figure 7: The concurrence (circles) are plotted for (a) $5$ qubit and (b) $7$ qubit (c) $9$ qubit and (d) $11$ qubit for the initial state $\otimes^{N}{\ket{0}}_{y}$.

Figure 8: The concurrence (circles) are plotted for (a) $5$ qubit and (b) $7$ qubit (c) $9$ qubit and (d) $11$ qubit for the initial state $\otimes^{N}{\ket{+}}_{y}$.

For the cases $5$, $7$, $9$, and $11$ qubits, due to complexity and fairly large expressions, we restrict ourselves to numerical results instead of analytical ones. We observed that the entanglement dynamics for both the initial states do not exhibit periodic behavior (checked for $n$ up to $1000$), which is shown in Figs. 5, 6, 7 and 8 (plotted for $n$ up to $500$). We also study the Floquet operator dynamics in which we observe that it does not exhibit periodic behavior for large values of $n$, as shown in Fig. 12. The periodic behavior of entanglement dynamics, and floquet operator dynamics, and degenerated spectrum are the signature of integrability (as explained in the introduction). In the previous study Ref. [74], they have shown, that this model is integrable for the parameters $J=1$ and $\tau=\pi/4$ using these signatures, whereas we observe these signatures only for an even $N$ ranging from $6$ to $12$. Thus, for the parameter $J=1/2$ and $\tau=\pi/4$, the model shows integrability for $6$, $8$, $10$, and $12$ qubits and absence of integrability for odd $N$ ranging from $5$ to $11$.

Table 1: The eigenvalues for the initial state $\otimes^{N}\ket{0}$ and various value of $N$.

$N$	Eigenvalues ($\lambda$)
8	$\lambda=\frac{1}{2}\pm\frac{1}{96}\left\{\left(5+7\cos\left(\frac{2n\pi}{3}\right)\right)^{2}\left[\cos^{2}\left(\frac{3n\pi}{8}\right)+8\cos\left(\frac{3n\pi}{8}\right)\cos\left(\frac{5n\pi}{8}\right)+4\cos^{2}\left(\frac{5n\pi}{8}\right)\right]+\sin^{2}\left(\frac{n\pi}{2}\right)\left[9\cos^{2}\left(\frac{1}{8}(3+4n)\pi\right)+147\right.\right.$ $\left.\left.\left.~{}~{}~{}~{}(\cos^{2}\left(\frac{1}{24}(9+5n)\pi\right)+\cos^{2}\left(\frac{1}{24}(9+13n)\pi\right)\right)+\sin\left(\frac{1}{24}(9+5n)\pi\right)\left(-98\sqrt{3}\cos\left(\frac{1}{24}(9+13n)\pi\right)+49\sin\left(\frac{1}{24}(9+5n)\pi\right)-\right.\right.\right.$ $\left.\left.\left.~{}~{}~{}~{}280\sin\left(\frac{1}{8}(3+7n)\pi\right)\right)+49\sqrt{3}\sin\left(\frac{1}{12}(9+5n)\pi\right)+280\sqrt{3}\cos\left(\frac{1}{24}(9+37n)\pi\right)\sin\left(\frac{1}{8}(3+7n)\pi\right)+400\sin^{2}\left(\frac{1}{8}(3+7n)\pi\right)-\right.\right.$ $\left.\left.~{}~{}~{}~{}\sin\left(\frac{1}{24}(9+13n)\pi\right)\left(98\sin\left(\frac{1}{24}(9+5n)\pi\right)-280\sin\left(\frac{1}{8}(3+7n)\pi\right)-49\sin\left(\frac{1}{24}(9+13n)\pi\right)\right)-14\cos\left(\frac{1}{24}(9+5n)\pi\right)\right.\right.$ $\left.\left.~{}~{}~{}\left(21\cos\left(\frac{1}{24}(9+37n)\pi\right)+\sqrt{3}\left(20\sin\left(\frac{1}{8}(3+7n)\pi\right)+7\sin\left(\frac{1}{24}(9+13n)\pi\right)\right)\right)-49\sqrt{3}\sin\left(\frac{1}{12}(9+13n)\pi\right)+\sin\left(\frac{1}{8}(\pi+4n\pi)\right)\right.\right.$ $\left.\left.~{}~{}~{}\left(18\cos\left(\frac{1}{8}(3+4n)\pi\right)+9\sin\left(\frac{1}{8}(\pi+4n\pi)\right)\right)\right]\right\}^{\frac{1}{2}}$.
10	$\lambda=0.5\pm 0.5\cos(n\pi)\left[0.2509765625+0.426776695296637\cos(\frac{n\pi}{4})+0.073223304703364\cos(\frac{3n\pi}{4})-0.176776695296637\right.$ $~{}~{}~{}~{}~{}~{}\left.\left(\sin(\frac{n\pi}{4})+\cos(\frac{3n\pi}{4})\right)+0.2490234375\left(\cos(\frac{n\pi}{2})-\sin(\frac{n\pi}{2})\right)\right]^{\frac{1}{2}}$.
12	$\lambda=0.5\pm\left[0.041748046875+0.040916791185\cos\left(\frac{n\pi}{12}\right)+0.035634186961\left(\cos\frac{n\pi}{4}\right)+0.06243896484375\cos\left(\frac{n\pi}{3}\right)\right.$ $~{}~{}~{}~{}\left.+0.026199786034\cos\left(\frac{5n\pi}{12}\right)+0.01542619052857\cos\left(\frac{7n\pi}{12}\right)+0.02081298828125\cos\left(\frac{2n\pi}{3}\right)+0.0061138599142\cos\left(\frac{3n\pi}{4}\right)\right.$ $~{}~{}~{}~{}\left.+0.00070918537814\cos\left(\frac{11n\pi}{12}\right)-0.0053867977527\left(\sin\left(\frac{n\pi}{12}\right)+\sin\left(\frac{11n\pi}{12}\right)\right)-0.0147601635233\left(\sin\left(\frac{n\pi}{4}\right)+\sin\left(\frac{3n\pi}{4}\right)\right)\right.$ $~{}~{}~{}~{}\left.-0.036049153161\left(\sin\left(\frac{n\pi}{3}\right)+\sin\left(\frac{2n\pi}{3}\right)\right)-0.020103802903\left(\sin\left(\frac{5n\pi}{12}\right)+\sin\left(\frac{7n\pi}{12}\right)\right)\right]^{\frac{1}{2}}$.

Table 2: The linear entropy ($S(n)$) for the initial state $\otimes^{N}\ket{0}$ and various values of $N$.

$N$	$S(n)$
8	$S_{(0,0)}^{8}(n)=\frac{1}{2}\left\{1-\frac{1}{576}\left[5+7\cos\left(\frac{2n\pi}{3}\right)\right]^{2}\left[\cos\left(\frac{3n\pi}{8}\right)+\cos\left(\frac{11n\pi}{8}\right)\right]^{2}+\sin^{2}\left(\frac{n\pi}{2}\right)\left[-\left(18-9\sqrt{2}\right)\cos(n\pi)+378\sin\left(\frac{1}{12}(3+5n)\pi\right)\right.\right.$ $\left.\left.\hskip 42.67912pt+560\left(\cos\left(\frac{2n\pi}{3}\right)+\sqrt{3}\sin\left(\frac{2n\pi}{3}\right)\right)+196\left(\cos\left(\frac{4n\pi}{3}\right)-\sqrt{3}\sin\left(\frac{4n\pi}{3}\right)\right)+378\sqrt{3}\sin\left(\frac{1}{12}(9+13n)\pi\right)+9\left(-90+\sqrt{2}\right)\right.\right.$ $\left.\left.\hskip 42.67912pt-378\sqrt{3}\sin\left(\frac{1}{12}(9+5n)\pi\right)+378\sin\left(\frac{1}{12}(3+13n)\pi\right)-792\sin\left(\frac{1}{4}(\pi+7n\pi)\right)\right]\right\}$
10	$S_{(0,0)}^{(10)}(n)=0.37451171875-(0.21338834765)\cos\left(\frac{n\pi}{4}\right)-0.036611652352\cos\left(\frac{3n\pi}{4}\right)+0.088388347649\left(\sin\left(\frac{n\pi}{4}\right)+\sin\left(\frac{3n\pi}{4}\right)\right)$ $~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-0.12451171875\left(\cos\left(\frac{n\pi}{2}\right)-\sin\left(\frac{n\pi}{2}\right)\right)$.
12	$S_{(0,0)}^{(12)}(n)=\frac{1}{16}\left\{8-0.1953125e^{-in\pi}\left(1\,+e^{in\pi}\right)^{2}\left[\cos\left(\frac{3n\pi}{8}\right)+1.1\cos\left(\frac{25n\pi}{24}\right)+1.1\cos\left(\frac{41n\pi}{24}\right)\right]^{2}-0.1667096466\left[\cos\left(\frac{n\pi}{8}\right)-\right.\right.$ $\left.\left.\hskip 45.52458pt\cos\left(\frac{7n\pi}{8}\right)-0.4142135623731\left(\sin\left(\frac{n\pi}{8}\right)+\sin\left(\frac{7n\pi}{8}\right)\right)+(1-\cos(n\pi))\left(-0.944591414368\cos\left(\frac{29n\pi}{24}\right)-0.15540858564\right.\right.\right.$ $\left.\left.\left.\hskip 44.10185pt\cos\left(\frac{13n\pi}{24}\right)-1.180445403469\sin\left(\frac{13n\pi}{24}\right)+0.724810484858\sin\left(\frac{29n\pi}{24}\right)\right)\right]^{2}\right\}$.

Table 3: The eigenvalues and linear entropy ($S(n)$) for the initial state $\otimes^{N}\ket{+}$ and various values of $N$.

$N$	Eigenvalues	$S(n)$
8	$\frac{1}{2}\pm\frac{1}{48}\left[10\cos\left(\frac{n\pi}{4}\right)+7\left(\cos\left(\frac{5n\pi}{12}\right)+\cos\left(\frac{11n\pi}{12}\right)\right)\right]$	$S_{(\pi/2,-\pi/2)}^{(8)}(n)=\frac{1}{2}\left[1-\frac{1}{144}\cos^{2}\left(\frac{n\pi}{4}\right)\left(5+7\cos\left(\frac{2n\pi}{3}\right)\right)^{2}\right]$.
10	$\frac{1}{2}\pm\frac{1}{64}\left[\sqrt{258+510\cos\left(\frac{n\pi}{2}\right)+256\cos(n\pi)}\right]$	$S_{(\pi/2,-\pi/2)}^{(10)}(n)=\frac{1}{2}\left\{1-\frac{1}{4096}\left[17\left(1+\cos(n\pi)\right)+\sqrt{2}\left(1-\cos(n\pi)\right)\right.\right.$ $\left.\left.\hskip 56.9055pt+30\cos\left(\frac{n\pi}{2}\right)\right]^{2}\right\}$.
12	$\frac{1}{2}\pm\frac{1}{64}\left[11\cos\left(\frac{n\pi}{12}\right)+10\cos\left(\frac{3n\pi}{4}\right)+11\cos\left(\frac{17n\pi}{12}\right)\right]$	$S_{(\pi/2,-\pi/2)}^{(12)}(n)=\frac{1}{2}-\frac{1}{2048}\left[11\cos\left(\frac{n\pi}{12}\right)+10\cos\left(\frac{3n\pi}{4}\right)+11\cos\left(\frac{17n\pi}{12}\right)\right]^{2}$.

V Results for general $N$ qubits.

Using our procedure, in principle one can solve for the case of any finite $N$, but various expressions for the eigensystem, floquet operator, and its $n^{th}$ power become more cumbersome. Therefore, for $N>12$ we restrict ourselves only to the numerical methods. Here, we use the same signatures used to establish integrability for $N=6,~{}8,~{}10$ and $12$ qubits. Surprisingly, we find exactly the same behavior of these signatures for the said parameters ($J=1/2$ and $\tau=\pi/4$). Particularly, we find that the entanglement dynamics is periodic such that for the initial states $\otimes^{N}{\ket{0}}_{y}$ ($\otimes^{N}{\ket{+}}_{y}$) and $N=4m+2$ the period is $8$ ($4$) respectively, whereas for $N=4m+4$ it is $24$($12$), where $m\in\left\{0,1,2,...\right\}$. The results are shown in Fig. 9. We numerically find that $C_{\mbox{max}}$ tends to zero as $N$ increases for both the initial states. Thus, the entanglement becomes multipartite in nature with $N$. For any odd $N>11$ with the same parameter, we find that the entanglement dynamics is not periodic at all, which is shown in Fig. 10. We observed the same non-periodic behavior for $n$ up to $1000$ (results not shown here).

Figure 9: Correlations using linear entropy (circles), entanglement entropy (squares), and concurrence (upper triangle) are plotted for ((a) $1002$ qubit, (b) $1000$ qubit) for the initial states $\otimes^{N}{\ket{0}}_{y}$ and ((c) $1002$ qubit, (d) $1000$ qubit) for the initial states $\otimes^{N}{\ket{+}}_{y}$

Figure 10: (Color online) Correlations using linear entropy (black) and entanglement entropy (red) are plotted for (a) $101$ qubit, (b) $203$ qubit for the initial states $\otimes^{N}{\ket{0}}_{y}$ and (c) $101$ qubit, (d) $203$ qubit for the initial states $\otimes^{N}{\ket{+}}_{y}$

Figure 11: The deviation $\delta(n)$ for various values of $N$ (integrable case).

Figure 12: The deviation $\delta(n)$ for various values of $N$ (nonintegrable case).

We also study the operator dynamics using a numerical method. We define $\delta(n)=\sum_{i,j}(A_{i,j}A^{*}_{i,j})/2N$, where, $A_{i,j}={\mathcal{U}}^{n}_{i,j}-{\mathcal{U}}_{i,j}$ [74]. Numerically it is observed that the division by $2N$ ensures the average of $\delta(n)$ is one. For any $n$, $\delta(n)=0$ if and only if ${\mathcal{U}}^{n}={\mathcal{U}}$, which confirms the periodic nature of $\mathcal{U}$ implying integrability. The definition above makes sure that $\delta(n)$ is invariant under the global unitary transformation (for proof see Supplemental Material [121]). A similar definition in previous work was lacking this invariance, although it can detect the periodic behavior correctly [74]. For our parameter results are plotted in Fig. 11. We have checked periodicity for $N$ as large as $1000$ and $n$ up to $10000$ (results are not shown here). The observed period is $48$ for all even $N>6$. The same periodic nature could not be observed for (a) even $N$ and $J$ $\notin$ $\left\{1,1/2\right\}$, (b) odd $N$ and $J$ $\neq 1$. Which is shown in the Figs. 11 and 12 (results for odd $N$ are shown in Fig. 12).

Figure 13: Degeneracy of the quasienergies of $\mathcal{U}$ for (a) integrable case and (b) nonintegrable case.

Figure 14: Probability distribution of quasi-energies

Additional signature of integrability can be found from the eigenangle spectrum of $\mathcal{U}$ [94, 80, 98]. We find that it is highly degenerate to take values from the finite set $\{0,\pm\pi/4,\pm\pi/8,\pm 3\pi/4,\pm 3\pi/8,\pm 5\pi/8,\pm 7\pi/8\}$ (see (a) of Fig. 13), whereas for odd $N$, $J=1/2$ and $\tau=\pi/4$, we find that degeneracy disappears (see (b) Fig. 13). In fact, we find that the eigenangles are almost uniformly distributed (see Fig. 14). Thus, with these signatures, we can readily conjecture that the system is quantum integrable for any even-$N>12$, whereas for an odd $N$, it isn’t.

VI Conclusions

In this paper, we have generalized the recently introduced many-body model [74]. The model consisted of qubits with all-to-all ising interaction. Keeping other parameters same, we have generalized this interaction to $J$, with $J=1$ reducing to the recent case [74]. We have studied QI in the generalized version. We have found that apart from $J=1$, QI also exists for $J=1/2$ and only for an even number of qubits. We have analytically calculated the eigensystem, time evolution of the unitary operator, reduced density matrix, and entanglement dynamics for an even number of qubits ranging from $6$ to $12$. We have used linear entropy, von Neumann entropy, and concurrence to measure the entanglement dynamics for various initial unentangled states. The behavior of entanglement dynamics and unitary operator dynamics are in accordance with QI systems with ising interaction. Due to fairly large expressions and complexity, we restricted ourselves to the numerical simulations for any $N>12$ qubits. We have introduced a quantity, $\delta(n)$, to investigate the periodic nature of the time-evolution operator numerically. We have shown that it remains invariant under the global unitary transformation. If this quantity is zero then it indicates the periodic nature of the unitary operator. We have observed that the entanglement dynamics and floquet operator dynamics show periodic behavior for any even $N$, while for any odd $N$, it does not. We have also observed the high degeneracy in the spectra for an even number of qubits. While for odd $N$, the degeneracy in the spectra disappears. We observed that its pairwise concurrence tends to zero with $N$, indicating the multipartite nature of entanglement. Thus, we can very well conjuncture, that the system also exhibits quantum integrability for parameters $J=1/2$ and $\tau=\pi/4$ for any even number of qubits, whereas it is absent for any odd number of qubits.

In recent work, a generalization of the Hubbard-Stratonovich transformation was employed to get an exact analytical solution for quantum-strong-long-range Ising chains [122]. Our model can be explored further in this direction. The experimental verification of our results, (for the smaller number of qubits) can be conducted in various setups such as NMR [123], superconducting qubits [124] and laser-cooled atoms [125], where the QKT has been implemented, whereas for a large number of qubits (of the order of $100$s), one can use ion traps [126, 15]. The validity of our conjecture can be established in this setup. While our study has successfully identified integrability in this model for the specified values of parameters $J=1,1/2$ and $\tau=\pi/4$. We hope our work raises several open questions. Few of them are as follows: 1) As the number of qubits increases, entanglement becomes multipartite in nature. One such measure is the Meyer and Wallach Q measure, which is very much related to linear entropy and we have studied it here. Further analysis using various multipartite entanglement measures can be an interesting direction to pursue [127, 128, 129, 130];  2) Are there other possible values of $J$ that exhibit integrability within this framework?;  3) Our findings encourage further exploration to search for other integrable systems in which all-to-all interaction is present.

VII Acknowledgment

The authors are grateful to the Department of Science and Technology (DST) for their generous financial support, making this research possible through the sanctioned project SR/FST/PSI/2017/5(C) to the Department of Physics of VNIT, Nagpur.