Atoms in separated resonators can jointly absorb a single photon

The coherent nonlinear process where a single photon simultaneously excites two or more two-level systems (qubits) in a single-mode resonator has recently been theoretically predicted. Here we explore the case where the two qubits are placed in different resonators in an array of two or three weakly coupled resonators. Investigating different setups and excitation schemes, we show that this process can still occur with a probability approaching one under specific conditions. The obtained results provide interesting insights into subtle causality issues underlying the simultaneous excitation processes of qubits placed in different resonators.

One of the most interesting nonlinear optical effects predicted in the USC regime consists of the simultaneously excitation of two or more spatially separated atoms by a single photon 23,[84][85][86] . This last puzzling result, which has been studied in a quantum system constituted by two qubit ultrastrongly coupled to a single-mode resonator, provides new insights into the various quantum aspects of the interaction between light and matter and can find useful applications for the development of novel quantum technologies. Although this effect clearly demonstrates a relevant role of the counter-rotating terms and virtual processes in the USC regime, the singlemode approach does not allow to fully understand some subtle causality issues underlying this process, since the resonator mode is completely delocalized along the cavity. Specifically, a drawback of this simplified description of the electromagnetic field is that any information about the spatial separation between the two atoms is lost. Hence, the question arises if it is possible to observe this effect in the presence of natural or artificial atoms which are actually spatially separated. Here we provide a positive answer to this question, even though further work will be required for a full understanding of the impact of the spatial separation of the atoms on the joint absorption and emission of single photons. Moreover, it has been recently pointed out 87 that the description of a cavity-QED system in terms of the single-mode quantum Rabi model in the USC and deep-strong coupling (DSC) regimes can lead to the violation of relativistic causality, and a multi-mode version of the quantum Rabi model is required in order to fully capture the propagation properties of the light field necessary to comply with causality. Here, although we do not adopt a multimode approach, we still assume the atoms to be spatially separated. Specifically, we introduce a simplified description of spatial separation, considering the two atoms embedded in different single-mode resonators. For example an array of N nearest-neighbour single-mode coupled resonators corresponds to a system with N distinct sites, like in tight-binding models in solid state physics, where propagation effects are taken into account. In these models the propagation speed depends on the interaction rate J between nearest neighbours cavities. For example the complete transfer of a photon from one resonator to the adjacent one requires a time T = π/(2J).
In the present work, we show that the simultaneous excitation process of two qubits by a single cavity photon, as described in Ref. 23 , can take place also in a cavity-array system of two or three cavities, where the qubits are placed in different resonators and ultrastrongly interact with them. We observe that this effect can be achieved by probing the system via two different excitation mechanisms: (1) by exciting one of the normal modes of the coupled-cavity array or (2) by selectively exciting only a single cavity. The substantial difference between the two cases is that, while the first case corresponds to the interaction of two qubits with a delocalized field and leads to a deterministic simultaneous excitation process as in Ref. 23 , the second case constitutes one of the simplest examples of a localized system in which the excitation is initially stored in a single resonator. This makes this effect even more counterintuitive, since the excitation is continuously transferred between the nearest-neighbor resonators, while the effect requires both atoms to feel the photon field at the same time in order to take place. In both cases, we study the temporal evolution of the system within both theoretical and numerical approaches, providing a clear and physically intuitive description of the propagation and causality mechanisms underlying this simultaneous excitation phenomena. The process described here could be experimentally realized in stateof-the-art circuit QED systems. Moreover, these effects can find useful applications in quantum information processing and quantum communication protocols, where reliable and controllable entanglement between distant qubits in a quantum network is of fundamental importance 88,89 .

