Entanglement through qubit motion and the dynamical Casimir effect

TL;DR

This paper investigates how qubit motion and the dynamical Casimir effect can generate entanglement in superconducting quantum systems, highlighting the role of photon correlations and relativistic concepts.

Contribution

It introduces a novel analysis of entanglement generation via the DCE considering different qubit motions, connecting quantum optics with relativistic effects.

Findings

Photon absorption and emission are key to entanglement.

Entanglement can be related to the concept of simultaneity.

The study bridges superconducting quantum technologies with relativistic quantum effects.

Abstract

We explore the interplay between acceleration radiation and the dynamical Casimir effect (DCE) in the field of superconducting quantum technologies, analyzing the generation of entanglement between two qubits by means of the DCE in several states of qubit motion. We show that the correlated absorption and emission of photons is crucial for entanglement, which in some cases can be linked to the notion of simultaneity in special relativity.

Figures

Equations

$$H$$

$$+\hslash \sum _{i=1}^2{g}_i\cos (k_i{x}_{q_i}(t)){\sigma }_i^x(a_i^{†}+a_i)+$$

$$+\hslash g_{1,2}(t)(a_1^{†}+a_1)(a_2^{†}+a_2),$$

$$g_{1,2}(t)(a_1^{†}+a_1)(a_2^{†}+a_2)\to g_0/2(a_1^{†}{a}_2^{†}+a_1{a}_2).$$

$$H$$

$$+\hslash \sum _{i=1}^2{g}_i\cos (k_i{x}_{q_i}(t)){\sigma }_i^x(a_i^{†}+a_i)+$$

$$+\hslash g_0/2(a_1^{†}{a}_2^{†}{e}^{i{\omega }_{d}t}+a_1{a}_2{e}^{-i{\omega }_{d}t}),$$

$${\rho }_{qubits}^{(3)}=\left(\begin{array}{cccc}0& 0& 0& {\rho }_{14}\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ {\rho }_{41}& 0& 0& 1\end{array}\right)$$

$$C(\rho )=2|{\rho }_{14}|=2|⟨ee|\rho |gg⟩|=2|⟨00ee|{\rho }_{total}|00gg⟩|$$

$$=g_1{g}_2{g}_0\times$$

$$\displaystyle\Biggl{|}\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}dt_{2}\int_{0}^{t_{2}}dt_{3}\cos(f_{1}(t_{1}))\cos(f_{2}(t_{2}))e^{i\omega_{d}t_{3}}+$$

$$\displaystyle+\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}dt_{2}\int_{0}^{t_{2}}dt_{3}\cos(f_{1}(t_{2}))\cos(f_{2}(t_{1}))e^{i\omega_{d}t_{3}}\Biggr{|},$$

$$C_{\text{rest}}=\frac{g_0{g}_1{g}_2}{{\omega }_d}{t}^2+O(t).$$

$$C_{\text{v const}}$$

$$\displaystyle\Biggl{|}\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}dt_{2}(e^{ik_{1}v_{1}t_{1}}+e^{-ik_{1}v_{1}t_{1}})\times$$

$$(e^{i{k}_2{v}_2{t}_2}+e^{-i{k}_2{v}_2{t}_2})(e^{i{\omega }_d{t}_2}-1)+$$

$$+{\int }_0^{t}d{t}_1{\int }_0^{t_1}d{t}_2(e^{i{k}_1{v}_1{t}_2}+e^{-i{k}_1{v}_1{t}_2})\times$$

$$\displaystyle(e^{ik_{2}v_{2}t_{1}}+e^{-ik_{2}v_{2}t_{1}})(e^{i\omega_{d}t_{2}}-1)\Biggr{|}$$

$${\omega }_d$$

$$k_1|v_1|$$

$$c/L_1+c/L_2$$

$$|v_1|/L_1$$

$$C_{\text{v const, 1}}=\frac{g_0{g}_1{g}_2}{{\omega }_d^2}|\sin ({\omega }_{d}t)|t+O(t^0),$$

$$C_{\text{v const, 2}}=\frac{g_0{g}_1{g}_2}{2{\omega }_{d}kv}|\sin (kvt)|t+O(t^0).$$

$$x_1(t)$$

$$x_2(t)$$

