Exact relaxation dynamics and quantum information scrambling in multiply quenched harmonic chains
Abstract
The quantum dynamics of isolated systems under quench condition exhibits a variety of interesting features depending on the integrable/chaotic nature of system. We study the exact dynamics of trivially integrable system of harmonic chains under a multiple quench protocol. Out of time ordered correlator of two Hermitian operators at large time displays scrambling in the thermodynamic limit. In this limit, the entanglement entropy and the central component of momentum distribution both saturate to a steady state value. We also show that reduced density matrix assumes the diagonal form long after multiple quenches for large system size. These exact results involving infinite dimensional Hilbert space are indicative of local thermal behaviour for a trivially integrable harmonic chain.
I Introduction
The behaviour of isolated quantum systems under nonequilibrium conditions is of great interest from both theoretical Rigol et al. (2007); Marcos et al. (2008); Polkovnikov et al. (2011); Srednicki (1994); D’Alessio et al. (2016); Deutsch (2018) as well as experimental point of view Islam et al. (2015); Vidmar et al. (2015); Kaufman et al. (2016). A generic method to reach a nonequilibrium regime is by a single or multiple quenches of the system parameters. For nonintegrable systems, the quench protocol often leads to thermalization Srednicki (1994); Marcos et al. (2008), although there can be subtleties depending on the choice of initial states Barthel and Schollwöck (2008); Cramer et al. (2008a).
For integrable systems, the situation is not as unequivocal. It is generally known that an isolated integrable system does not thermalize Marcos et al. (2008); Deutsch (2018), but can be described by a generalized Gibbs ensemble (GGE) Rigol et al. (2007); Vidmar and Rigol (2016); Essler and Fagotti (2016); Fioretto and Mussardo (2010).However, it has recently been found that the Heisenberg spin chain, integrable under the Bethe ansatz scheme, can exhibit weak eigenstate thermalization for a typical, but not necessarily all eigenstates Alba (2015); Alba and Calabrese (2017); Grisins and Mazets (2011a). It is also known that certain finite dimensional systems exhibit comparable statistical relaxation regardless of whether they are integrable or not Santos et al. (2012). Furthermore, an integrable JaynesCummings model can show rapid decoherence analogous to chaotic dynamics and the inability to recover the purity at large times increases in the thermodynamic limit Angelo et al. (1999). For an integrable system with infinite dimensional Hilbert space, Bogoliubov and Krylov showed that under certain assumption about the thermal reservior, the system relaxes to thermal state in thermodynamic limit (see Chirikov (1986) and references therein). In 1dbosonic systems, a class of special initial states thermalizes under integrable dynamics governed by GrossPitaveskii equation Grisins and Mazets (2011b).
Various physical quantities have been used to analyze the onset of statistical relaxation and thermalization following a quench. In a quantum manybody system, the reduced density matrix (RDM) carries important information regarding the relaxation dynamics Peschel and Chung (1999); Chung and Peschel (2000); Peschel (2004); Peschel and Eisler (2009); Peschel (2003); Fagotti and Essler (2013); Calabrese et al. (2012). In particular, the vanishing of the offdiagonal matrix elements of the reduced density matrix indicates the emergence of a thermal state. Similarly, the entanglement entropy calculated from the RDM measures the loss of information and its dynamics describes how the quantum information is spread leading to local thermal behaviour Islam et al. (2015); Kaufman et al. (2016). The one body momentum distribution, also obtained from the RDM, carries the signature of the thermal behavior of the system Marcos et al. (2008).
Recently, the out of time order correlator (OTOC) Larkin and Ovchinnikov (1969) has gained prominence in the context of scrambling of quantum information in nonequilibrium systems Maldacena et al. (2016); Rozenbaum et al. (2017); Swingle (2018); Heyl et al. (2018); Lin and Motrunich (2018a, b); Li et al. (2017); Kukuljan et al. (2017); Luitz and Bar Lev (2017); Cotler et al. (2018). Although information scrambling is usually a property of chaotic systems, the OTOC of certain nonlocal operators exhibit scrambling even in an integrable Ising chain Lin and Motrunich (2018a). It has been further argued that scrambling could be independent of the integrability of the Hamiltonian Iyoda and Sagawa (2018) and that mutual information in an integrable spin chain can exhibit weak scrambling Alba and Calabrese (2019). The above list of examples, which is by no means exhaustive, indicates that there are still many open questions in the context of thermalization and information scrambling in integrable systems.
One of the key points of this work is to address the question of dynamical relaxation and information scrambling in an integrable system. The main differences of our analysis with the existing approaches in the literature are