Models and results
Here we study a quantum system consisting of an array of N weakly coupled single-mode resonators, each of them interacting with a two-level atom (e.g., a superconducting qubit). The total Hamiltonian of the system can be written as 23 (hereafter, = 1): − describe, respectively, the qubit and cavity Hamiltonians in the absence of interaction, â † n (â n ) is the bosonic creation (annihilation) operator for the nth resonator mode with frequency ω (n) c , and σ q . The cavity-cavity interaction Hamiltonian is given by where J is the next-neighbour hopping rate. We are assuming that the coupling strength J between the two cavities is weak, thus we can apply the RWA to the cavity-cavity interaction Hamiltonian Ĥ cc . Finally, the last term of Eq. (1), describing the interaction between the qubits and the cavity modes, reads.
where g n ≡ |g n |e iϕ n denotes the coupling rate of the nth qubit to the corresponding cavity field X n ≡â † n +â n , the angles θ n parametrize the relative contribution of the transverse and longitudinal couplings, while σ (n) x and σ (n) z are the Pauli matrices for the qubits. In circuit QED systems, the angle θ n and the transition frequency ω (n) q can be continuously tuned by changing the magnetic field externally applied to, e.g., a flux qubit (see, e.g., Ref. 43 ). An important feature of the interaction Hamiltonian is that it contains terms that do not conserve the total number of excitations. Specifically, the transverse coupling ∝σ x contains terms like âσ − and â †σ + which create or Scientific Reports | (2020) 10:21660 | https://doi.org/10.1038/s41598-020-78299-x www.nature.com/scientificreports/ annihilate two excitation simultaneously; whereas the σ z -coupling changes the resonator photon number by one, while leaving the number of qubit excitations unchanged. Notice that the parity of qubit n is conserved only for θ n = 0 , corresponding, for a flux qubit, to a zero external flux offset 43 . The terms in the total Hamiltonian which do not conserve the number of excitations in the system can be safely neglected in the weak-coupling regime, where the rotating-wave approximation (RWA) is valid. However, these terms become relevant for systems entering the USC regime, where the coupling strength g n reaches an appreciable fraction of the unperturbed frequencies (here, ω (n) c and ω (n) q ) of the bare systems. One of the most interesting consequences of the presence of these counter-rotating terms is the possibility to coherently couple quantum states with different numbers of excitations. These unconventional couplings determine new intriguing physical processes as, for example, multiphoton Rabi oscillations 21 and the possibility to excite two or more spatially separated atoms with a single photon 23 . Here, we investigate the latter process, considering a fundamentally different setup, consisting of weakly coupled cavity arrays, where each cavity ultrastrongly interacts with a single qubit. Our results show that, owing to the ultrastrong light-matter interaction and the parity-symmetry breaking, it is still possible to simultaneously excite two qubits, even if they are placed in different cavities.
Two cavities, two qubits and one photon. We first focus on the study of an array of two identical cavity-qubit systems ω where the two cavities are weakly coupled together and each one of them ultrastrongly interacts with a single qubit (see Fig. 1). In this case, the Hamiltonian of Eq. (1) becomes The Hamiltonian in Eq. (4) can be conveniently rewritten in terms of the bosonic symmetric and antisymmetric normal modes, also referred to as supermodes defined via the operators â S(A) = (â 1 ±â 2 )/ √ 2 , which diagonalise the Hamiltonian describing two weakly coupled harmonic oscillators.
In this case, we obtain The Hamiltonian in Eq. (6) describes two bosonic modes, one symmetric and one antisymmetric with corresponding frequencies ω S = ω c + J and ω A = ω c − J , both interacting with two qubits. We now diagonalise numerically Ĥ , indicating the resulting energy eigenvalues and eigenstates as ω i and |E i � , with i = 0, 1, . . . , . We label the states such that ω k > ω j for k > j . In our analysis, we use the notation |N A , N S , q 1 , q 2 � = |N A � |N S � |q 1 � |q 2 � for the eigenstates |E i � , where q = g, e denotes the qubit ground or excited states, respectively, and |N S(A) � = {|0�, |1�, |2�, . . . } represents the Fock states with a photon occupation number N in the symmetric (antisymmetric) normal mode.
We set the normalized hopping rate and the light-matter coupling strength as J/ω c = 0.05 and η ≡ |g|/ω c = 0.3 , respectively. Moreover, we set the phases for the cavity-qubit coupling strengths to ϕ 1 = 0 and ϕ 2 = π , respectively. This specific phase difference does not affect the dynamics, however it will play an important role for the case of three coupled cavities (see Sect.""Three cavities, two qubits and one photon). We also consider a mixing angle θ = π/6 , such that both the longitudinal and transverse contributions to the − + g e iϕ nX n cos θσ (n) x + sin θσ (n) Figure 1. Sketch of the system. Two identical spatially separated optical resonators, with resonance frequency ω c , are weakly coupled together, each one of them ultrastrongly interacting with a single two-level system (qubit) with transition frequency ω q . The photon hopping rate J between the two resonators and the light-matter coupling strength g are indicated.
Scientific Reports | (2020) 10:21660 | https://doi.org/10.1038/s41598-020-78299-x www.nature.com/scientificreports/ interaction have comparable values. Figure 2a shows the energy differences ω i0 = ω i − ω 0 for the lowest-energy states as a function of the normalized qubit frequency ω q /ω c , which can be experimentally tuned by changing the external magnetic flux acting on the qubit. We observe that, when ω q ≃ ω A /2 , the spectrum displays an avoided-level crossing between the states |E 3 � and |E 4 � . It is worth noticing that the resonance condition is quite different from the expected one, i.e., ω q = ω A /2 , because ω q and ω A are bare resonance frequencies of the matter and light components. The actual physical frequencies are significantly dressed by the interaction (see, e.g., 21,23,25,79,90 ). Note that, just outside this avoided-level crossing region, one level remains flat as a function of the qubit frequency with energy ω ≈ ω A , while the other shows a linear behaviour with ω ≈ 2ω q . The origin of this splitting is due to the hybridization of the states |1 A , 0 S , g, g� and |0 A , 0 S , e, e� . When the splitting is at its minimum, these states are well approximated by the superposition states The (numerically calculated) normalized minimum splitting has a value 2� eff /ω c = 16 × 10 −3 , where eff is the effective coupling rate between two qubits and one photon. It is important to observe that the coherent coupling between these two states would not be allowed within the RWA, since they have a different number of excitations. Moreover, as the excitation number difference between the two states is odd, parity-symmetry breaking ( θ = 0 ) is required in order to observe this splitting. As reported in Ref. 23 , the effective coupling between the states |E 3 � and |E 4 � can be analytically described by an effective Hamiltonian. Moreover, this coherent coupling is not direct, but can only occur via virtual transitions which are enabled by the counter-rotating terms in Ĥ . In this way, the initial state |1 A , 0 S , g, g� evolves to virtual intermediate states that eventually do not conserve the energy, but it finally evolves to a real energy-conserving state, i.e., |0 A , 0 S , e, e� . It is interesting to observe that, for ω q ≃ ω S /2 , the energy spectrum displays a crossing between the levels |0 A , 1 S , g, g� and |0 A , 0 S , e, e� , showing that the two qubits do not interact with the symmetric normal mode of the coupled cavities. Indeed, if the coupling strengths g 1 and g 2 have opposite signs it can be shown that all the possible intermediate virtual transitions for the process (a) Energy differences ω i0 = ω i − ω 0 for the lowest-energy dressed states of Ĥ as a function of the normalized qubit frequency ω q /ω c (which can be experimentally tuned by changing the external flux bias acting on the qubits). We consider a normalized coupling rate η ≡ |g|/ω c = 0.3 between the qubit and the resonators, while the normalized photon hopping rate between the two resonators is J/ω c = 0.05 . The phases for the cavity-qubit coupling strengths are set to ϕ 1 = 0 and ϕ 2 = π , respectively, and the longitudinal interaction coupling term is included by considering a mixing angle θ = π/6 . (b) Enlarged view of the inset in (a). When the cavity-qubit coupling strengths have opposite phases ( ϕ 1 = 0, ϕ 2 = π ), the avoided-level crossing (red solid curves) results from the coupling between the states |1 A , 0 S , g, g� and |0 A , 0 S , e, e� due to the presence of counter-rotating terms in the system Hamiltonian. The energy splitting reaches its minimum at ω q ≃ ω A /2 . The crossing between the states |0 A , 1 S , g, g� and |0 A , 0 S , e, e� at ω q ≃ ω S /2 indicates that the qubits do not couple with the symmetric normal mode. The complementary result is obtained when the coupling strengths have the same phase (blue dashed curves). In this case, this splitting at ω q ≃ ω A /2 disappears, while the energy spectrum displays an avoided-level crossing around ω q ≃ ω S /2 , arising from the coherent coupling between the states |0 A , 1 S , g, g� and |0 A , 0 S , e, e�.
Scientific Reports | (2020) 10:21660 | https://doi.org/10.1038/s41598-020-78299-x www.nature.com/scientificreports/ |0 A , 1 S , g, g� → |0 A , 0 S , e, e� , induced by the interaction term proportional to X S , lead to an intermediate state proportional to |e��e| − |e��e| , thus giving a vanishing contribution. Figure 2b shows the comparison of the avoided-level crossing behavior for the cases ϕ 1 = 0, ϕ 2 = π (solid red curves) and ϕ 1 = ϕ 2 = 0 (dashed blue curves). It can be observed that the two choices lead to complementary results. Indeed, when the two coupling strengths have same signs, the splitting at ω q ≃ ω A /2 disappears, while we observe the presence of an avoided-level crossing around ω q ≃ ω S /2 arising from the coherent coupling between the states |0 A , 1 S , g, g� and |0 A , 0 S , e, e� . This result shows that, depending on the relative signs of the coupling strengths, for N = 2 the simultaneous excitation of two qubits placed in different resonators can be achieved by coupling the qubits either to the symmetric or antisymmetric normal mode.
In order to fully understand and characterise this process, we fix the qubit frequency at the value where the splitting, as shown in Fig. 2a, between the energy levels corresponding to the eigenstates |E 3 � and |E 4 � is minimum and consider the system initially prepared in the one-photon state |1 A � ≡ |1 A , 0 S , g, g� . As we will see later, the preparation of this state can be experimentally achieved by sending an appropriate electromagnetic Gaussian pulse to the first cavity. Figure 3a displays the numerically calculated time evolution of the occupation probabilities P (k) (t) ≡ �P k � , with P k = |k��k| , for the one-photon states |1 A � (i.e., one photon in the antisymmetric normal mode) and |1 1 � (one photon in the first cavity), together with the probability P (ee) (t) of having both qubits simultaneously excited. In terms of the dressed energy eigenstates of the system, these states can be expressed, respectively, as Figure 3. (a) Time evolution of the occupation probabilities P (k) (t) ≡ �P k � , with P k = |k��k| , for the singlephoton states |1 A � (red solid curve) and |1 1 � (blue dot-dashed curve), together with the probability P (ee) (t) of having both qubits simultaneously excited (black dashed curve) for the system initially prepared in the state Vacuum Rabi oscillations showing a reversible excitation exchange process between the qubits and the resonators are clearly visible. The joint absorption of an antisymmetric cavity photon by the two qubits is achieved after a Rabi half period � eff t = π/2 , with the excitation probability P (ee) approaching one, even if they are placed in different resonators. Here, the effects of dissipation have not been included.
(b) Temporal evolution of the same occupation probabilities P (k) (t) considered in (a), but after the arrival of a narrow Gaussian pulse exciting the first cavity when the system is initially prepared in its ground state. The amplitude and the central frequency of the pulse are A/ω c = 4.2 × 10 −2 and ω d = (ω 30 + ω 40 )/2 , respectively.
Scientific Reports | (2020) 10:21660 | https://doi.org/10.1038/s41598-020-78299-x www.nature.com/scientificreports/ where for concision, we omitted in the last equation the photonic states, which are intended to be in the ground state. Here, |e, e� stands for |0 A , 0 S , e, e� . As expected, since |1 A � = (|1 1 , 0 2 � − |0 1 , 1 2 �)|g, g�/ √ 2 , at the initial instant of time P (1 A ) (0) = 1 , while P (1 1 ) (0) = 1/2 . As time evolves, vacuum Rabi oscillations, showing the reversible excitation exchange process between the qubits and the resonators, are clearly visible. Specifically, we observe that, after a Rabi half period � eff t = π/2 , one photon in the antisymmetric cavity mode is jointly absorbed by the two qubits, even if they are placed in different resonators. Moreover, the excitation probability P (ee) approaches one, showing that the multiatom absorption of a single photon can essentially be deterministic. Notice that, in order to provide a clearer description of this counter-intuitive excitation mechanism, the effects of dissipation have not been taken into account. This approximation becomes experimentally reasonable when the system loss rates are smaller than the frequency splitting between the levels involved at the avoided-level crossing, so that the first Rabi cycles are almost not affected by dissipation. A similar oscillating dynamics can be obtained for the system initially prepared in the state |e, e� . In this case, the two qubits will jointly and coherently release their energy to the cavity. The time evolution of the system will be as shown in Fig. 3a, but with the initial time t = π/(2� eff ).
As mentioned before, instead of starting from the ideal initial state |1 A � , we now consider a more realistic case where the system is initially in its ground state |E 0 � = |0 A , 0 S , g, g� and study a direct excitation of the first cavity by an electromagnetic Gaussian pulse. The corresponding driving Hamiltonian is where X 1 =X S +X A and with A and τ the amplitude and the standard deviation of the Gaussian pulse, respectively. The central frequency of the pulse has been chosen to be in the middle of the two split transition energies ω d = (ω 30 + ω 40 )/2 . The pulse bandwidth must be sufficiently narrow in order to ensure that only the states |E 3 � and |E 4 � are excited, so that the pulse can directly excite the state |1 A � = (|E 3 � + |E 4 �)/ √ 2 , and the symmetric mode |1 S � is not excited. This corresponds to a pulse duration ∼ τ significantly larger than the transfer time ∼ π/(2J) of photons from one cavity to the other. In this way, even if the system is fed through a single cavity only, both cavities are actually excited simultaneously without any causality issue, because the excitation time τ is larger than the transfer time ∼ π/(2J) . With this excitation scheme, the resulting dynamics is very similar to the case of two atoms in a single cavity 23 . Figure 3b shows the dynamics of the occupation probabilities P (1 A ) , P (1 1 ) , and P (ee) after the arrival of the π -like Gaussian pulse initially exciting the first cavity described by the Hamiltonian in Eq. (11). We observe that, since the pulse time width is not much narrower than the Rabi period, after the arrival of the pulse, the antisymmetric normal mode is not completely populated and the excitation is partially transferred to the qubits. Therefore, the first peak of the antisymmetric normal mode occupation probability in Fig. 3b is slightly lower than the second one. Once the antisymmetric mode is completely populated, the dynamics of vacuum Rabi oscillations, showing the reversible excitation exchange between one photon in the antisymmetric normal mode and the two qubits, is the same as in Fig. 3a. It is also worth noticing that, in the absence of any system nonlinearity as, e.g., in the case of a system empty (without atoms) resonators, coherent excitation, as described by the Hamiltonian in Eq. (11), would give rise to a coherent intra-cavity field, not to a single-photon Fock state. However, the very strong interaction of the cavity array with the two atoms induces an anharmonicity to the level structure which is able to prevent the resonant excitation of higher-photon states (photon-blockade, see, e.g., 91,92 ). In the case of lower atom-cavity coupling strengths, the photonic system could be complemented by additional Kerr nonlinearities, which are able to induce photon blockade 23 .
The simultaneous excitation process of the two qubits with a single photon can also be achieved by initially exciting only the first cavity. This can be realized experimentally by feeding the cavity with a fast Gaussian pulse (with respect to the transfer time between the cavities). Specifically, in order to completely populate the first cavity before the excitation is transferred to the second one, the pulse bandwidth Ŵ = π/τ has to be much larger than the energy splitting = 2J between the two cavity normal modes, i.e., Ŵ ≫ � . Figure 4a shows the time evolution of the occupation probabilities P (1 1 ) (t) , P (1 2 ) (t) , and P (ee) (t) with the system initially prepared in the state Besides the expected photon hopping between the two weakly-coupled cavities, as the time evolves, we observe that this excitation-transfer process is accompanied by the simultaneous excitation of the two qubits. However, in contrast to the previous case, where the probability for the qubits to be simultaneously excited reached one, here we observe the maximum probability P (ee) = 1/2 at t = π/(2� eff ) . This difference can be explained by considering that a direct excitation of the first cavity only corresponds to an equal weight superposition of both the symmetric and antisymmetric normal modes (see Eq. 13), each of them carrying half of the total initial excitation. However, since the qubits do not interact with the symmetric mode, only the antisymmetric excitation contribution can be transferred to the two qubits, resulting into a maximum joint qubit excitation probability P (ee) = 1/2 . This result indicates that the simultaneous excitation of two qubits with one photon is still possible, but the qubits cannot be excited with probability P (ee) = 1 . Moreover, the corresponding values P (1 1 ) = P (1 2 ) = 1/4 for the cavity occupation probabilities can be explained by directly following the dynamics of the system. Indeed we observe that, when the system is prepared in the superposition state |ψ 1 � , then the state |0 A , 1 S , g, g� evolves freely, as an eigenstate of Ĥ . On the contrary, since the qubits are coupled with the antisymmetric normal mode, once Scientific Reports | (2020) 10:21660 | https://doi.org/10.1038/s41598-020-78299-x www.nature.com/scientificreports/ again we observe the coherent energy exchange process |1 A , 0 S , g, g� ↔ |0 A , 0 S , e, e� so that at t = π/(2� eff ) the system is in the state In terms of the energy eigenstates of the Hamiltonian of the uncoupled system ( J = g = 0 ), the state |ψ 2 � can be expressed as thus explaining the observed values, in Fig. 4a, for the occupation probabilities P (1 1 ) , P (1 2 ) and P (ee) . Finally, even for this case we consider the system initially prepared in its ground state and study the system dynamics after the arrival of a Gaussian pulse exciting the first cavity and described by Eq. (11). Unlike the previous case, in order to excite the states |E 3 � , |E 4 � and |E 5 � we need to apply a broad-bandwidth Gaussian pulse with central frequency ω d ≃ (ω 34 + ω 50 )/2 , where ω 34 = (ω 30 + ω 40 )/2 . The temporal evolution of the occupation probabilities, after the arrival of the Gaussian pulse feeding the first cavity, is shown in Fig. 4b. It is interesting to observe that, in the absence of the qubits, the excitation would simply be transferred from one cavity to the another. Indeed, the simultaneous excitation process can take place because, except for some specific instants of time, in which the excitation is totally localized in one cavity only, the field is delocalized over the two adjacent cavities. For this reason, both qubits feel the electric field simultaneously reaching the maximum excitation probability when the field is equally distributed between the two cavities, even if they are detuned from the cavity mode.
The causal mechanisms underlying the possibility for a single photon to be jointly absorbed by the two qubits under different excitation processes can be explained by considering a simpler system constituted by the two weakly coupled cavities only. If the system is initially prepared in the state |1 1 , 0 2 � , the complete population transfer |1 1 , 0 2 � → |0 1 , 1 2 � will take place after a period T = π/(2J) . From an experimental point of view, if we consider the system to be prepared in its ground state |0 1 , 0 2 � , the direct excitation of the antisymmetric normal (14) |ψ 2 � = |0 A , 0 S , e, e� + |0 A , 1 S , g, g� / √ 2.
In contrast to this, the realization of a localized cavity mode (e.g., the state |1 1 , 0 2 � ) can be achieved by feeding the first cavity with a fast optical Gaussian pulse whose bandwidth has to be large enough to excite the superposition state (|+� + |−�)/ √ 2 . The condition Ŵ ≫ 2J ensures that, being δt ≪ T , the pulse duration is short enough to localise the electric field in the first cavity. Then, due to the cavity-cavity interaction the excitationtransfer process will take place and the electric field will be completely delocalized over the two cavities only at the instants of time when P (1 1 ) = P (1 2 ) = 1/2.