$$C_{\text{arccos}}=\frac{4g_0{g}_1{g}_2}{{\omega }_d(n+1{)}^2}\frac{v^{2n}}{L^{2n}}{t}^{2n+2}+O(t^{2n}).$$

$$f_n(x)$$

$$x_1(t)$$

$$=-L_1+L_1{f}_n\left(\frac{t}{\tau }-\lfloor \frac{t}{\tau }\rfloor \right)$$

$$x_2(t)$$

$$=L_2-L_2{f}_n\left(\frac{t}{\tau }-\lfloor \frac{t}{\tau }\rfloor \right)$$

$${\rho }_{qubits}=\left(\begin{array}{cccc}{\rho }_{11}& 0& 0& {\rho }_{14}\\ 0& {\rho }_{22}& 0& 0\\ 0& 0& {\rho }_{33}& 0\\ {\rho }_{41}& 0& 0& {\rho }_{44}\end{array}\right)$$

$${\rho }_{11}=⟨ee|\rho |ee⟩$$

$${\rho }_{22}=⟨eg|\rho |eg⟩$$

$${\rho }_{33}=⟨ge|\rho |ge⟩$$

$${\rho }_{44}=⟨gg|\rho |gg⟩$$

$${\rho }_{14}=⟨ee|\rho |gg⟩$$

$${\rho }_{qubits}=\left(\begin{array}{cccc}0& 0& 0& c_{00ee}^{(3)\ast }\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ c_{00ee}^{(3)}& 0& 0& 1\end{array}\right)$$

Full text

I Introduction

Superconducting quantum technologies together with circuit quantum electrodynamics (cQED) Wallraff et al. (2004) form one of the most promising candidates for processing quantum information as well as the experimentation on the foundations of quantum mechanics. The advantages of this technology are, on the one hand, the strong coupling between the superconducting qubits and the resonant cavities, as well as the ability to widely tune this coupling and all the parameters of the system, allowing the investigation of new phenomena Astafiev et al. (2010); Niemczyk et al. (2010) such as, among many others, the Dynamical Casimir Effect (DCE) Wilson et al. (2011) or the Unruh effect Nation et al. (2012).

The DCE is a member of a large family of effects linked to the quantum fluctuations of the vacuum, among which are the Lamb displacement Lamb and Retherford (1947), the magnetic moment of the electron or the Casimir (static) effect Casimir (1948); Lamoreaux (1997); Mohideen and Roy (1998). The latter is produced by a reduction in the density of modes imposed by certain boundary conditions, which leads to a pressure of radiation exerted by the vacuum. Its dynamical counterpart shares a similar origin, except for the fact that boundary conditions must be time-dependent, which can be achieved by means of SQUIDs Wilson et al. (2011) . This effect has also been measured by modulating the effective speed of light Lähteenmäki et al. (2013). On the other hand, the DCE can be considered a resource to generate quantum correlations, including quantum entanglement Andreata and Dodonov (2005); Dodonov and Dodonov (2012); Felicetti et al. (2014); Sabín et al. (2015); Samos-SaénzdeBuruaga and Sabín (2017). In this paper, we will focus on this feature.

The Unruh effect is another member of the family of aforementioned phenomena, since it consists in the measurement of thermal radiation by a detector moving with constant proper acceleration through the quantum vacuum Davies (1975); Unruh (1976). Unfortunately, this has never been observed since it requires unreachable accelerations to generate detectable signals, although it can be increased by several orders of magnitude when two-level systems are accelerated through resonant cavities by means of non-adiabatic boundary conditions Scully et al. (2003). However, this modification of the original effect is still difficult to reach experimentally and has not been observed. In this report we will analyse another scenario: it is possible to simulate the movement of a qubit in a cavity by modulating the qubit-cavity coupling in the same way said movement would Felicetti et al. (2015), leading in general to the excitation of qubit and cavity from vacuum as in both variants of the effect. In the same way as with the DCE, this simulation of the Unruh effect can be used as a mechanism for the generation of entanglement García-Álvarez et al. (2017), only this time between a qubit and the photons of a cavity, not just in the field as in the first effect.