We consider a series of multiple global quenches with no restriction to their number, in contrast to a single quench which is typical in the existing literature. As we shall show, the multiple quench protocol leads to important differences compared to a single quench.

Our analysis involves systems with infinite dimensional Hilbert space at each site and local Hermetian operators. This is particularly relevant for information scrambling and OTOC where majority of the exisiting analyses involves finite spin systems and unitary operators.

The time development is treated in an exact analytical fashion by solving TDSE using nonlinear Ermakov equation. The total number of particles and the number of quenches can be arbitrarily large and the analysis is valid for indefinitely large time.
In order to achieve the above, we study the exact relaxation dynamics following multiple global quenches in a harmonic chain with oscillators which is trivially integrable. For this system, the exact entanglement dynamics following a single quench has been studied earlierGhosh et al. (2017). Under a multiple quench protocol, we show that the offdiagonal terms in the oneparticle reduced density matrix Peschel and Chung (1999); Peschel (2003) vanish exponentially fast in time in the limit of large , while the diagonal terms saturates to a steady value. The mixing of a large number of incommensurate and irrational normal mode frequencies is responsible for this feature. In contrast, for a single quench, the relaxation to the steady state is much slower. This feature indicates that the system exhibits locally thermal behaviour under multiple quenches. Under the same conditions, the momentum distribution and the entanglement entropy also show the signatures of a steady state.
Universality of the small time behaviour of OTOCs in chaotic systems and its relation with Lyapunov exponents has been the focus of most of the studiesRozenbaum et al. (2017); Maldacena et al. (2016); Swingle (2018); Lakshminarayan and Prakash (2019). Here we want to focus on long time behaviour much after the completion of the quench protocol. The multiple quench protocol has important consequences in the quantum information scrambling and OTOC. While recurrences are characteristic feature of integrability in OTOC, in contrast for this system OTOC saturates to a nonzero value under multiple quenches in thermodynamic limit, which is indicative of scrambling Alba and Calabrese (2019). As our exact analysis is valid for arbitrarily large number of particles, we can clearly demonstrate the finite size effect by varying the number of oscillators.
The harmonic chain under consideration can be experimentally realized in various systems such as in optical tweezers Barredo et al. (2016) and the tuning of individual coupling parameters can be done using ultracold atoms Hunger et al. (2011) or Rydberg states Buchmann et al. (2017); Macrì and Pohl (2014). In particular, it can be simulated using BoseHubbard model in the strong superfluid phaseCramer et al. (2008b). Various quench protocols have already been experimentally realized for the BoseHubbard model using cold atoms in optical lattices Vidmar et al. (2015); Kaufman et al. (2016); Islam et al. (2015). It is thus plausible that the predictions of our work can be empirically verified.
This paper is organized as follow. In Sec. 2, we setup the formalism for multiple quench protocol and obtain the solutions of TDSE using nonlinear Ermakov equations. In Sec. 3, we discuss the time dependence of RDM and show that the offdiagonal elements vanish for large time and in the thermodynamic limit. In Sec. 4, we obtain the momentum distribution and enetanglement entropy using the RDM and discuss their steady state properties. Quantum information scrambling is discussed in Sec. 5 and it is shown that OTOC saturates to nonzero steady value for the case of multiple quenches in the same limit. We conclude the paper in Sec. 6 with a summary of results and remarks.
Ii Harmonic chain and the quench protocol
We consider an isolated harmonic oscillator chain with oscillators and with periodic boundary condition. The Hamiltonain is given by
(1)  
where = and is an real symmetric matrix. Here the frequency and the nearest neighbour coupling are explicit functions of time. A solution of the corresponding time dependent Schrödinger’s equation (TDSE) Lewis and Riesenfeld (1969); Lohe (2009); Ghosh et al. (2017)
(2) 
can be written as
(3)  
where , , ,, is an orthogonal transformation which transforms the matrix to its diagonalized form and is a diagonal matrix with elements . Here satisfies the the nonlinear Ermakov equation Lewis and Riesenfeld (1969); Lohe (2009); Pinney (1950) given by
(4) 
where ’s are the normal mode frequencies of the Hamiltonian, which have the form
(5) 
with . We choose the initial condition as the ground state of the time independent prequenched Hamiltonian of the oscillator chain, which requires that and . Note that finding the solutions of the TDSE is equivalent to finding the solutions of the Ermakov equations (4).
Solution for series of quenches: Our quench protocol is shown in Fig. (1). At time , the frequency is quenched to and is left unchanged. This defines a single quench. After a time , the frequency is changed back to , defining the second quench. This sequence is repeated till the required number of quenches is achieved. Now we shall discuss the solution of Ermakov equations under such quenches. In terms of the variable
(6) 
the Ermakov equations (4) can be written as
(7) 
The corresponding solutions are given by
(8) 
where , and are the time independent constants which are obtained from the boundary conditions on or . Furthermore, using the above equation we could express and its derivatives in a matrix form as
(9) 
where
(10)  
The continuity of wavefunction across the quenches implies the continuity in and . and defined just after and before any quench respectively are related by
(11) 
where the form of B matrix could be found using Eq.(7) as
(12) 
The coefficient vectors with elements across the quench are related as
(13) 
with the boundary condition,
(14) 
This gives a complete solution of the Ermakov equations for the given quench protocol.
Iii Reduced density matrix
In this Section we shall discuss the timedependent onebody RDM obtained under a series of quenches in the coordinate representation. This RDM would subsequently be used to obtain the onebody momentum distribution and the entanglement entropy.
In order to obtain the time dependent onebody RDM, we trace out the rest of the system except the site. The RDM has the form
(15) 
Using (3) in (15) and carrying out the integration, we get
(16)  
where , , are given by
(17)  
A generic offdiagonal element of the RDM has the form
(18)  
where denotes the distance from the diagonal axis of the RDM and corresponds to the diagonal elements of the RDM. The equation (18) shows that the off diagonal elements are shifted Gaussians with a time dependent exponent. The ratio between the diagonal and a generic offdiagonal element can be written as
(19) 
The exponents and are plotted in Fig. 2 and are found to be real positive for all time. This also follows from the reality of the eigenvalues of the RDM, which is Hermitian Ghosh et al. (2017). For a series of multiple quenches we find that rapidly saturates to a small positive value irrespective of . However has very different behaviour depending on the value of . For higher , has a large value which oscillates very little with time. For smaller , the mean value of is also high but the fluctuations are very big as well. On the other hand, for a single quench, both and show appreciable oscillations and the value of is much less irrespective of the value of .
In Eq.(III) we have given the ratio of the offdiagonal to the diagonal elements of the RDM. We thus find that for multiple quenches and in the thermodynamic limit, the offdiagonal matrix elements of the RDM tend to zero. For a single quench, the suppression of the offdiagonal matrix elements with time is much less.
Hence, for multiple quenches with the large system size, the RDM assumes a diagonal form at a time large compared to the duration of the quenches. This can be qualitatively understood as follows. The solution of the Ermakov equation for each normal mode contains an irrational number given by the square root of corresponding normal mode frequency . As the number of oscillators increases, a large number of irrational and incommensurate frequencies start contributing to the wave function and the RDM. The mixing of a large number of such modes is responsible for the statistical relaxation of the quantities and with time, which in turn reduces the RDM to a diagonal form.
The diagonal form of the RDM is generally indicative of a thermal behaviour. For a quenched integrable harmonic chain, the thermal behaviour is governed by the generalized Gibbs ensemble Biroli et al. (2010) defined in the number representation. In our analysis, the RDM has been obtained in the coordinate representation. Moreover, in order to establish the GGE for a bipartite system, it is necessray to first keep the subsystem size fixed and make the complement infintely large and then to take a thermodynamic limit for the subsystem size as well Calabrese et al. (2012). In our analysis, the subsystem size has been kept fixed while only the size of the complement has been increased. Thus in the limit of subsystem size becoming large, it is expected that reduced density matrix will be described by GGE.
Iv momentum distribution and entropy
We derive the analytical expression of one body momentum distribution by taking the Fourier transform of (16), which is given by
(20)  
After the integration, the momentum distribution takes the form
(21) 
At time just before the quench the distribution (21) is a gaussian with mean at . In Fig. 3 the central component of the momentum distribution is plotted with time for both single and multiple quenches and with . For a single quench, in the long time limit the momentum distributions reaches a steady value but still shows appreciable oscillations. For multiple quenches the value in the long time limit is lower and the fluctuations are negligible. Thus in the long time and in the thermodynamic limit, the higher the number of quenches, the better is the relaxation of the system to a steady state.
In order to study the entanglement entropy, we bipartite the system into two parts of one oscillator versus oscillators. The entanglement entropy of the smaller subsystem calculated using the reduced density matrix (16) has the form Ghosh et al. (2017)
(22) 
where the has the form given by
(23) 
, are given in eq. (17).
The von Neumann entropy as a function of time is plotted in Fig. 4. In a finite dimensional system, the maximum entropy is proportional to the Hilbert space dimension. The bound on entanglement entropy is saturated for generic chaotic systems. Here the system being infinite dimensional, there is no finite limit to the maximum entanglement entropy. As in the case of the momentum distribution, higher number of quenches leads to smaller fluctuations in the entanglement entropy and a larger steady state value. This can be qualitatively understood in terms of reduced variation in eigenvalues of RDM given by Ghosh et al. (2017). The variation in reduces as or in other words as (see Fig. 2).
V Quantum information scrambling
The OTOC between two operators and separated by a lattice distance is defined as
(24) 
where . For well separated local operators, the OTOC starts from a zero value and then increases as the information propagates with LiebRobinson velocity.
Here we consider the OTOC between position and momentum operators, which are local and Hermitian. We choose , which are labelled by the site index. The expression for in the Heisenberg picture for a single time dependent quantum harmonic oscillator with frequency is given by Ji et al. (1995)
(25) 
where is the solution of corresponding Ermakov equation. Using Eq.(V) with the canonical commutation relations we get,
(26) 
In Fig.(5), OTOC is plotted for a fixed system size and with different number of quenches. We observe that OTOC saturates to a nonzero value with very small fluctuation for the case of multiple quenches in contrast to a constant slow decrease with large fluctuation for single quench. A strong system size dependence on the OTOC is clearly seen from Fig. 6. We find that even with multiple quenches, the OTOC fluctuates enormously for smaller number of particles and does not saturate to any steady value. In Fig.(5), the quasirecurrence in the case of multiple quenches for occurs at a time where is the LiebRobinson velocity for this system with the given set of parameters after quench. We also note that the quasirecurrences decrease in the thermodynamic limit as well as with the higher number of quenches.
The saturation of the OTOC to a nonzero steady value is indicative of information scrambling. Even though the system under consideration is integrable, we find that under multiple quenches and in the thermodynamic limit, OTOC saturates to a nonzero steady value. In this formalism, the saturation can be attributed to the mixing of large number of modes with incommensurate frequencies. In the conventional picture of expanding the state as superposition of eigenfunctions of postquench Hamiltonian, the multiple quenches would amount to accessing larger proportions of Hilbert space. This is also consistent with larger entanglement entropy with increasing number of quenches.
Vi conclusion
In this paper we have analyzed the relaxation dynamics and quantum information scrambling in an isolated harmonic chain under multiple quenches. The various physical quantities show remarkably different nonequilibrium behaviour under the multiple quench protocol compared to a single quench. The RDM has been shown to assume a diagonal form exponentially fast compared to a single quench. The entanglement dynamics and the momentum distribution also shows relaxation to a steady state. The exact analytical results obtained here are valid for arbitrary number of particles. In order to demonstrate the finite size effects, we have graphically exhibited our results for and . It is clearly seen that the quasirevivals of the physical observables, characteristic of the finite size effects, reduce remarkably as is increased.
The conclusions are similar for the quantum information scrambling and OTOC. For a single quench and for a low value of , the OTOC shows almost complete revival and the scrambling is practically nonexistent. On the other hand, for five quenches and with , the OTOC saturates to a finite steady value with very little revival. It may be noted that we have considered Hermitian operators to evaluate the OTOC and the system has an infinite dimensional Hilbert space. This is very different compared to the usual finite spin systems where the OTOC is normally evaluated using unitary operators.
The above result indicates that for multiple quenches and in the thermodynamic limit, the integrable harmonic chain relaxes to a steady state, the RDM assumes a thermal form and the OTOC exhibits quantum information scrambling. In our formalism this has been achieved using an exact solution of the TDSE valid throughout the quench protocol and all the observables have been evaluated using the exact solution. As the number of the oscillators is increased, a large number of irrational and incommensurate normal mode frequencies start contributing to the observables whose number is of the order of . The mixing of these incommensurate frequncies leads to the steady state and saturating behaviour of the various physical quantities.
It is customary to describe the thermal behaviour of integrable systems in terms of GGE. As discussed in Section III, for a fixed subsystem size, it is difficult to make a direct comparison of our results with the GGE. A larger question is whether the system begins to loose integrability under the multiple quench protocol. For a smooth time dependence of the system parameters, a single particle harmonic oscillator remains integrable Angelo et al. (2012). Here we consider sudden quenches which are certainly not smooth in time. Hence, the status of the integrabiliy under the considered quench protocol is not quite clear.
References
 Rigol et al. (2007) M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, Phys. Rev. Lett. 98, 050405 (2007).
 Marcos et al. (2008) R. Marcos, D. Vanja, and O. Maxim, Nature 452, 854 (2008).
 Polkovnikov et al. (2011) A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Rev. Mod. Phys. 83, 863 (2011).
 Srednicki (1994) M. Srednicki, Phys. Rev. E 50, 888 (1994).
 D’Alessio et al. (2016) L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, Advances in Physics 65, 239 (2016), https://doi.org/10.1080/00018732.2016.1198134 .
 Deutsch (2018) J. M. Deutsch, Reports on Progress in Physics 81, 082001 (2018).
 Islam et al. (2015) R. Islam, R. Ma, P. M. Preiss, M. Eric Tai, A. Lukin, M. Rispoli, and M. Greiner, Nature 528, 77 (2015).
 Vidmar et al. (2015) L. Vidmar, J. P. Ronzheimer, M. Schreiber, S. Braun, S. S. Hodgman, S. Langer, F. HeidrichMeisner, I. Bloch, and U. Schneider, Phys. Rev. Lett. 115, 175301 (2015).
 Kaufman et al. (2016) A. M. Kaufman, M. E. Tai, A. Lukin, M. Rispoli, R. Schittko, P. M. Preiss, and M. Greiner, Science 353, 794 (2016), https://science.sciencemag.org/content/353/6301/794.full.pdf .
 Barthel and Schollwöck (2008) T. Barthel and U. Schollwöck, Phys. Rev. Lett. 100, 100601 (2008).
 Cramer et al. (2008a) M. Cramer, A. Flesch, I. P. McCulloch, U. Schollwöck, and J. Eisert, Phys. Rev. Lett. 101, 063001 (2008a).
 Vidmar and Rigol (2016) L. Vidmar and M. Rigol, Journal of Statistical Mechanics: Theory and Experiment 2016, 064007 (2016).
 Essler and Fagotti (2016) F. H. L. Essler and M. Fagotti, Journal of Statistical Mechanics: Theory and Experiment 2016, 064002 (2016).
 Fioretto and Mussardo (2010) D. Fioretto and G. Mussardo, New Journal of Physics 12, 055015 (2010).
 Alba (2015) V. Alba, Phys. Rev. B 91, 155123 (2015).
 Alba and Calabrese (2017) V. Alba and P. Calabrese, Proceedings of the National Academy of Sciences 114, 7947 (2017), https://www.pnas.org/content/114/30/7947.full.pdf .
 Grisins and Mazets (2011a) P. Grisins and I. E. Mazets, Phys. Rev. A 84, 053635 (2011a).
 Santos et al. (2012) L. F. Santos, F. Borgonovi, and F. M. Izrailev, Phys. Rev. Lett. 108, 094102 (2012).
 Angelo et al. (1999) R. M. Angelo, K. Furuya, M. C. Nemes, and G. Q. Pellegrino, Phys. Rev. E 60, 5407 (1999).
 Chirikov (1986) B. V. Chirikov, Foundations of Physics 16, 39 (1986).
 Grisins and Mazets (2011b) P. Grisins and I. E. Mazets, Phys. Rev. A 84, 053635 (2011b).
 Peschel and Chung (1999) I. Peschel and M.C. Chung, Journal of Physics A: Mathematical and General 32, 8419 (1999).
 Chung and Peschel (2000) M. Chung and I. Peschel, Phys.Rev. B 62, 4191 (2000).
 Peschel (2004) I. Peschel, Journal of Statistical Mechanics: Theory and Experiment 2004, P06004 (2004).
 Peschel and Eisler (2009) I. Peschel and V. Eisler, Journal of Physics A: Mathematical and Theoretical 42, 504003 (2009).
 Peschel (2003) I. Peschel, Journal of Physics A: Mathematical and General 36, L205 (2003).
 Fagotti and Essler (2013) M. Fagotti and F. H. L. Essler, Phys. Rev. B 87, 245107 (2013).
 Calabrese et al. (2012) P. Calabrese, F. H. L. Essler, and M. Fagotti, Journal of Statistical Mechanics: Theory and Experiment 2012, P07016 (2012).
 Larkin and Ovchinnikov (1969) A. Larkin and Y. N. Ovchinnikov, Journal of Experimental and Theoretical Physics 28, 1200 (1969).
 Maldacena et al. (2016) J. Maldacena, S. H. Shenker, and D. Stanford, Journal of High Energy Physics 2016, 106 (2016).
 Rozenbaum et al. (2017) E. B. Rozenbaum, S. Ganeshan, and V. Galitski, Phys. Rev. Lett. 118, 086801 (2017).
 Swingle (2018) B. Swingle, Nature Physics 14, 988 (2018).
 Heyl et al. (2018) M. Heyl, F. Pollmann, and B. Dóra, Phys. Rev. Lett. 121, 016801 (2018).
 Lin and Motrunich (2018a) C.J. Lin and O. I. Motrunich, Phys. Rev. B 97, 144304 (2018a).
 Lin and Motrunich (2018b) C.J. Lin and O. I. Motrunich, Phys. Rev. B 98, 134305 (2018b).
 Li et al. (2017) J. Li, R. Fan, H. Wang, B. Ye, B. Zeng, H. Zhai, X. Peng, and J. Du, Phys. Rev. X 7, 031011 (2017).
 Kukuljan et al. (2017) I. Kukuljan, S. c. v. Grozdanov, and T. c. v. Prosen, Phys. Rev. B 96, 060301 (2017).
 Luitz and Bar Lev (2017) D. J. Luitz and Y. Bar Lev, Phys. Rev. B 96, 020406 (2017).
 Cotler et al. (2018) J. S. Cotler, D. Ding, and G. R. Penington, Annals of Physics 396, 318 (2018).
 Iyoda and Sagawa (2018) E. Iyoda and T. Sagawa, Phys. Rev. A 97, 042330 (2018).
 Alba and Calabrese (2019) V. Alba and P. Calabrese, arXiv:1903.09176 (2019), https://arxiv.org/abs/1903.09176 .
 Ghosh et al. (2017) S. Ghosh, K. S. Gupta, and S. C. L. Srivastava, EPL (Europhysics Letters) 120, 50005 (2017).
 Lakshminarayan and Prakash (2019) A. Lakshminarayan and R. Prakash, arXiv:1904.06482 (2019), https://arxiv.org/abs/1904.06482 .
 Barredo et al. (2016) D. Barredo, S. de Léséleuc, V. Lienhard, T. Lahaye, and A. Browaeys, Science 354, 1021 (2016).
 Hunger et al. (2011) D. Hunger, S. Camerer, M. Korppi, A. Jöckel, T. Hänsch, and P. Treutlein, C. R. Physique 12, 871 (2011).
 Buchmann et al. (2017) L. F. Buchmann, K. Mølmer, and D. Petrosyan, Phys. Rev. A 95, 013403 (2017).
 Macrì and Pohl (2014) T. Macrì and T. Pohl, Phys. Rev. A 89, 011402 (2014).
 Cramer et al. (2008b) M. Cramer, C. M. Dawson, J. Eisert, and T. J. Osborne, Phys. Rev. Lett. 100, 030602 (2008b).
 Lewis and Riesenfeld (1969) H. R. Lewis and W. B. Riesenfeld, J. Math. Phys. 10, 458 (1969).
 Lohe (2009) M. A. Lohe, J. Phys. A: Math. Theor. 42, 035307 (2009).
 Pinney (1950) E. Pinney, Proc. Amer. Math. Soc. 1, 681 (1950).
 Biroli et al. (2010) G. Biroli, C. Kollath, and A. M. Läuchli, Phys. Rev. Lett. 105, 250401 (2010).
 Ji et al. (1995) J.Y. Ji, J. K. Kim, and S. P. Kim, Phys. Rev. A 51, 4268 (1995).
 Angelo et al. (2012) R. M. Angelo, E. I. Duzzioni, and A. D. Ribeiro, Journal of Physics A: Mathematical and Theoretical 45, 055101 (2012).