Three cavities, two qubits and one photon.
Here we extend the previous analysis to a more complex system consisting of an array of three weakly-coupled cavities, where the two end cavities ultrastrongly interact with a single qubit while the central cavity is empty (as shown in Fig. 5a). The Hamiltonian describing the system is where the normalized hopping rate and light-matter coupling strength are, respectively, J/ω c = 0.05 and η ≡ |g|/ω c = 0.3 , and the presence of the longitudinal interaction term is taken into account by considering a mixing angle θ = π/6 . Due to the small value of the normalized hopping rate J/ω c , we apply the RWA to the cavity-cavity interaction term. Moreover, we consider the case in which the two end cavities are resonant (ω (1) , while the central cavity can be detuned by an amount . The interaction between these three cavities is described by the Hamiltonian www.nature.com/scientificreports/ which produces three normal modes: an antisymmetric flat state, and two symmetric modes whose energy splitting is = 8J 2 + 2 . The transformation, which diagonalizes Ĥ ′ C , can be written in matrix form: with N ± = 8J 2 + (� ± �) 2 1/2 . In Eq. (18), â A and â S 1 (S 2 ) are, respectively, the bosonic annihilation operators for the antisymmetric state and the two symmetric normal modes whose corresponding frequencies are ω A = ω c and ω S 1,2 = (2ω c + � ∓ �)/2. The effect of the detuning on the spatial profile of the three normal modes is displayed in Fig. 5b. We observe that, unlike the antisymmetric mode, which is not affected by the detuning, the spatial profile of the two symmetric normal modes changes significantly with increasing values of , and also exhibits an opposite behavior. Specifically, while the lower-energy mode â S 1 becomes totally delocalized in the end cavities, the higher-energy mode â S 2 completely localizes in the central cavity. This preliminary analysis suggests that the choice of coupling the two qubits to the antisymmetric mode, which is both detuning-independent and delocalized in the two end cavities where the qubits are placed, is the best strategy for achieving the desired effect of simultaneous excitation of two qubits with a single photon. Focusing on the resonant case = 0 , the Hamiltonian of Eq. (16) can be conveniently rewritten in terms of the normal mode operators as Figure 6a shows the energy differences ω i0 = ω i − ω 0 for the lowest-energy eigenstates |E i � (with i = 0, 1, . . . ) of Ĥ ′ , numerically calculated as a function of the normalized qubit frequency ω q /ω c by setting the phases of the two cavity-qubit coupling strengths to ϕ 1 = 0 and ϕ 3 = π , respectively. Here we use the notation |N S 1 , N S 2 , N A , q 1 , q 2 � = |N S 1 � |N S 2 � |N A � |q 1 � |q 2 � for the eigenstates |E i � , where q = g, e denote the qubit ground or excited states, respectively, and |N k � = {|0�, |1�, |2�, . . . } , with k ∈ [S 1 , S 2 , A] , represents the Fock state with photon occupation N k in the corresponding normal mode.
Note that the energy spectrum presents a more complicated structure with respect to the case of the two cavity-qubit array. Specifically, two energy-level crossings, both involving the state |0 S 1 , 0 S 2 , 0 A , e, e� (which displays a linear behavior with ω ≈ 2ω q ), can be observed in correspondence to the qubit frequencies ω q ≃ ω S 1 /2 and ω q ≃ ω S 2 /2 . The other states involved in the two energy-level crossings are, respectively, |1 S 1 , 0 S 2 , 0 A , g, g� and |0 S 1 , 1 S 2 , 0 A , g, g� , indicating that when the cavity-qubit coupling strengths have opposite phases, the qubits do not couple with the two symmetric normal modes. Interestingly, in the region between these two-level crossings, an apparent additional one between the levels |E 4 � and |E 5 � appears at ω q ≃ ω A /2 . Actually, what appears as a crossing on this scale turns out to be an avoided-level crossing on an enlarged view, as in Fig. 6b. This splitting, which has a normalized value 2� eff /ω c = 2 × 10 −3 at its minimum, clearly originates from the hybridization of the states |0 S 1 , 0 S 2 , 0 A , e, e� and |0 S 1 , 0 S 2 , 1 A , g, g� . The resulting states are well approximated by It is important to observe that, similarly to the two cavity-qubit array case, the coherent coupling between these two states would neither be allowed within the RWA, nor in the absence of the longitudinal interaction term (θ = 0) . Moreover, the states |0 S 1 , 0 S 2 , 0 A , e, e� and |0 S 1 , 0 S 2 , 1 A , g, g� do not couple directly, but the process occurs via intermediate energy non-conserving processes enabled by the counter-rotating terms in Ĥ ′ . It is interesting to observe that, when the coupling strengths are opposite in phase, all the possible intermediate virtual transitions for the process |0 S 1 , 0 S 2 , 1 A , g, g� → |0 S 1 , 0 S 2 , 0 A , e, e� are induced by the interaction term proportional to X A , while the term proportional to X S 2 −X S 1 gives vanishing contributions. The choice of same coupling strength phases (ϕ 1 = ϕ 3 ) would lead to the complementary situation with the two qubits decoupled from the antisymmetric mode and simultaneously interacting with the two symmetric modes â S 1 and â S 2 . In this case, the simultaneous excitation of the two qubits with a single photon can occur via two different processes ( |1 S 1 , 0 S 2 , 0 A , g, g� → |0 S 1 , 0 S 2 , 0 A , e, e� and |0 S 1 , 1 S 2 , 0 A , g, g� → |0 S 1 , 0 S 2 , 0 A , e, e� ). However, since the symmetric normal modes are mainly localized in the central empty cavity, the energy splittings of the avoided-level crossings between these states, as well as the corresponding effective couplings, are much smaller.
Dynamics after exciting the antisymmetric mode. In order to fully understand the excitation transfer between a single photon and two qubits in a three cavity-qubit array, we now study the dynamics of the system initially prepared in the one-photon state |1 A � ≡ |0 S 1 , 0 S 2 , 1 A , g, g� , fixing the qubit frequency at the value where the splitting in Fig. 6b between levels |E 4 � and |E 5 � is minimum. Moreover, for the sake of simplicity, we consider the system loss rates to be significantly smaller than the frequency splitting between the levels involved in the avoided-level crossing, so that the effect of dissipation can be neglected.
Scientific Reports | (2020) 10:21660 | https://doi.org/10.1038/s41598-020-78299-x www.nature.com/scientificreports/ photon in the first and third cavities, respectively), together with the probability P (ee) (t) of having both qubits simultaneously excited, are displayed in Fig. 7a. As expected, P (1 A ) (0) = 1 at the initial instant of time, and since |1 A � = (|0 1 , 0 2 , 1 3 � − |1 1 , 0 2 , 0 3 �)/ √ 2 , we observe that the end cavities are equally populated P (1 1 ) (0) = P (1 3 ) (0) = 1/2 , while the central cavity is empty P (1 2 ) (0) = 0 . As time evolves, the excitation is progressively transferred to the two qubits at the same time, until the maximum simultaneous qubit excitation [with P (ee) (0) = 1 ] is reached at t = π/(2� eff ) . It is interesting to observe that, since we are exciting the antisymmetric mode, the central cavity remains empty during the whole process. This fact remains valid even if the central cavity is detuned (� = 0) , so that the system displays the same dynamics in the non-zero detuning case, the only difference being that the simultaneous excitation of the two qubits is achieved at a different instant of time.
The excitation of the antisymmetric normal mode, which is the most effective way to achieve the desired effect, can be experimentally achieved by sending a suitable narrow Gaussian pulse to the first cavity, whose central frequency of the pulse has to be chosen to be in the middle of the two split transition energies ω d = (ω 40 + ω 50 )/2 . For the sake of simplicity, here we do not present numerical calculations for the dynamics of the system excited by a Gaussian pulse, since the results would not add any additional physical information.
Dynamics after exciting the first cavity. We now turn to the study of the system dynamics when, instead of exciting the antisymmetric mode, we directly excite only the first cavity. Unlike the previous case, this process strongly depends on the detuning .
The dynamics of the occupation probabilities P (k) for a three coupled cavity-qubit array system initially prepared in the one-photon state |1 1 � can be studied for arbitrary values of the detuning by using a semi analytical approach. Choosing the energy eigenstates |E j � of Ĥ ′ as basis, and considering that only the states |E 3 �, |E 4 �, |E 5 � ,  Figure 6. Energy differences ω i0 = ω i − ω 0 for the lowest-energy dressed states of H ′ as a function of the normalized qubit frequency ω q /ω c in the absence of detuning ( = 0 ). We set the normalized qubit-resonator coupling rate and the inter-cavity photon hopping rate to η ≡ |g|/ω c = 0.3 and J/ω c = 0.05 , respectively. The phases for the cavity-qubit coupling strengths are ϕ 1 = 0 and ϕ 3 = π , while the longitudinal interaction coupling term is included by considering a mixing angle θ = π/6 . (b) Enlarged view of the inset in (a). When the cavity-qubit coupling strengths are opposite in phase ( ϕ 1 = 0, ϕ 3 = π ), the avoided level crossing results from the coupling between the states |0 S 1 , 0 S 2 , 1 A , g, g� and |0 S 1 , 0 S 2 , 0 A , e, e� , due to the presence of counterrotating terms in the system Hamiltonian. The energy splitting reaches its minimum at ω q ≃ ω A /2 . The presence of two energy-level crossings corresponding to the qubit frequencies ω q ≃ ω S 1 /2 and ω q ≃ ω S 2 /2 , both involving the state |0 S 1 , 0 S 2 , 0 A , e, e� , indicates that the two qubits do not couple with the two symmetric normal modes.
Scientific Reports | (2020) 10:21660 | https://doi.org/10.1038/s41598-020-78299-x www.nature.com/scientificreports/ and |E 6 � are involved in the process, the time evolution of the initial state |1 1 (t)� in the Schrödinger picture can be written as: where ω j are the numerically evaluated eigenvalues and the coefficients c j , which depend on J and , are given by the elements of the transformation matrix in Eq. (18). The time evolution of the occupation probability P (k) (t) for a generic state |k� = 6 j=3 d j |E j � are simply given by:  Figure 7. (a) Time evolution of the occupation probabilities P (k) (t) ≡ �P k � , with P k = |k��k| , for the onephoton states |1 A � (red solid curve) and |1 1 � (blue dot-dashed curve), together with the probability P (ee) (t) of having both qubits simultaneously excited (black dashed curve). Here, the system is initially prepared in the one-photon state |1 A � ≡ |0 S 1 , 0 S 2 , 1 A , g, g� and in the absence of detuning ( = 0 ). As the times evolves, the excitation is progressively transferred to the two qubits at the same time, until the maximum simultaneous qubit excitation ( P (ee) = 1 ) by a single photon is reached at t = π/(2� eff ) . (b) Temporal dynamics of the occupation probabilities P (1 1 ) (blue solid curve), P (1 2 ) (green dashed curve) and P (1 3 ) (red dot-dashed curve), together with the probability P (ee) (black dashed curve) of having both qubits simultaneously excited when the system is initially prepared in the one-photon state |ψ 1 � = |1 1 � of the first resonator with = 0 . (c) Temporal dynamics of the same occupation probabilities P (k) (t) considered in (b) but in the presence of the detuning �/ω c = 0.5.
Scientific Reports | (2020) 10:21660 | https://doi.org/10.1038/s41598-020-78299-x www.nature.com/scientificreports/ Figure 7b shows the time evolution of the occupation probabilities P (1 1 ) (t) , P (1 2 ) (t) , P (1 3 ) (t) , and P (ee) (t) for the resonant case = 0 , when the system is initially prepared in the one-photon state As the time evolves, we observe that the excitation, continuously propagating between the three cavities, is progressively transferred to both qubits at the same time, even if they are detuned from the central cavity. We observe that initially the system dynamics is not affected by the presence of the qubits and the system behaves like an array of three weakly-coupled resonant cavities. Once the excitation starts to be transferred to the qubits, the system undergoes a more complex dynamics and the probability for the qubits to be simultaneously excited reaches its maximum P (ee) = 1/2 at t = π/(2� eff ) . As expected, at this instant of time the electric field is completely delocalized only in the end cavities with P (1 1 ) = P (1 3 ) = 1/4 , while the central cavity is depopulated. This result can be easily understood by considering that, after the system is prepared in the state |ψ 1 � , the superposition |0 S 1 , 1 S 2 , 0 A , g, g� − |1 S 1 , 0 S 2 , 0 A , g, g� /2 evolves freely since the qubits are not coupled to the symmetric modes.
In contrast to this, due to the interaction between the qubits and the antisymmetric mode, the above-described coherent-energy-exchange process |0 S 1 , 0 S 2 , 1 A , g, g� → |0 S 1 , 0 S 2 , 0 A , e, e� takes place so that at t = π/(2� eff ) the system will be in the state The observed values for the occupation probabilities P (k) in Fig. 4b can be explained by considering that, in terms of the energy eigenstates of the Hamiltonian of the uncoupled system (J = g = 0) , this state can be expressed as Finally, in Fig. 7c we present numerical results for the dynamics of the system initially prepared in the state |ψ 1 � for �/ω c = 0.5 . In this case, the excitation of the first cavity can be obtained by exciting a superposition of the lowest-energy symmetric mode and the antisymmetric mode (see Fig. 5b). This could be experimentally realized by sending a suitable broad Gaussian pulse to the first cavity, able to excite the energy levels |E 3 � , |E 4 � , |E 5 � and |E 6 � simultaneously.
We observe that the system dynamics displays a different trend with respect to the resonant case. Indeed, due to the strong detuning, the central cavity acts like a high-potential barrier and the excitation is transferred back and forth via photon tunneling only between the end cavities, with the system effectively behaving like a two coupled cavity-qubit array. During the whole process, the central cavity remains very low-populated and the excitation is progressively transferred to the qubits, and the maximum simultaneous excitation probability P (ee) = 1/2 is reached at t = π/(2� eff ) , where � eff /ω c = 4.5 × 10 −4 is the effective coupling for �/ω c = 0.5 between two qubits and a single photon.
The processes described here could be experimentally observed by placing two superconducting artificial atoms at opposite ends of an array of capacitively-coupled superconducting waveguides. These anomalous multiatom excitation and emission processes can find applications for the development of novel quantum technologies for quantum information and communication as, for example, the realization of new effective methods for quantum information transfer between photons and qubits in quantum networks.

Conclusions
When the light-matter coupling strength increases, the vacuum fluctuations of the electromagnetic field become able to efficiently induce virtual transitions, replacing the role of the intense applied fields in nonlinear optics 10,11,81 . In this way, higher-order processes involving counter-rotating terms can create an effective coupling between two states of a system with different numbers of excitations. One of the most interesting examples is the process where a single photon in an electromagnetic resonator can jointly excite two atoms interacting with the same resonator. We have investigated this intriguing nonlinear optical process in the case where each of the two atoms is coupled to a distinct resonator. Specifically we studied: (1) the case of two coupled resonators, each of them interacting with a single atom, and (2) the case of two resonator-atom systems weakly coupled through a central resonator (resonant or detuned with respect to the other two). We studied the dynamics of these systems under different excitation conditions, showing that a coherent energy transfer between a single photon and two spatially separated atoms can still occur with a probability approaching one, under specific excitation conditions.
The results obtained provide interesting insights into subtle causality issues underlying the simultaneous excitation of two-level systems placed in spatially separated resonators. The processes described here could be experimentally realized in state-of-the-art circuit QED systems. Specifically, in our calculations we used a normalized coupling strength η = |g|/ω c = 0.3 . It is not a problem to experimentally reach such a value in circuit QED systems. Using a flux qubit coupled to an LC oscillator, values of η of the order of ≃ 0.65 , and even beyond it, have been obtained. In order to observe the dynamics shown, e.g., in Fig. 3, the broadenings should be significantly lower than the splitting shown in Fig. 2b. Such a splitting is of the order �/ω c ∼ 10 −3 . For a typical resonance frequency ω c /2π ≃ 6 GHz, the splitting is thus of the order of 0.35 GHz. In circuit QED systems, typical damping rates are orders of magnitude lower than this value 5 . Notice also that such a splitting can be increased by increasing η ; it scales as η 323 . Recently, the possibility to coherently couple distant ( ∼ 1 m) cavity-QED systems, (23) |ψ 1 � = |0 S 1 , 1 S 2 , 0 A , g, g� − |1 S 1 , 0 S 2 , 0 A , g, g� /2 − |0 S 1 , 0 S 2 , 1 A , g, g�/ √ 2 = |1 1 �.
Scientific Reports | (2020) 10:21660 | https://doi.org/10.1038/s41598-020-78299-x www.nature.com/scientificreports/ connecting them by a microwave cable with negligible loss 88,89 has been experimentally demonstrated. These results opens the possibility to implement the process proposed here involving distant artificial atoms. It would be interesting and useful for applications to extend these studies to the cases of spatially separated atomic ensembles 86 , and distant qubits coupled to open waveguides 93 . Moreover, further insights on the processes here analysed can be gained by considering two or more spatially separated atoms in a multimode cavity 87 . It would also be interesting to extend the description of other vacuum-boosted nonlinear optical processes in the USC regime 25,81 , including coupled and/or multi-mode resonators.