In this paper we consider an scenario where both the DCE and the simulated enhanced Unruh radiation can take place (see Fig. 1). It consists of two superconducting resonators sharing a common SQUID and each of them coupled to a superconducting qubit. DCE radiation can be generated by means of the modulation of the magnetic flux threading the SQUID, while the coupling of the qubits can be tuned to simulate their motion. We show that the correlated absorption and re-emission of the DCE radiation by the qubits is crucial for the generation of entanglement, as in the case where the qubits are static Felicetti et al. (2014). If the motion of the qubit preserves the correlated nature of the absorptions and emissions, entanglement is preserved as well. However, in general, uncorrelated motion of the qubits will result in the vanishing of entanglement, even for low simulated velocities. In the case of equal-length cavities, this physics can be linked to the notion of simultaneity in special relativity: breaking down simultaneity in the absorption would make entanglement vanish.

The rest of the paper has been structured as follows. First, we present the superconducting setup, its Hamiltonian and some relevant features. Then in the next section, we discuss the results of, on the one hand, applying perturbation theory to the calculation of the concurrence, and on the other hand, calculating this same magnitude solving the master equation of the system numerically. Finally, we summarize our conclusions.

II Setup: DCE and simulated qubit motion

We consider a system composed of two superconducting qubits -in particular, modifications of the usual design of a transmon qubit Srinivasan et al. (2011)- whose coupling to the electromagnetic field can be controlled by the magnetic flux threading the SQUIDs that compose them, which offers an opportunity to simulate the generation of acceleration radiation in the cavity-enhanced Unruh effect. In addition, each of these qubits will interact directly with only one of the modes of only one of the two cavities or transmission line resonators of the system. These resonators interact in turn with each other by means of another SQUID, which allows to produce time-dependent boundary conditions in the resonators, giving rise to the DCE. An outline of this system can be found in Fig 1.

The Hamiltonian is

[TABLE]

where the sum runs over both cavities and both qubits. The first line contains the static Hamiltonian, where $\omega_{c_{i}}$ is the frequency of cavity $1$ or $2$, depending on the subscript, as well as $\omega_{q_{i}}$ is that of the qubits. On the other hand, $a^{\dagger}_{i}$ and $a_{i}$ are the creation and annihilation operators of the corresponding cavity and $\sigma^{z}_{i}$ the third Pauli operator of each qubit. The second line of the equation contains the interaction of the qubits with their cavities, being $k_{i}$ the wave vector of the cavity, $g_{i}$ the maximum intensity of its coupling, $\sigma^{x}_{i}$ the first Pauli operator of the qubit and $x_{q_{i}}(t)$ simulated trajectory of qubit motion. Experimentally, the product $k_{i}x_{q_{i}}(t)$ is actually $f=\phi(t)/\phi_{0}$ with $\phi_{0}$ the quantum of magnetic flux and $\phi(t)$ the flux through the SQUID Felicetti et al. (2015); García-Álvarez et al. (2017) that may be controlled from the outside with a typical nanosecond resolution. The third line contains the interaction between the cavities, being $g_{1,2}(t)$ the time-dependent coupling assuming that the boundary conditions produced by the SQUID do not destroy the structure of normal field modes or make resonances of new modes with the qubits. This is the case when $g_{1,2}(t)=g_{0}\cos(\omega_{d}t)$ with $\omega_{d}=\omega_{c_{1}}+\omega_{c_{2}}$ matching the sum of the cavity mode frequencies. Moreover, in this case the interaction Hamiltonian can be approximated by a two-mode squeezing Hamiltonian Felicetti et al. (2014):

[TABLE]

So, finally:

[TABLE]

Previous work Felicetti et al. (2015) has shown that, even for moderate values of the coupling, the modulation of the coupling strength might resonate with the counterrotating terms of the Hamiltonian. Therefore, we will not perform the rotating wave approximation (RWA).

Finally, notice that if either or both of the trajectories $x_{i}(t)$ are changed by $L_{i}-x_{i}(t)$ and the relevant coupling constant $g_{i}\rightarrow-g_{i}$ is inverted, then the Hamiltonian does not undergo any change, which is quickly deduced using the expression of the cosine of the sum and substituting $k_{i}=\pi/L_{i}$. This symmetry can be interpreted as a mechanism by which a path $x_{q_{i}}(t)$ that passes between the two ends of the cavity in a finite time can be extended beyond that time reflecting it with respect to the center of the cavity. That is, if for a time $\tau$ we have $x_{i}(\tau)=L_{i}$ and at a later time $x(\tau+\delta t)$ exceeds $L_{i}$, then the path can be modified as follows $x_{i}^{\prime}(\tau+\delta t)=L_{i}-x_{i}(\delta t)$ , producing the same Hamiltonian which governed the evolution up to $\tau$, except for the sign of $g_{i}$. In other words, this symmetry offers a natural bounce condition to continue trajectories that reach the ends of the cavities, so natural that it will prove useful throughout the work.

