Dynamics of Entanglement and `Attractor' states in The Tavis-Cummings Model

We study the time evolution of $N_q$ two-level atoms (or qubits) interacting with a single mode of the quantised radiation field. In the case of two qubits, we show that for a set of initial conditions the reduced density matrix of the atomic system approaches that of a pure state at $\sfrac{t_r}{4}$, halfway between that start of the collapse and the first mini revival peak, where $t_r$ is the time of the main revival. The pure state approached is the same for a set of initial conditions and is thus termed an `attractor state'. The set itself is termed the basin of attraction and the features are at the center of our attention. Extending to more qubits, we find that attractors are a generic feature of the multi qubit Jaynes Cummings model (JCM) and we therefore generalise the discovery by Gea-Banacloche for the one qubit case. We give the `basin of attraction' for $N_q$ qubits and discuss the implications of the `attractor' state in terms of the dynamics of $N_q$-body entanglement. We observe both collapse and revival and sudden birth/death of entanglement depending on the initial conditions.

the phenomena of 'collapse and revival' of the Rabi oscillations in the qubit system when the field is started in a coherent state |α [7], Here the state |n is the eigenstate of the photon number operatorn = a † a with eigenvalue n and θ is the initial phase of the radiation field. |α| determines the size of the field, withn = α|n|α = |α| 2 .
In short what is found is that the initial Rabi-oscillations of the probability of being in a given qubit state decay on a time scale called the collapse time, t c , but then revive after a much longer time, t r , due to the quantum 'graininess' of the radiation field [8,9,10]. Whilst the nature of this remarkable phenomenon in the case of one qubit is now well understood [6], both experimentally [2] and theoretically [11], the multi qubit case [5,12] has not been explored to the same extent.
For instance, the occurrence between t c and t r , of what we shall call the 'attractor' state has only been fully investigated in the one qubit limit. This intriguing aspect of quantum dynamics was discovered and highlighted by Gea-Banacloche [13,14] in his study of the one qubit case. In this paper we report on our investigation of analogous 'attractors' in cases of more than one qubit.
The organization of this paper is as follows. In Sec. II, the dynamics of the resonant JCM (single qubit coupled to a cavity system) is presented. This sets up our notation, analytic tools and methods of data presentation. In Sec. III we extend the model to two qubits and show that the idea of an 'attractor' state still exists. We investigate the dynamics in depth and quantify the set of initial states that approach the 'attractor', giving the 'basin of attraction'. We show that the atomic system exhibits interesting entanglement properties, in particular that the entanglement can vanish and then return at a later time. Furthermore, the manner in which this occurs depends on whether we are in or outside this basin of attraction. In Sec. IV we develop these ideas further from those of the two-qubit case and find a set of initial states (the basin of attraction) from which the system always goes to the 'attractor' state for an arbitrary number of qubits. In the limit of large qubit number we rewrite our qubit equations in terms of spin coherent states and demonstrate some intriguing phenomena associated with the N q -qubit JCM, especially the collapse and revival of N q -body entanglement. Whilst our method of analysis, and some of our results, are similar to those of Chumakov et al. [15] and Meunier et al. [16], our focus is on the quantum information aspects of the problem and in particular the dynamics of entanglement.

II. FEATURES OF THE JAYNES-CUMMINGS MODEL
We first review the JCM and appearance of 'attractor' states in the one qubit case [17,18,19].
The JCM model consists of a single qubit coupled to a single mode of a quantised electromagnetic field. The two possible states of the qubit are the ground state |g , with energy ǫ g , or its excited state |e , with energy ǫ e . We consider the following Hamiltonian: σ z = |e e| − |g g|,σ + = |e g|,σ − = |g e| whereâ † andâ are the creation and annihilation operators of photons with frequency ω, Ω = ǫ e −ǫ g and λ the cavity-qubit coupling constant. The model is only valid close to resonance and we shall only consider the resonant case of ω = Ω.
We have an initial state: where the state of the qubit is normalised so |C e | 2 + |C g | 2 = 1 and the field is started in a coherent state. In the interaction picture (i.e. in a rotating reference frame) we obtain a solution to the Hamiltonian: C e C n cos λ √ n + 1t − iC g C n+1 sin λ √ n + 1t |e, n C g C n+1 cos λ √ n + 1t − iC e C n sin λ √ n + 1t |g, n + 1 + C g C 0 |g, 0 .
In spite of its relative simplicity, Eq. (5) exhibits a wide range of interesting phenomena. A much studied example is the collapse and revival of Rabi oscillations. If the qubit is started in the ground state then initially the probability of the qubit being in the ground state oscillates from one to zero at the Rabi frequency. These oscillations then decay on a time scale called the collapse time, t c = 2 λ , but then revive after a much longer time: Whenn ≫ λt we can approximate |Ψ 1 (t) in Eq. (5) by Ψ 1 (t) : where β ± 1 2 (t) is a normalisation factor, D ± 1 2 (t) is a state of the qubits and Φ ± 1 2 (t) is a state of the field.
With the solution written in this form it is clear there are two distinct parts to the wavefunction.
In fact the field just consists of two different coherent states that are out of phase by λt/ √n . In a penetrating study of this model Gea-Banacloche [14] noted that at a time tr 2 , in the very largen approximation, the qubit disentangles itself completely from the field and the wavefunction factors into a product state of the radiation field and the qubit.
To make this more apparent we rewrite Eq. (8) -Eq. (10) at the time t = tr 2 .
With the wavefunction written in this way it can be seen that the two states D + 1 2 ( tr 2 ) and D − 1 2 ( tr 2 ) are in fact the same at the time tr 2 . We label this state ψ + 1 att : This state is independent of the initial coefficients C e and C g and only depends on the phase θ of the initial coherent state, Eq. (1). Throughout this paper we follow Phoenix and Knight [20] and call this state an 'attractor' state. Note that we use the term 'attractor' to refer to ψ ± 1 att even though in the theory of non-linear differential equations the term attractor is reserved for solutions which, when reached, after evolving from an arbitrary initial condition, stop changing with time. As a consequence the qubit state can be factorised out of the wavefunction meaning the qubit and the field are in a product state: It is clear that all the information of the initial state about the qubit is stored in the radiation field which happens to be in a "Schrödinger cat" state. At a later time t = 3tr 2 , the qubit again disentangles itself completely from the field and goes to a pure state. The qubit state at this time is orthogonal to ψ + 1 att and we label it ψ − 1 att . This remarkable behaviour is most unlike the consequences of a linear Schrödinger equation satisfied by a qubit without the cavity field or with a classical field and, as was stressed by Gea-Banacloche, is a natural route to quantum state preparation. All initial qubit states tend to the same 'attractor' state (which is determined by the initial phase of the field state) and the quantum information about the initial qubit state has effectively been swapped out and encoded into the field at the 'attractor times'.
Before we move on to describe the results of numerical simulations we consider the field, Eq.
(10), at two revival times t = t r and t = 2t r . At both these times the field states Φ ± 1 As a consequence the system is once again in a product state as the field can be factorised out of the wavefunction for the whole state. The information about the initial state of the qubits has been transferred from the field back to the qubits.
By considering the average number of photons in the field,n, to be large many interesting observations have been made of the JCM. It was shown that the system repeatedly evolves from a product state to an entangled state and back again. Although the analytical equations for large values ofn clearly show this behaviour, we now consider the exact solution |Ψ 1 (t) for modestn (n = 50), shown in Eq. (5) and comment on the results of numerical simulations for the qubit and the photon sector separately.
We consider the state of the qubits after tracing out the field state. The probability of the qubits being in the state |g is, where ρ q (t) = Tr f (|Ψ 1 (t) Ψ 1 (t)|) is the reduced density matrix at a certain time for the qubits when the field has been traced out. Throughout this paper all instances of the letter q refer to the qubit system and f refers to the field. The results are depicted in Fig. 1 for the initial state |Ψ 1 (t = 0) = |g |α and the well known phenomenon of collapse and revival can be seen. For numerical reasons the sum over n is truncated at n max = 200 as this is adequate for the average number of photons in the mode ofn = 50.
An entropy measurement of the system can be made at any point in time. Entropy is a well defined quantity that measures the disorder of a system and the purity of a quantum state. In cases where the time-dependent Schrödinger equation (as opposed to some master equation for reduced dynamics) determines the time evolution of a system, the value of the entropy remains constant. The complete systems (qubits plus field) we study are always started in a pure state and follow Schrödinger dynamics; therefore the entropy of the complete systems always remains zero as we have a pure state. In these systems the individual entropies, of the qubits or the field, are of more interest. We calculate the partial entropy from the reduced density matrix ρ q using the rescaled von Neumann entropy [21,22], which is rescaled such that the entropy of the qubits ranges from zero, for a pure state, to one, for a maximally mixed state. For the bipartite case (the one qubit N q = 1 JCM) this reduces to the standard von Neumann entropy. As was discovered by Araki and Lieb [23], if the system is started in a pure state then S q (t)=S f (t) for all time, where S f (t) is the entropy of the field system calculated from reduced density matrix ρ f (t) = Tr q (|Ψ 1 (t) Ψ 1 (t)|). Therefore for the rest of the paper we shall refer to the value of S q (t) only when we mention entropy.
In quantum information these partial entropies provide a measure of the entanglement present between the two systems, known as the entropy of entanglement [24,25]. If the entropy is zero then the state being measured is pure and there is no entanglement present. If the value of the entropy is non zero then there is some entanglement present between the two parts of the system and the reduced density matrix is mixed. Returning to the JCM we note that the 'attractor' state, Eq. (14), is a pure state and hence the value S q ( tr 2 ) will be zero. So the qubit has disentangled itself from the field and there is no entanglement present at this time. The reduced entropy is a powerful tool in understanding entanglement in the JCM and it is worthwhile recalling an example of its power to elucidate the quantum dynamics [20,26].
We have numerically evaluated S q for qubits started in the initial state |Ψ 1 (t = 0) = |g |α using the exact solution for |Ψ 1 (t) and the results are shown in Fig. 1 by the curve labelled (a).
We also plot the probability P 1 att (t) = ψ + 1 att ρ q (t) ψ + 1 att that the qubit is in the 'attractor' state ψ + 1 att , as defined in Eq. (14), shown by the curve labelled (c). Evidently at the time tr 2 the reduced entropy goes very close to zero and P 1 att ( tr 2 ) goes very close to one. As mentioned above, this means there is very little field-qubit entanglement present and the qubit is essentially in the 'attractor' state, ψ + 1 att , as defined in Eq. (14). At 3tr 2 the probability of being in this 'attractor' state approaches zero, as the qubit is in the state ψ − 1 att , which is orthogonal to ψ + 1 att . These results all agree with the discussions of Eq. (15). At future 'attractor' times the system alternates between ψ + 1 att and ψ − 1 att . A big difference between numerical results and the largen solution is noticed at the time t = t r .
Earlier we showed that the solution in the largen approximation will return to a pure state at this time, but in Fig. 1 the entropy has a large value of 0.7, wedged between two maxima (indicating large entanglement between qubit and field). This is because the width of this dip is on a time scale similar to the collapse time, much narrower than the entropy dip at tr 2 . As the value ofn is increased in numerical simulations, the value of the entropy at t = t r gets closer to zero, but we see that even for fairly large coherent states, the largen approximation can deviate significantly from the full solution.
So far we have reviewed the time evolution of the reduced qubit density matrix in the one-qubit and one mode system. Before we move on to address the two-qubit case we calculate the Q function [27] to investigate the radiation field, This is shown in Fig. 2 in the rotating reference frame for the exact solution as contour plots.
We give an example of a 3d plot of Fig. 2(b) in Fig. 3, but for the rest of the paper we will use the 2D contour plots.
The Q function shows a set of two wavepackets of the cavity field, each of which can be thought of as a Gea-Banacloche state [14] in the largen solution. As a function of time the two 'blobs' move around the complex plane and follow the circular path of radius √n . These 'blobs' start to smear out and become distorted from the original coherent state. Although the wavepackets all begin in the same place, they evolve with different frequencies, so the states begin to separate.
After a period of time, depending on the differences between the frequencies, the different 'blobs' a : t 0  are separated by more than their diameter and can be easily distinguished. Eventually these two wavepackets will again overlap and this occurs at the revival time, t = t r . When these 'blobs' have some overlap in phase space there will be Rabi oscillations in the qubit probability, giving the periodic sequence of revivals seen in Fig. 1. However, as can be seen in Fig. 2(d), at the revival time, t r , the 'blobs' have undergone some distortion and so the wave packets are no longer simple coherent states. This distortion elongates the wave packets, and so they overlap for a longer time than in the ideal coherent state case. We therefore understand the lack of a full amplitude revival, P 1g (t r ) = 1, as due to the spreading, and hence distinguishability, of the field states. When the value ofn gets very large there is less distortion and the Gea-Banacloche states stay exact coherent states throughout their evolution.
Also included in Fig. 2 is an arrow for each 'blob' which represents the average atomic po- The N q qubit Hamiltonian, which in this section we treat for N q = 2 so i = 1, 2, is a generalisation of Eq. (2). It is called the Tavis-Cummings Hamiltonian [5]: For brevity we only consider the solution for the case λ 1 = λ 2 = λ and ω = Ω 1 = Ω 2 . This model has been studied in detail by Chumakov et al. [28] and the largen solution has been found [15,16] with similar results to our analysis. Adding another qubit to the system yields a more complex system which can exhibit interesting entanglement properties as not only can the two qubits be entangled with the field, but they can also be entangled between themselves.
We have investigated the evolution of the system over time, with analytical solutions in the largen approximation [29] and the numerically evaluated exact solution. To summarise these we start with the initial state: where the state is normalised so |C ee | 2 + |C eg | 2 + |C ge | 2 + |C gg | 2 = 1. Following the same procedure as in the one qubit case we solve the Hamiltonian and approximate for largen. The time evolution of this state takes the following form: where k = −1, 0 or 1 and: |Φ k (t) = e 2ikπt/tr α .
Note that the state |D 0 (t) = |D 0 is time-independent and is thus a maximally entangled qubit state for all times if C eg = C ge [29]. The field states are called Gea-Banacloche states and once again have a specific qubit state associated to each of them.
By making the observation shown in Eq. (26) it is easy to see that the two qubit states |D ±1 (t) are direct products of the one qubit states D ± 1 2 (t) up to a mapping t → 2t. This can be used to predict the behaviour of |D ±1 (t) . As stated in Sec. II D ± 1 2 ( tr 2 ) = ψ + 1 att therefore we can conclude that |D ±1 (t) will go to a direct product of ψ + 1 att at t = tr 4 .
So we have found a term that is similar to the one qubit 'attractor' state, ψ + 1 att , which we shall label the two qubit 'attractor' state, ψ + 2 att : In the one qubit case the qubit system went to a second 'attractor' state ψ − 1 att at the time 3tr 2 and means D ±1 ( 3tr

A. Basin of Attraction
The two qubit JCM is more complicated than the original one qubit case and shows some very interesting properties. Firstly we rewrite the wavefunction for the whole system as the time t = tr 4 to understand the two qubit case in more detail.
With the wavefunction written in this form we can see which initial conditions result in the entropy going to zero. Indeed when β 0 = 0 the qubit and field parts are each in a product state, so there is no qubit-resonator entanglement present. There is also no qubit-qubit entanglement as these are in the unentangled 'attractor' state and this means there is absolutely no entanglement in the system at t = tr 4 . Interestingly, the information on the initial state of the qubits is again stored in the radiation field and the qubits are in a spin coherent 'attractor' state (see appendix A).
The initial conditions of the qubit that results in β 0 = 0 are: e iθ C ee = e −iθ C gg = a and C eg = C ge = 1 2 − |a| 2 . So the initial state of the qubits will be of the form: where θ is the initial phase of the radiation field and 0 ≤ |a| ≤ 1 √ 2 . These states lie in the symmetric subspace and have some interesting properties that will be discussed in section III B.
While the qubit state is within this basin of attraction its state is parameterized by the single complex variable a and the wavefunction only contains two field states, |Φ ±1 (t) . At the time tr 4 these two field states are macroscopically distinct coherent states: We now see that the initial state of the qubits alters the relative proportions of |α with respect to |−α . This is also seen in the one qubit case and is the mechanism whereby the quantum information of the initial qubit state is encoded in the field. State preparation of a Schrödinger cat state utilizing this phenomenon has been discussed [14,30]. With the two qubit case the same protocol can be used, but is not so useful as the initial qubit state has to be within the basin of attraction in contrast to the one qubit case where all initial qubit states lead to a Schrödinger cat state at the 'attractor' times.
There is also a connection between the field states |Φ ±1 (t) in the two qubit case and the field states Φ ± 1 2 (t) from the one qubit case, |Φ ±1 (t) = Φ ± 1 2 (2t) . Again there is a mapping t → 2t. In Sec. II it was shown that the field and qubit become unentangled at t r and 2t r as . Using this information from the one qubit case we can predict that Φ +1 ( tr 2 ) = Φ −1 ( tr 2 ) and |Φ +1 (t r ) = |Φ −1 (t r ) which is indeed the case. As a consequence the field part of the wavefunction can be factorised out and there is no entanglement present between the qubits and the field at the times tr 2 and t r . We now consider the numerical results to test the analytical predictions and see how the system behaves whenn = 50, shown in Fig 4. We consider the probability of the qubits being in the state |gg and the probability P 2 att (t) = ψ + 2 att ρ q (t) ψ + 2 att , that ψ + 2 att occurs at a particular time for two sets of initial conditions, one satisfying β 0 = 0 (Fig 4(b)) and one not (Fig 4(a)). Evidently the pattern of collapse and revival is more complicated than in the one qubit case [31], but nevertheless we still see characteristic dips in the entropy. The value of the entropy at tr 4 is not the same in the two cases shown in Fig. 4 and in fact they show two completely different behaviours. In Fig.   4(a) we notice that the dip in the entropy goes to a value of 0.35, very different to the value of the entropy in Fig. 1, but in Fig 4(b) we do indeed find that the entropy is very small at tr 4 . The initial state used for Fig. 4(a) is not in the basin of attraction, whereas the initial state used for A second dip in the entropy can be seen in Fig. 4 at the time t = 3tr 4 which, in the largen approximation, goes to zero. The probability of being in the 'attractor' state ψ + 2 att goes to zero at this time. The field has disentangled itself from the field again, so at this point the qubits are again in a pure state, the orthogonal 'attractor' state ψ − 1 att . It can be seen from the analytic solution that, in the largen approximation, the field and qubits also disentangle themselves at the time t r , for all initial conditions of the qubits, and also at tr 2 if the initial qubit state is in the basin of attraction. This is not seen in the diagrams in Fig. 4 and instead the value to the entropies at these times are fairly high which means the numerically evaluated exact system is far from the analytic results in largen approximation. This is similar to the one qubit case and can be explained by considering the dynamics of the field.
We can once again plot the Q function for the numerical solution to pictorially display what is occurring in the cavity field. The values of Q(α, t) in the complex α plane at fixed times are shown in Fig. 5 when the qubits are initially in the state |gg , which is outside the basin of attraction.
This time the Q function shows a set of three wavepackets of the cavity field (compared to the two for the one-qubit case, Fig. 2). As a function of time the three 'blobs' separate, but only two move around the complex plane and follow the circular path of radius √n . Again the arrows show the dipole moment of the qubit state associated with each wavepacket, The 'blob' that does not move over time has a dipole moment of zero, so only the two moving Gea-Banacloche states have an arrow.
For the one-qubit case the revival time occurred when the two wavepackets overlapped. In this case we have two 'blobs' that move around phase space and overlap at the time t = tr 2 , leaving the stationary wavepacket at the initial position. This causes a rather weak revival peak as only two ically different and when β 0 = 0. Moving from the one qubit case to the two qubit case involves a mapping of t → 2t and a similar mapping can be used to understand the equations for more qubits.

B. Collapse and Revival of Entanglement
An interesting new feature of the two-qubit case, as opposed to the one-qubit case, is that there can be entanglement present between the two qubits themselves, not just between the qubits and the field. So far we have used the entropy to measure the amount of entanglement between the field and the qubit system as the whole system is in a pure state. Unfortunately this measure can not be used to calculate the entanglement between the two qubits as they are sometimes in a mixed state themselves. To find the amount of entanglement between the two qubits we use the mixed state tangle τ , defined as [32,33]: The λ ′ s are the square root of the eigenvalues, in decreasing order, of the matrix 3) and ρ * q represents the complex conjugate of ρ q in the |e , |g basis. A tangle of unity indicates that there is maximum entanglement between the two qubits. If the tangle takes the value of zero, then there is no entanglement present between the two qubits. To clarify, the entropy will be used to measure the amount of entanglement between the field and the qubit system whereas the tangle will be used to measure the entanglement between the qubits.
In this subsection we focus on the initial states inside the basin of attraction. The value of the tangle for the states given by Eq. (33) can be seen in Fig. 6. This diagram shows only two points where τ = 0, for a = ± 1 2 , which demonstrates that there are only two product states in the 'basin of attraction'. In contrast, there are many points in which τ = 1, for a = e iφ √ 2 where φ is an arbitrary phase, or a = ir where r is a real number. This means that all possible levels of entanglement are represented in the basin of attraction, Eq. (33).
Whilst almost all states in the basin of attraction given in Eq. (33) describe entangled states, these all evolve into ψ + 2 att at t = tr 4 and ψ − 2 att at t = 3tr 4 , at which times they are pure and unentangled. So if the system is started off in one of these entangled states defined by Eq. (33), at the time tr 4 the system goes to a state which has zero entanglement. Therefore we have gone from having entanglement in the system to having absolutely no entanglement, either between the qubits and the resonator or between the qubits themselves. At a later time this entanglement returns to the system and can be said to 'revive'. The tangle returns to its initial value at tr 2 and t r , regardless of the phase of the coherent state. This can be seen when we consider the wavefunction for the system at this time: As stated in Sec. III A this state shows no entanglement between the qubits and the field and therefore the entropy of the reduced qubit density matrix will be zero. The value of the tangle for the qubit state at t = tr 2 is found to be exactly the same as the tangle for the initial state of the qubits, therefore τ (0) = τ ( tr 2 ). We have called this phenomenon 'Collapse and Revival of Entanglement' [34]. Analogous to the single qubit case, the quantum information in the initial state of the qubits (for states within the basin of attraction) is swapped out and encoded into the field state at the 'attractor' times. Calculations demonstrate that the value of the tangle remains near zero for long periods between revival. Therefore we are dealing with a phenomenon which has been described as the 'death of entanglement' by Yu and Eberley [35,36] which forms the centre of much current interest [37]. Qing et al. [38] have found a similar collapse and revival for the same model we have studied here but for very different initial conditions. Even just within the system we consider here, detailed studies reveal interesting links between the dynamics of the disappearance (and reappearance) of entanglement and the basin of attraction, as we discuss in this and the subsequent subsection.
So if the qubit system starts in the basin of attraction and τ = 0, over time the initial entanglement in the qubits is exchanged for entanglement between the qubits and the field, and then for a superposition of field states. At the time of the 'attractor' state, tr 4 , there is absolutely no entanglement present in the system. Then the process starts to reverse. Firstly there is again entanglement between the qubits and the field and at the two qubit revival time, tr 2 , τ returns to its initial value and only qubit entanglement is present.
The actual extent of the entanglement revival for the exact solution can be seen in Fig. 7 (a) where we show the tangle for different initial states within the basin of attraction. On first glance we see that the amount of tangle present at t = tr 2 is nowhere near the initial amount of tangle present at t = 0, which is far from the prediction in the largen approximation. The value of the entropy at tr 2 is around 0.5 for the case when a = 1 √ 2 , rather than the predicted value of zero, Fig. 7 (b). In order for there to be maximum possible entanglement in the system there has to be the least amount of entropy, so there is a trade off in tangle for entropy [25].
Another very interesting feature illustrated in Fig. 7 (a) is the manner in which the entanglement vanishes. For all states within the basin of attraction, although the initial entanglement decays rapidly (before reviving at tr 2 and later times) it does so with a Gaussian envelope (as per the Rabi oscillation collapse). It therefore goes smoothly to zero-this is reflected by the further observation that the eigenvalue expression in Eq. (35) never goes negative, so the max operation never needs to be implemented for the tangle. This can be explained by considering the rank of the matrix used to calculate the tangle, ρ q (σ y ⊗ σ y ) ρ * q (σ y ⊗ σ y ). The rank of a matrix is equal to the number of non-zero eigenvalues. The matrix used to calculate the tangle for this system has a maximum of rank= 2 when the qubits are started in the basin of attraction, so λ 3 = λ 4 = 0. Due to the definition of the tangle, λ 1 ≥ λ 2 and Eq. (35) becomes: This characteristic 'Collapse and Revival' behaviour is to be contrasted with behaviour where entanglement vanishes or appears with finite gradient and where the max operation is needed in Eq. (35) because the eigenvalue expression does go negative. To distinguish these behaviours, we use the terminology 'sudden death/birth' of entanglement [35] for the latter. This phenomenon does not occur for initial qubit states within the basin of attraction (which all show collapse and revival of entanglement) but can for those outside, as will be demonstrated in the next subsection. (Bottom)'Collapse and Revival of Entanglement' for many different initial qubit states in the basin of attraction withn = 50. These are shown in many different colours depending on the initial entanglement measure. It shows that the entanglement in the qubits is initially lost then returns at the first revival peak, t = tr 2 . As the initial state is in the basin of attraction we know the whole system has zero entanglement at t = tr 4 . So there is even no entanglement between the field and the qubits. Therefore we know that over time the amount of entanglement varies a great deal.

C. Outside the Basin of Attraction
The previous subsections discussed interesting features of the entanglement dynamics for the two-qubit JCM when considering initial conditions inside the basin of attraction. In this subsection we highlight interesting features that emerge when considering all other cases-when the qubit Hilbert space is not restricted and, for example, the terms (b|eg + 1 − |b| 2 |ge )/ √ 2, where b is a complex number, are present in the initial state, as well as unequal superpositions of |gg and |ee .
To date, in other works it has been found that the JCM for two qubits can exhibit an interesting phenomenon called 'Sudden Death of Entanglement', for certain initial conditions of the qubits.
Some investigations start with the qubits and the field in a mixed state and some others do not [38].
In our work we consider the case when the radiation field is started in a coherent state. Qubits states inside the basin of attraction demonstrate collapse and revival of entanglement, but not sudden death or birth. However, if we venture outside the basin of attraction and include other initial conditions of the two qubit system, we find that there is indeed sudden death of entanglement.
This means that we do not restrict the condition β 0 = 0 in Eq. (23). Now, without this restriction, the qubit system will not go towards the 'attractor' state at tr 4 and 3tr 4 in the largen approximation as β 0 = 0. Furthermore the entropy will be non zero at these times as the qubit state can not be factorised out of the wavefunction. However, the system will still have zero entropy at the 'full' revival time t r as |Φ −1 (t r ) = |Φ 0 (t r ) = |Φ +1 (t r ) = |α and the field state can be factorised out of the wavefunction. Like before τ (0) = τ (t r ) so the amount of entanglement between the qubits is the same at the main revival time and t = 0. For n = 50, this revival will not be complete and so the tangle will not return to its initial value. Fig. 8 shows examples of the qubit entanglement dynamics for initial states which are outside the basin of attraction for the numerically evaluated exact solution. In this figure a new measure of qubit entanglement has been used called the concurrence ζ = √ τ and has been used as a direct comparison to previous work [36,37,39].
Remarkably, Fig. 8 shows that a field initially in a coherent state can give sudden death and birth of entanglement, with initial qubit states of various forms lying outside the basin of attraction.
A case where the initial conditions are in the basin of attraction is shown for reference. Clearly peaks in the entanglement measure at both tr 4 and 3tr 4 arise (with death/birth at finite gradient)times when there would be no entanglement if the initial qubit state was in the basin of attraction.
There are still peaks at the times tr 2 and t r , but these have a lower maximum and also now exhibit birth/death behaviour rather than Gaussian collapse and revival. One initial condition that yields 'Sudden Death of Entanglement' has been shown in Fig. 8 alongside an entropy measurement.
To recap, we find that our system can display both a sudden birth/death of entanglement as well as a Gaussian collapse and revival of entanglement, depending on initial conditions. Furthermore, we understand that, in this system, the sudden death occurs when λ 1 − λ 2 − λ 3 − λ 4 becomes negative and so, according due to the max operation used in Eq. 35 the tangle goes to zero with a finite gradient. In contrast, within the basin of attraction, there are only two non-zero eigenvalues λ 1 > λ 2 , so the tangle can only display a smooth collapse.
These additional peaks in the qubit entanglement seen in Fig. 8 are due to the inclusion of |D 0 and |Φ 0 (t) in the wavefunction. From the analysis in Sec. III we showed that D ±1 ( tr 4 ) = ψ + 2 att , the unentangled two qubit 'attractor' state, so all entanglement present at t = tr 4 is due to |D 0 . If |D 0 is a product state then there will be no entanglement in the qubit system. However if |D 0 is entangled it still does not mean there will be some qubit entanglement as a mixture of entangled and product states need not be entangled. For example curve (b) in the top diagram of Fig. 8 shows there is zero concurrence at tr 2 , even though |D 0 contains some entanglement.

IV. N q QUBIT JAYNES-CUMMINGS MODEL
Having found a sense in which the one qubit 'attractor' state has analogues for two qubits, it is natural to inquire whether similar generalisation holds for the case of more than two qubits. We use the Hamiltonian shown in Eq. (19) for N q qubits. We have found that there are indeed cases where the system goes to a state similar to the one qubit 'attractor' state and have a general form for the basin of attraction for N q qubits.
In the largen approximation the state of the system takes the following form [16]: |D k (t) has a complicated form and occurs in the paper by Meunier et al. [16]. For this work we will be most concerned with the states that appear in the expression for the basin of attraction.
These take the following simple form: Note that: and we can use this equation to understand the simple dynamics of the N q solution. The field initially starts in a coherent state that then splits into N q + 1 different 'blobs' and a revival peak occurs every time at least two of these wavepackets overlaps. The strength of the revival is determined by the number of wavepackets that overlap and these revival peaks occur at definite intervals. Using Eq. (41) we can write a general form for the revival times: k + 1 p t r , k + (p − 1) p t r , p = 1, 2, ..., N q and k = 1, 2, ..., ∞.
Eq. (41) also highlights that when there are an even number of qubits the field returns to its initial state |α after a time of t r . The time till the wavepackets return to their original position is determined by the slowest moving wavepackets, which are e ±2iπt/tr α in the even qubit case and e ±iπt/tr α in the odd qubit case. When N q is odd the wavepackets return to their start position after a time 2t r and there will always be no stationary wavepacket at the start position.
A sketch of the Q function of the field for three qubits is shown in Fig. 9 for large values of n and we once again show the dipole moments as arrows. Using |D k in reference [16] we can compute a general term for the equation of the dipole moments in the largen approximation and is shown below.
where − → σ Nq =σ x Nqî +σ y Nqĵ +σ z Nqk , is the Pauli vector for N q qubits and We notice again the pattern of the arrows changes over time and that the wavepackets come to their original starting position after a time 2t r .

A. Basin of Attraction
Using the above experience with one and two qubits, we investigate whether a basin of attraction exists for N q > 2. Using the expansion pioneered by Chumakov et al. [15] and developed a : t 0 Re Α Im Α

e : t tr 2
Re Α Im Α

f : t 2tr 3
Re Α Im Α g : t 5tr 6 Re 5t independently by Meunier et al. [16] we have indeed found a basin of attraction for N q qubits. As D ± 1 2 ( tr 2 ) = ψ + 1 att , and using the observation made in Eq. (40), we can therefore conclude that D ± Nq 2 (t) will go to a direct product of ψ + 1 att at t = tr 2Nq . then only D ± Nq 2 (t) are in the wavefunction and the qubits will go to ψ + Nq att at t = tr 2Nq . Using this information we have found a basin of attraction for N q qubits.
As before φ is an arbitrary phase and θ is the initial phase of the radiation field with 0 ≤ |a| ≤ 1 √ 2 Nq −1 . The basin of attraction is in the symmetric subspace and |N q , m is the symmetrised state with N e excited qubits, N g qubits in the ground state and m = Ne−Ng 2 . For example, for the three qubit case, N q = 3, m = 1 2 = 1 √ 3 (|eeg + |ege + |gee ) as N e = 2 and N g = 1 etc. We note that although the state space of the qubits increases with increasing N q , the basin of attraction is still parameterised by the single complex value a. The significance of this will become apparent in Sec.

IV B.
We emphasise that when the qubits start in the basin of attraction then there are only two wavepackets in phase space of the photon field for all time. The initial conditions have restricted the system to behave like the one qubit case as is highlighted by considering the wavefunction at In the one qubit case the qubit system also goes to a second 'attractor' state ψ − 1 att at the time = ψ − 1 att ⊗Nq = ψ − Nq att . Qubits started in the basin of attraction, Eq. (48), will go to these 'attractor' states and hence the entropy will go to zero.  (t) = ψ ± Nq att . Note that only revival and 'attractor' times in the interval 0 < t ≤ t r are shown. Further times can be obtained by adding integer multiples of t r .
We can construct a series for when these 'attractor' times will occur for N q qubits, shown in table I with the revival times from Eq. (42). The general equation for these times is: .., N q and k = 1, 2, ..., ∞.
In Fig. 10 we show the 3 and 4 qubit evolution starting in a state in the basin of attraction (Eq. (48)). In general we do see there are dips in the entropy that go close to zero and the probability of being in the 'attractor' state goes to close to one. These features get more pronounced asn is increased. As for the two-qubit case, when this probability goes to zero, the system is in the second 'attractor' state, ψ − Nq att , which is orthogonal to the first.

B. Collapse and Revival of 'Schrödinger Cat' States
Surprisingly, the basin of attraction has a simple form when written in terms of spin coherent states, which are analogous to coherent states of a photon field. As described in Appendix A, such states are given by: where |N q , m is defined in Sec. IV A. A spin coherent state is a product state and this is easily seen when it is written in Majorana representation [40] in Eq.
From Eq. (52) it is clear that z = e −iθ , N q is orthogonal to z = −e −iθ , N q and ψ Nq is a collection of spin 'Schrödinger cat' states. We also write the state of the field at the time of the first dip in the entropy for an arbitrary number of qubits, t = tr 2Nq : the initial qubit state has moved from the qubit system to the field system [34].
To illustrate this migration we use a function similar to the field Q function, Q q , using the completeness relation in appendix A, Q q (z, t) corresponding to ψ Nq in Eq. (53) is depicted in Fig 11. There are indeed two separate peaks that correspond to two macroscopically different spin coherent states.
The 'attractor' states for N q qubits in Eq. (48) are also spin coherent states (see appendix A).
In Fig 12 we plot a spin Q q function diagram for the 'attractor' state.
So far we have shown, for largen, that the qubit states in the basin of attraction are spin Schrödinger cat states when a = 1 √ 2 Nq . At a later time, t = tr 2Nq , the qubits are in the spin coherent 'attractor' state and the field is in a Schrödinger cat state. Then the qubits return once again to a spin Schrödinger cat state at the time t = tr Nq . This is shown in Table II. So we conclude that the 'Schrödinger cat' state migrates from the qubits to the field and back again.

C. Collapse and Revival of Entanglement
Whilst the basin of attraction for the N q qubit case is a very small part of the Hilbert space of all initial states it shows many different degrees of entanglement, just like the two qubit case.
Some values of a and φ give the qubit state to be in a product state and other values give some degree of entanglement between the qubits. When considering more than two qubits there is no universally accepted measure of entanglement. Indeed, for N q ≥ 4 and multilevel systems there exist infinitely many inequivalent kinds of entanglement [41]. GHZ type entanglement exists for all numbers of qubits and, up to local unitaries, a maximally entangled GHZ state has the form: where b|b ⊥ = 0. In the previous section it was stated that the basin of attraction can be rewritten as a spin 'Schrödinger cat' state which is a state that is a superposition of two macroscopically distinguishable states. However expanding states in the basin of attraction using Eq. (52) we see that it has a GHZ form: although it is only a maximal GHZ state when a = 0 and a = For three qubits there are two different types of entanglement, GHZ [42] and W [43]. |GHZ = 1 √ 2 (|eee + |ggg ) can be regarded as a maximally entangled state, but if one of the qubits is traced out the remaining state has absolutely no entanglement present. The entanglement of |W = 1 √ 3 (|eeg + |ege + |gee ) is maximally robust and is not so fragile to particle loss (i.e. there is still entanglement present after tracing out any individual qubit). It is not possible to transform |GHZ into |W even with a small probability of success and vise versa, and thus the two states represent fundamentally different types of entanglement. In particular, the GHZ state represents a truly 3-body entanglement. Similarly, any N q -qubit GHZ state contains specifically N q -body entanglement.
We have used a pure state GHZ entanglement measure introduced by Coffman et al. [44]: where τ ABC is the amount of GHZ entanglement, C 2 A(BC) = 2 √ detρ A and C 2 AC is the entanglement after qubit B has been traced out. The results for different values of a for the three qubit version of Eq. (48) are shown in Fig. 13. This measurement only works for three qubit states that are pure, so we only perform the calculation on the initial state. This plot shows the same shape as that of Fig. 6, but values of a are now limited to a smaller range. On the basis of Figs. 6 and 13, and using Eq. (58), we conjecture that for more qubits the basin of attraction will always contain maximally entangled states that involve all the qubits and two unentangled states, but the size of the basin will get smaller as the number of qubits increases.
For the largen approximation it was shown in Sec. III B that if a state is started in the basin of attraction τ (0) = τ ( tr 4 ). This can also be extended to the three qubit case and τ ABC (0) = τ ABC ( tr 6 ).
Therefore if the initial qubit state is in the basin of attraction and τ ABC = 0 then the whole entanglement will return to the qubit system at t = tr 6 . To see if 'Collapse and Revival' of entanglement also occurs in the exact solution for three qubits we start the qubits in an initial state, |ψ 3 , that is a maximally entangled GHZ state (τ ABC = 1).
We then calculated the probability P init (t) = ψ 3 |ρ q (t)|ψ 3 , which is the probability the qubits are in the initial state. If the probability collapses and revives over time then we can conclude that the three qubit system also exhibits 'Collapse and Revival' of entanglement. The results are shown in  The qubits are started in a maximally entangled GHZ state, (τ ABC = 1), that lies in the basin of attraction and are then allowed to evolve. At a later time the entropy goes close to zero and the qubits are close to the product form shown in Eq. (47). This is a completely unentangled state and as the entropy is also zero in the largen approximation, there is no entanglement in the entire system. After the 'attractor' time the entropy increases and is no longer close to zero. We follow the green curve in Fig. 14 and note that P init (t) increases to a value near 0.8 and therefore we propose there is some three qubit mixed state entanglement present in the qubit system. In the largen approximation this value does go to unity, so the system is in the initial state, up to a phase. In this approximation we will therefore necessarily get a full revival of three-body GHZ entanglement.
As further confirmation that the entanglement is of 3-body form, we have also calculated the entanglement present between just two of the qubits after tracing out one of the qubits. It was found that the tangle between any two qubits remained zero throughout the evolution. Therefore we can conclude the entanglement present is just GHZ entanglement and there is no W component present (as a W component would lead to non-zero 2-qubit entanglement).
We have demonstrated the revival of entanglement for the numerically evaluated exact solution for two and three qubits, and shown it to be true for all values of N q in the largen approximation.
Of particular note is that the type of entanglement we observe is of the GHZ form and is thus a true N q -qubit entanglement that cannot be reduced to a simpler form. Remarkably, this implies that the quantum information encoded in an N q -qubit GHZ entangled state-for any value of N q -can be encoded into the state of a single field mode.

V. CONCLUSION
We have studied the 'Collapse and Revival of Entanglement' in the two, three and many qubit JCM. We have shown that the idea of an 'attractor' state discovered by Gea-Banacloche [13] for one qubit does indeed exist in the multi qubit cases described by the Tavis-Cummings model.
Unlike the one qubit case, the qubits have to be within a basin of attraction in order for the qubit system to disentangle itself from the field at certain times. When the qubits are started in a state in the basin of attraction then the system behaves in a similar way to the one qubit case.
The basin of attraction contains mostly entangled states and there are some interesting entanglement properties associated with the whole system when the qubits are started in this basin of attraction. The system evolves and there is then entanglement between the field and the qubit system, but when the system goes to the 'attractor' state all entanglement has completely left the system. The field and the qubits are both in pure states with the qubit state being in the completely unentangled 'attractor' state. At a later time it has been found that the entanglement does return to the system and then the process repeats. We represented explicit examples of entanglement present between the qubits and also between the qubit system and the field. We have called this phenomenon 'Collapse and Revival of Entanglement'. For initial qubit states outside the basin of attraction, sudden death and rebirth [37] of entanglement can occur.

APPENDIX A: SPIN COHERENT STATES
There exist spin states analogous to the coherent state of the harmonic oscillator. For the one dimensional harmonic oscillator the coherent states are functions of a variable α which runs over the entire complex plane: where the normalisation factor N = e |α| 2 /2 .
These states from a complete set, in the sense that [7] 1 π |α α|d 2 α = ∞ n=0 |n n| = 1 (A4) where the right hand side is the unit matrix, or identity operator.
To find the analogue for spin half particles we consider a single particle of spin S = N q /2. We define the ground state |0 as the state such thatŜ z = S|0 whereŜ z = Nq i=0 σ z i is the operator of the z component of spin. Then the operatorŜ + =Ŝ x +iŜ y 2 creates spin deviations. In a similar way to the field coherent state the spin coherent state for many spin half particles is [45,46,47]: In our work we present the 'attractor' states for many qubits in Eq. 56 as spin coherent states.
These states themselves are indeed spin coherent states. If we set z = ie iθ in Eq. A7 then we get ψ + Nq att and if we set z = −ie iθ then we get ψ − Nq att . When looking at the coherent state of a harmonic oscillator, there exists a set of states called Schrödinger cat states that are the addition of two coherent states that are macroscopically different, e.g. 1 √ 2 (|α + |−α ). There also exist spin coherent states that are analogous to these Schrödinger cat states. They are also defined as the addition of two macroscopically different states, i.e. 1 √ 2 (|z, N q + |−z, N q ). In our work we have plotted Q function diagrams for the cavity field state using the coherent state. There also exists an analogous function for the qubits, defined in terms of spin coherent states in a similar way, using the completeness relation Eq. A8. This spin Q function is given by Q q (z, t) = z|ρ q (t)|z 1 + |z| 2 2 . (A9) We have used this in our analysis to make interesting plots which demonstrate pictorially that the states in the basin of attraction for multi-qubit systems are indeed Schrödinger cat states of the qubits.

ACKNOWLEDGMENTS
The work of C.E.A.J. was supported by UK HP/EPSRC case studentship, and D.A.R. was supported by a UK EPSRC fellowship. We thank W. J. Munro for helpful discussions.