III Results

We will use two different methods to analyze entanglement generation in the setup discussed in the previous section, namely perturbation theory to obtain an approximate expression of the concurrence and then numerical simulations to integrate the master equation governing the system by means of the Python package QuTiP Johansson et al. (2013). These two methods make up the following two subsections. As a final section, dissipation and temperature are addressed.

III.1 Perturbative results

The global state $\ket{\Psi(t)}$ under Hamiltonian (II) is expanded up to third order in the couplings $g_{1}$, $g_{2}$, $g_{0}$ from the ground state $\ket{00gg}$, where $\ket{0}$, $\ket{1}$ and so on are the Fock number states of each cavity and $\ket{g}$, $\ket{e}$ are the ground and excited states of each qubit. We aim to compute the entanglement dynamics of the qubits, so it is necessary to trace over the cavity field states, which leads to the following density matrix:

[TABLE]

where $\rho_{14}$ the matrix element $\bra{ee}\rho\ket{gg}$ and the not zero entry in the diagonal corresponds to $\bra{gg}\rho\ket{gg}$. For a more detailed derivation refer to Appendix A. Then, using the concurrence $C(\rho)$ Hill and Wootters (1997) as entanglement measurement, we find:

[TABLE]

where $f_{i}=k_{i}x_{i}(t)$. It should be noted that in this system the concurrence can be interpreted in a very intuitive fashion: it is proportional to the probability that two photons are emitted, one in each cavity, and that each qubit absorbs one. At first sight, Eq. (III.1) seems to conclude that the qubit motion will only reduce concurrence, since integrating cosines will not provide any further contribution. However, these can cause a resonance, either with each other or with the term due to the emission. In the following paragraphs, we will study these resonances produced by different trajectories, paying attention to the conditions that must be fulfilled for their existence.

III.1.1 Static qubits

If the qubits are static Eq. (III.1) further simplifies, since the cosine functions are all constant. By performing the first integration, we get $\int_{0}^{t_{2}}dt_{3}e^{i\omega_{d}t_{3}}=(e^{i\omega_{d}t_{2}}-1)/i\omega_{d}$ and the constant term will eventually give rise to:

[TABLE]

This quadratic behaviour of the concurrence seems to be in agreement with Felicetti et al. (2014) for the moderate values of time where the perturbative approach is valid -the perturbative approximation will eventually break down, as we will see in detail below.

III.1.2 Constant velocity

Particularizing Eq. (III.1) for the case of qubits moving with constant velocities $v_{i}$ with initial positions at $x=0$, we get :

[TABLE]

By inspection of Eq. (III.1.2), we find that resonances will appear if either or both following conditions are verified:

[TABLE]

Note that $\omega_{d}=\omega_{1}+\omega_{2}$, so the above conditions turn into:

[TABLE]

The first one implies a superluminal velocity of at least one of the qubits, and was already found in Sabín et al. (2017). It is related with the emission of Ginzburg radiation at superluminal constant velocities. The second one -since both qubits start in the same position- means that the distance of the qubits to $x=0$ -in units of the corresponding cavity length- is always the same for both qubits $x_{1}/L_{1}=-x_{2}/L_{2}$. Then the hamiltonians of both qubits are equivalent at any time, entailing that the absorption and re-emission of DCE photons is perfectly correlated. In particular, if the lengths of the cavities are equal $L_{1}=L_{2}$ then the distances are exactly the same, which means that absorptions and emissions occur always simultaneously. This suggests an interesting link between the generation of entanglement in a quantum setup and a key notion of special relativity, such as simultaneity.

When both conditions in Eq. (III.1.2) are met at the same time, namely $|v_{1}|/L_{1}=|v_{2}|/L_{2}=c/L_{1}+c/L_{2}$, then we find:

[TABLE]

while if only the first condition is met $|v_{1}|/L_{1}=c/L_{1}+c/L_{2}\neq|v_{2}|/L_{2}=kv/\pi$, then

[TABLE]

In the latter case, it is the velocity of the non-resonant qubit and not the mirror frequency what modulates the generation of concurrence. Even if we assume $kv\approx\omega_{d}$, then the concurrence in Eq. (11) is reduced by a factor 1/2 with respect to the concurrence in Eq.(10), which highlights again the importance of the correlations among the absorptions and emissions of photons.

Interestingly, we can use the symmetry of the Hamiltonian described at the end of Section III, that in this type of trajectory translates into an inversion of the speed of the qubits when they try to leave the cavities. The full symmetry would change the sign of the relevant coupling, but these signs can be canceled if both trajectories arrive at the ends of the cavities simultaneously and their speeds are then inverted. With these bounces, trajectories of constant velocity can be extended in time. The concurrence (III.1.2) inherits this Hamiltonian symmetry: after $n$ bounces the generated entanglement is $n$ times the entanglement after the first bounce. We see that in this case the simultaneity in the bounce plays a crucial role by simplifying the calculation of the concurrence and extending their relevance. Finally, if the second condition in Eq. (III.1.2) is met but not the first it is not possible to obtain a closed analytical expression of the concurrrence. However, we will present numerical results in the next section.

III.1.3 Other trajectories

A convenient family of qubit trajectories would be given by:

[TABLE]

$\tau$ being the flight time of the qubits, namely the time that it takes for each qubit to traverse its cavity. These trajectories can be seen in Fig. (2). They exhibit divergences in the velocity and the acceleration when time gets close to the flight time, and also at $t=0$ for the particular case $n=1$.

We find the following concurrence:

[TABLE]

Therefore, with these trajectories we are able to produce resonances with arbitrary powers of time. Of course, all the above is restricted by the perturbative approximations adopted, which will eventually break down. We now proceed to show results obtained by a numerical integration of the master equation of the system, which enables the exploration of the long-time dynamics.

III.2 Numerical results

III.2.1 Static qubits and constant velocities

We consider four trajectories (see Fig. (3)).

The first one is the case where the qubits are static, already discussed in Felicetti et al. (2014). Then, we consider the trajectory analyzed perturbatively in the previous section, where the qubits start at $x=0$ and move with opposite velocities, giving rise to correlated -simultaneous for equal cavities- absorption and emission of photons. The third case is related to the second by the symmetry relation discussed throughout this work, since one qubit starts at the other end of its cavity. Finally, we consider a trajectory which breaks the correlations among the absorptions and emissions of photons, since one qubit starts out at the center of its cavity while the other starts at 0.

In Fig. 4 we show the numerical results for the concurrences at long times. We reproduce the results for the static case in Felicetti et al. (2014) up to the point of maximum concurrence, where they propose to switch off the coupling in order to optimize the entanglement generation. We show that in the two trajectories which preserve the correlated absorptions and emissions the high concurrence is indeed preserved, even achieving larger maximum values. Finally, in the asymmetric case, the concurrence is significantly reduced, as expected. Interestingly, this effect occurs already at low non-relativistic velocities, which highlights the role of the correlation/simultaneity in the generation of entanglement in this setup.

In order to achieve the regime of large velocities, it is convenient to use trajectories similar to the ones in Section III.1.3. We will explore them in the next subsection.

III.2.2 Other trajectories

In this case, we have considered the following trajectories (see Fig. 5). We first use similar trajectories as in Eq. (12), but extended by means of the bounce symmetry:

[TABLE]

$\tau=v_{1}/L_{1}=v_{2}/L_{2}$ being again the time of flight which takes the qubits to traverse their cavities and $\lfloor{x}\rfloor$ being the floor function. In this case, as it can be seen in Fig. 5, the trajectories are synchronized in such a way that the absorption and emission of photons is correlated.

Another case that we consider is obtained by replacing in Eq. (14) $x_{2}(t)$ by $x_{2}(t+\tau/2)$, which preserves the symmetry between the qubits. Finally, making instead the replacement with, for instance, $x_{2}(t+0.1\tau)$, then the qubits are out of phase. We show the concurrences for these trajectories in Fig. 6, which again shows that the correlations/simultaneity between the qubits are crucial to understand the magnitude of entanglement generation. In order to further illuminate this point, it is interesting to discuss the population in the Bell basis. In the static case and all the cases where the qubit positions are correlated in the way discussed above, photons are emitted and absorbed in pairs, and therefore it is expected that all the population is in the Bell states $\ket{\phi_{\pm}}=1/\sqrt{2}(\ket{gg}\pm\ket{ee})$. However, if the positions of the qubits are uncorrelated it is possible that one of the qubits emits or absorbs a photon while the other does not. This enables the population of the other Bell states $\ket{\psi_{\pm}}=1/\sqrt{2}(\ket{ge}\pm\ket{eg})$. In Fig. 7 we confirm that this is indeed the case: only in the low-concurrence case a significant population eventually appears in $\ket{\psi_{\pm}}$. Comparing Fig. 6 with 7, we see that jumps in the population of $\ket{\psi_{\pm}}$ are correlated with falls in the value of the concurrence, as expected.

III.3 Dissipation and Temperature

As a final remark to this analysis, dissipation and temperature must be taken into account. The typical temperature in superconducting circuits experiments is of the order of $10-50$ mK, giving a thermal photon number of $10^{-9}-10^{-2}$. This is clearly negligible and hasn’t been considered in the calculations. Analogous simulations to those presented in this letter have been run considering relaxation and dephasing times up to $10^{4}$ ns and resonator decoherence times of $10^{5}$ ns, values at experimental reach by superconducting circuits technology Devoret and Schoelkopf (2013), as are the values of the rest of parameters such as characteristic frequencies and couplings. This dissipation proved to be irrelevant in the dynamics of the system, as initially supposed.

IV Conclusions

We have analyzed the entanglement dynamics between two qubits in a system where each one interacts with a resonant cavity with tunable coupling, which allows to simulate their motion. The cavities interact in turn with each other through a SQUID, which implements a boundary condition that can be modulated by the magnetic flow threading it. This results in a two-mode squeezing hamiltonian which is the source of entanglement. We show that a high degree of entanglement can be generated both in the case where the qubits are static -previously discussed in Felicetti et al. (2014)- and in the cases where the motion preserves the fact that photons are absorbed and re-emitted in pairs -one by each qubit- populating only the Bell states $\ket{\phi_{\pm}}$. Otherwise, if the motion of the qubits is such that photons can be emitted or absorbed by only one qubit, we find that the states $\ket{\psi_{\pm}}$ are also populated and the concurrence is dramatically reduced. If the cavities are equal in length, this means that high-concurrence trajectories are characterized by simultaneous absorption and emission of photons, which suggests an interesting link with the notion of simultaneity in special relativity.

These results pave the way for the exploration of special relativistic effects in a quantum setup. For instance, we envision the quantum simulation of the gedanken textbook experiments where trains moving at relativistic speeds are used to illustrate the relativity of simultaneity. In our setup, the magnitude of entanglement could be used, in principle, as a witness of simultaneity, and vice versa. Our results are fully within reach of current technology.

Acknowledgements

A.A. and C. S have received financial support through the Junior Leader Postdoctoral Fellowship Programme from “la Caixa” Banking Foundation (LCF/BQ/LR18/11640005). E.S. acknowledges support from the projects QMiCS (820505) and OpenSuperQ (820363) of the EU Flagship on Quantum Technologies, Spanish MINECO/FEDER FIS2015-69983-P and Basque Government IT986-16.

Appendix A

As pointed out at the beginning of section III.A, the global state of both qubits and resonator modes is expanded up to third order in the couplings which gives non-zero projections onto the states: $\{\ket{00gg},\ket{11gg},\ket{10ge},\ket{01eg},\ket{00ee},\ket{22gg},\ket{21ge},$ $\ket{12eg},\ket{33gg}\}$. When performing the partial trace over the fields many of this states (those with the qubits in the same state) will mix, leading to a density matrix with the following expression:

[TABLE]

Where the notation for subscripts is the same as for equation (4). The expression of these matrix elements in terms of the perturbative coefficients of $\ket{\Psi}$ are:

[TABLE]

Where the superscripts indicate the order in the perturbative expansion at which the coefficient appears -or in the case of multiple indices, the orders at which new contributions must be taken into account. Notice that since $\ket{\Psi}$ has been expanded up to third order, any term in the matrix elements with a greater power in the couplings (that is, adding their superscripts) must be considered negligible. With this in mind, the qubits’ density matrix reads:

[TABLE]

Which proves equation (4).