Dynamics of mesoscopic qubit ensemble coupled to cavity: role of collective dark states
A. A. Zhukov D. S. Shapiro W. V. Pogosov Yu. E. Lozovik





^1N. L. Dukhov All-Russia Research Institute of Automatics, 127055 Moscow, Russia
^2National Research Nuclear University (MEPhI), 115409 Moscow, Russia
^3V. A. Kotel'nikov Institute of Radio Engineering and Electronics, Russian Academy of Sciences, 125009 Moscow, Russia
^4Institute for Theoretical and Applied Electrodynamics, Russian Academy of
Sciences, 125412 Moscow, Russia
^5Institute of Spectroscopy, Russian Academy of Sciences, 142190 Moscow region,
Troitsk, Russia
^6Moscow Institute of Physics and Technology, 141700 Moscow region, Dolgoprudny, Russia


We consider dynamics of a disordered ensemble of qubits interacting with single mode photon field, which is described by exactly solvable inhomogeneous Dicke model. In particular, we concentrate on the crossover from few-qubit systems to the system of many qubits and analyze how collective behavior of coupled qubits-cavity system emerges despite of the broadening. We show that quantum interference effects survive in the mesoscopic regime – dynamics of an entangled Bell state encoded into the qubit subsystem remains highly sensitive to the symmetry of the total wave function. Moreover, relaxation of these states is slowed down due to the formation of collective dark states weakly coupled to light. Dark states also significantly influence dynamics of excitations of photon subsystem by absorbing them into the qubit subsystem and releasing quasiperiodically in time. We argue that predicted phenomena can be useful in quantum technologies based on superconducting qubits. For instance, they provide tools to deeply probe both collective and quantum properties of such artificial macroscopic systems.


02.30Ik, 42.50.Ct, 03.65.Fd


Dynamics of mesoscopic qubit ensemble coupled to cavity: role of collective dark states
    
    December 30, 2023
=======================================================================================





§ INTRODUCTION








Controllable manipulation by quantum states of spin-photon coupled systems attracts a considerable current interest, since ensembles of spins/atoms interacting with the quantized electromagnetic field are promising candidates for implementation of quantum information and computation devices <cit.>. One of the most perspective applications of such systems is storage of quantum information <cit.>. There are various physical realizations of spin-photon coupled systems, which range from superconducting artificial atoms (qubits) coupled to microwave resonators to nitrogen-vacancy (NV) centers in diamond and include even hybrid circuits combining two or more physical systems <cit.>. Effective parameters of state-of-art spin-photon systems can be very different as well as numbers of coherently interacting spins (qubits). For instance, a typical number of NV centers in the ensemble is macroscopically large, while the coupling of a single center to the cavity mode is very weak. In contrast, state-of-art superconducting quantum circuits are limited by tens of qubits, while interaction between a single qubit and the microwave radiation can be relatively strong <cit.>. Nevertheless, essentially any solid-state physical realization is characterized by the inhomogeneous broadening of the density of states, i.e., by the splitting between excitation frequencies of individual spins. Broadening is caused by fundamental mechanisms and therefore is poorly controlled: for example, an excitation frequency of superconducting flux qubits depends exponentially on Josephson energies of contacts embedded into their structure <cit.>, which makes it highly sensitive to characteristics of nanometer-scale Josephson junctions. For NV-centers, inhomogeneous broadening is induced by background disorder <cit.>.

Disorder in spin excitation energies is usually considered as a negative phenomenon, which prevents quantum information processing by introducing a decoherence <cit.>. In order to overcome this problem, various ideas for the spectral engineering of spin density profile <cit.> or spectral hole burning <cit.> have been proposed (mostly in the context of NV centers in diamond). However, inhomogeneous broadening can also play a positive role since it may be utilized for the construction of multi-mode quantum memories <cit.>. Theoretical description of such systems is usually provided within Tavis-Cummings (Dicke) model, while interpretations are developed in terms of so-called radiant and subradiant modes, which turn out to be coupled if broadening is present in the system, since in this case they do not match exact eigenstates of the Dicke Hamiltonian <cit.>. In other words, the excitation stored initially in the radiant mode is finally absorbed by a bath of subradiant states. A smart idea aimed to circumvent this problem is based on the hole burning technique: spins within certain 'dangerous' energy intervals, which are predominantly responsible for the interaction with subradiant modes, are neutralized for some time period by an external pulse <cit.>. As a result,  special light-matter quantum states can be engineered, these states being mostly localized within the spin subsystem. If energy dissipation in a cavity is much larger than in the spin subsystem, as usually valid for realizations based on natural quantum systems, by using this approach it is possible to efficiently protect quantum state from the cavity decay, which might lead to realizations of quantum memory prototypes.

The aim of the present paper is a general study of dynamics of inhomogeneously broadened spin ensembles of mesoscopic rather than macroscopic sizes. This is especially actual for possible realizations of such ensembles coupled to microwave resonators within superconducting platform (sometimes referred to as quantum metamaterials, see, e.g., Refs. <cit.>) and perhaps some other future solid-state circuits. In particular, we focus on the crossover from the system of just few qubits to the macroscopic system and study how collective dynamical behavior emerges along this crossover. This is done using an exact solution of Dicke Hamiltonian via Bethe ansatz <cit.>. We restrict ourselves to the regime of weak excitation, i.e., when there is no more that one excitation in the system. We also consider different initial conditions and show that they can result in qualitatively different dynamics. Our approach provides a simple and pictorial understanding for main features of system's evolution. It allows to obtain a direct access to Hamiltonian eigenstates, which can be classified as dark and bright, and their properties, as well as to study explicitly role of such states in the system's dynamics. In contrast to earlier studies <cit.>, we mostly concentrate on the dynamics starting from excitations within spin subsystem being motivated by recent experimental advances in hybrid systems <cit.> and consider mesoscopic qubit ensembles.

We show that in the limit of just few spins there appear Rabi-like oscillations between spin excitations and photon mode irrespective of the initial condition (excitation either in  spin or photon subsystem), as expected. However, as spin number increases, dynamics becomes highly sensitive to initial conditions. The most counterintuitive behavior is found for the initial condition of excitation in the spin subsystem – spin excited state becomes frozen through what we call Zeno-like effect <cit.>: its relaxation time grows with the number of spins, so that in the limit of infinite spin number the excited state does not decay at all. For finite systems, there appear periodic partial revivals of an excited state, while the period of revivals grows with the spin number. We also demonstrate that certain collective excitations encoded into the spin subsystem and characterized by a finite entanglement, such as an antisymmetric Bell state, become even more robust with respect to the influence of environment of remaining spins despite of the fact that entanglement is, in general, a very fragile entity. This collective state being constructed from a couple of spins, essentially does not decay at all even in presence of a bath of remaining spins, provided these two spins have excitation energies neighboring in the energy space. Nevertheless, for larger separation between spin excitation frequencies, the evolution remains highly sensitive to the symmetry of the wave function, i.e., to the minus or plus sign in the Bell state. This result highlights a nontrivial role of quantum interference effects for disordered ensembles in a mesoscopic regime. Since entanglement is a key resource for quantum technologies, while quantum interference effects are essential for the experimental demonstration of "quantumness" of artificial macroscopic systems, our conclusions might be important for applications. The Zeno-like effect we found is directly linked to the formation of a quasi-continuum of Hamiltonian eigenstates poorly coupled to light which we refer to as collective dark states. They have similarities with subradiant states of homogeneous model the latter states being totally uncoupled from light <cit.>. In the case of a single photon present in the system in the initial moment, dark eigenstates also significantly affect dynamics of a mesoscopic ensemble – they absorb the photon into the spin subsystem and then periodically release it giving rise to sharp revivals.

Our results are potentially useful for quantum states protection, storage and engineering. We believe that they can be also used to deeply probe quantum mechanical nature as well as collective properties of mesoscopic ensembles of artificial macroscopic spins, such as superconducting qubits, coupled to cavities.








§ HAMILTONIAN AND PRELIMINARIES


We consider an ensemble of two-level systems coupled to a single mode photon field. This coupled system is described by Dicke Hamiltonian of the form

    H=∑_j=1^L ϵ_j σ_j^+ σ_j^- + ω a^† a + g ∑_j=1^L (a^†σ_j^- + a σ_j^+),

where a^† and a correspond to the boson degree of freedom:

    [a, a^†]=1,

while σ_j^±, σ_j^z correspond to the paulion degrees of freedom and describe a set of L two-level systems:

    [σ_j^+, σ_j^-]=2σ_j^z, 
       
    [σ_j^z, σ_j^±] = ±σ_j^±.


The Hamiltonian (<ref>) commutes with the operator of the total pseudo-particle number, i.e., the number of bosons plus the number of excited two-level systems. Let us denote this number as M. Pseudo-particles of Dicke model are often referred to as excitations (of noninteracting system), but they should not be confused with excited states within a sector of given M (of interacting system). Namely, for any fixed M, there are different eigenstates of the Hamiltonian. At given M, the lowest energy state is the ground state, while others represent excited states. Note that ground state energies for different values of M can be also quite different. For example, if interaction constant g is large enough, a global ground state can be attained at some nonzero M. This behavior is closely related to the so called superradiant transition <cit.>.

Note that the Hamiltonian (<ref>) is based on rotating wave approximation, which neglects counterrotating terms of the form g(a σ_j^- + a^†σ_j^+). These terms do not conserve excitation number and it is known that they can be omitted provided the detuning between the cavity and spin is not too large, |ϵ_j - ω| ≪ω, see, e.g., Ref. <cit.>. However, counterotating processes have to be taken into account even in the resonance, but only in certain specific situations, such as parametric and periodic excitation of a coupled qubit-cavity system <cit.>. In the situation we consider in this article, counterrotating terms can indeed be safely neglected, since we are mostly interested in the interaction between the spin ensemble and cavity, which are close to the resonance, and do not treat such parametric excitations.

We also introduce an operator S^†(λ) defined as

    S^†(λ) = a^†+∑_j=1^Lg/λ-ϵ_jσ_j^+,

which is parametrized by the energy-like quantity (rapidity) λ. The state of the form

    ∏ _n=1^MS^†(λ_n) |↓↓…↓, 0 ⟩

is an eigenstate of the Hamiltonian, provided rapidities satisfy a set of Bethe equations
<cit.>

    λ_n-ω/g+∑_m≠ n^M 2g/λ_n - λ_m - ∑_j=1^L g/λ_n - ϵ_j=0,

while the eigenenergies E are expressed through the roots λ_n as

    E = ∑_n=1^M λ_n.

There are in general multiple solutions of Eq. (<ref>), i.e., many sets {λ_n }, which form an energy spectrum within a sector of a given M according to Eq. (<ref>).

We restrict ourselves to the regime of a single pseudo-particle, M=1. In this case, there is only single rapidity λ, which satisfies a single Bethe equation

    λ-ω/g - ∑_j=1^L g/λ - ϵ_j=0.

This is a polynomial equation of order L+1 and has the same number of solutions, which we refer to as λ^(α). It can be readily extracted from Eq. (<ref>) that all solutions are real. They correspond to different eigenstates of the Hamiltonian in the same sector M=1. The unnormalized eigenfunctions can be represented as

    |Φ_α⟩= S^†(λ^(α)) |↓↓…↓, 0 ⟩.

It is easy to find a norm as

    ⟨Φ_α|Φ_α⟩= 1+∑_j=1^Lg^2/(λ^(α)-ϵ_j)^2.

The normalized wave function thus reads as

    |φ_α⟩ = 1/√(⟨Φ_α|Φ_α⟩) |Φ_α⟩.


In Appendix A, we show how these results can be used to analyze system's evolution starting from different initial conditions, but corresponding to the same sector M=1. This number is of course conserved during the free evolution.

Note that, in absence of inhomogeneous broadening, Bethe states of Dicke model do not form a complete set (see, e.g., Ref. <cit.>), since degenerate nonradiating states, decoupled from light, cannot be obtained through Eq. (<ref>).



§ FLAT DISTRIBUTION WITH NEARLY CONSTANT DENSITY OF STATES




 §.§ Hamiltonian eigenstates


In the limit of a single spin L=1, the Hamiltonian (<ref>) reduces to the well known Jaynes-Cummings Hamiltonian. In this case, there exist only two solutions of Bethe equation (<ref>). These are shown schematically in Fig. 1, where a full resonance between spin excitation energy and photon frequency is assumed. In the case of many spins L, the set of Bethe roots becomes drastically different due to the splitting between spin energies. This situation is also illustrated schematically in Fig. 1 for distribution of spin energies having abrupt terminations. In this case, in addition to the two separated roots relevant for Jaynes-Cummings model, new roots do appear, which are confined between neighboring spin energies. Note that in the limit of strong interaction g, when splitting between spin energies is irrelevant, two separated roots become responsible for Rabi oscillations of the whole ensemble of spins with frequency g√(L).




Among physically meaningful distributions of the density of states induced by disorder are the Gaussian and Lorentzian  distributions <cit.> or q-Gaussian distribution relevant for NV-centers <cit.>. Furthermore, simplified equally-spaced distribution of ϵ between the two cut-offs is of fundamental importance, since it allows to grasp main features of system's dynamics produced by splitting of excitation frequencies. Physically, it might correspond to the broad distribution, for which a central part, most strongly interacting with the photon mode, has nearly constant density of states. In the next Section, we will briefly describe the effects due to distributions with smooth tails.

Now let us concentrate on such an equally-spaced distribution of spin energies ϵ_j with the difference between neighboring ϵ_j denoted as d. We assume that the width of the distribution Ω, i.e., the difference between maximum and minimum excitation energies, is independent on L, while the number of spins L is large. We then consider the limit L →∞ and find explicitly a leading order in L behavior as well as dominant corrections in 1/L essential for mesoscopic systems. In this limit, d → 0. In order to construct 1/L expansion, we assume that g√(L) is L-independent. Thus, d/g ∼ 1/√(L)≪ 1 in this limit, so that there are many spin excitation frequencies within the energy scale g. We also focus on the most interesting regime, when Ω and g√(L) are of the same order that results in a very rich phase diagram already for a static systems. While the former quantity is a natural scale to characterize broadening, the latter provides coupling energy between photon and spin subsystem in absence of broadening. Thus, at Ω∼ g√(L), there exists a pronounced competition between collective action of the whole ensemble of spins and their individual behavior  <cit.>. The situation we consider is of particular importance in the context of superconducting qubits-resonator systems, since state-of-art superconducting circuits seem to start entering such a regime, where collective properties become significant despite of the disorder <cit.>. Thus, we start from the limit of very few spins at g ≲Ω, which is addressed mainly numerically, and analyze the whole crossover to the limit g√(L)≫Ω with the particular focus on the intermediate mesoscopic regime g√(L)∼Ω.

In Appendix B, we address solutions to Bethe equation (<ref>) within our model. There are in total (L+1) solutions, (L-1) of them are confined between neighboring spin excitation frequencies, while two remaining roots are, in general, separated, as shown in Fig. 1. Physically, two additional roots and confined roots describe Hamiltonian eigenstates with quite different properties. Indeed, it follows from Eq. (<ref>) that each of these states |φ^(α)⟩ consists of a superposition of a single photon state and spin excited states with excitation energies detuned from λ^(α) by energy not too large, i.e., ≲ g. Therefore, eigenstates corresponding to confined roots are coupled essentially to each spin among a set of ∼ g/d ≫ 1 spins and to the photon mode. The coupling to the photon thus appears as quite weak. Therefore, such eigenstates can be characterized as dark states. Actually, they can be imagined as superpositions of many individual excited spins centered in energy space around a given spin, each superposition being only weakly coupled to light. In the limit of homogeneous model, these eigenstates should become fully decoupled from the light being gradually transformed to usual subradiant states. In contrast, two separated roots have smaller number of surrounding spin excitation energies, and coupling to the light for these two eigenstates is stronger, so that they can be refereed to as bright states. We would like to stress that collective dark states emerge only in the limit g/d ≫ 1, since each of them must be represented by a superposition of many individual spin states in order to be dark. Indeed in the regime of just few spins and at g ≲Ω coupling to the light for all eigenstates is significant.



 §.§ Dynamics of the system with single spin excited in the initial moment


Now we use our general results from Appendices A and B to study dynamics of the system with single spin excited in the initial moment, the excitation energy of this spin being ϵ_A. We rewrite the time dependent wave function (<ref>) in leading order as

    |ψ(t) ⟩≃d^2/g^2π^2∑ _α=1^L-1 e^-i λ^(α) tg/λ^(α)-ϵ_A1/1+1/π^2(lnϵ_L-ϵ_α/ϵ_α-ϵ_1)^2(a^†+∑_j=1^Lg/λ^(α)-ϵ_jσ_j^+)|↓↓…↓, 0 ⟩.

Let us stress that we omitted in Eq. (<ref>) a contribution from two separated roots (bright states), which is justified in this order. For the amplitude of the probability to find the initially excited spin still in this state we have

    ⟨ψ(t=0) |ψ(t) ⟩≃d^2/π^2∑ _α=1^L-1 e^-i λ^(α) t1/(λ^(α)-ϵ_A)^21/1+1/π^2(lnϵ_L-ϵ_α/ϵ_α-ϵ_1)^2.


Let us now focus on the situation, when all energies ϵ_j are centered around ω, while ϵ_A is also in a resonance with ω. It is easy to see that, under such conditions, δ_α≃ d/2 in a vicinity of ω. It is also clear that the logarithm in the right-hand side of Eq. (<ref>) can be omitted (at least for t not too large). We then obtain in leading order for t ≲ 2π/d

    ⟨ψ(t=0) |ψ(t) ⟩≃d^2/π^2 e^-iϵ_At∑ _α=1^L-1  e^-i (ϵ_α-ϵ_A) t1/(ϵ_α-ϵ_A+d/2)^2≃8/π^2 e^-iϵ_At∑_n=0^∞cos[(2n+1)td/2]/(2n+1)^2.

This function represents a simple periodic triangle wave of a period 4π / d. The initially excited spin decays on a time scale ∼ 1/d ∼ L. There is a revival after such a decay due to the finite size of the environment of other spins, i.e., the finiteness of the number of spins having frequencies close to ω. In other words, an initial excitation becomes redistributed between ∼ g/d spins most strongly interacting with an initially excited spin via the photon field. After some time, there occurs a refocusing of such a collective state into an initial state. At long times, an ideal periodic function (<ref>) becomes somehow smeared and less regular, since various subdominant contributions in 1/L come into play, which, in particular, results in certain finite occupations of strongly detuned spins.






The major conclusion deduced from Eq. (<ref>) is that the excited state decays on a time scale, which grows as the number of spins (density of states) in the system increases. Thus, the excited state becomes stabilized by the continuum of dark states giving rise to what we refer to as  Zeno-like effect. A similar effect of "radiation trapping" is known from literature for homogeneous Dicke model <cit.>. Thus, we show that such an effect also takes place for inhomogeneously broadened ensembles and we reveal how it emerges as the number of spins in such an ensemble grows.
Physically, this phenomenon might be attributed to the fact that the energy of the initial state with one of the spins excited is closer to the eigenenergies of dark states. The dynamics is then governed mostly by these dark states, which are weakly coupled to the light. Such an understanding has to be contrasted with naive expectations that an addition of extra 'parasitic' spins should only lead to a kind of a chaotization of Rabi oscillations between the given excited spin and the photon field. A chaotization occurs only in the limit of just few spins in the sense that in this limit dynamics is highly sensitive to precise detunings between spin excitation frequencies and the photon frequency. However, at L ≫ 1 it transforms to the universal triangle wave dependence sensitive only to the mean density of states in the vicinity of ω.

Let us now take into account a first order correction in 1/L. The most important contribution in this order stems from the fact that we have neglected two additional roots of Eq. (<ref>), i.e., the bright states. We can readily find that these states produce a correction to the amplitude of the probability (<ref>) given by

    δ⟨ψ(t=0) |ψ(t) ⟩≃1/L e^-iϵ_At∑_α=0,L  e^-i t (λ^(α)-ϵ_A)1/⟨Φ_0,L|Φ_0,L⟩g^2L/(λ^(α)-ϵ_A)^2.

Bright states give rise to Rabi-like oscillations contributing to the total dynamics.

Figure <ref> visualizes how dynamics changes as L grows. It shows the evolution of |⟨ψ(t=0) |ψ(t) ⟩| for L=4 (a), L=6 (b), L=10 (c), L=20 (d) and at g=0.05ω, Ω = 0.1ω, ϵ_A=ω. Two different results are compared – the first one is obtained numerically by solving Bethe equation and the second one is an explicit result, which contains a dominant contribution (<ref>) as well as a subdominant correction (<ref>). We indeed find quite good agreement between the numerical and explicit results starting from L≈10, while qualitative agreement exists even for smaller values of L. These plots visualize how Rabi oscillations at L=1 transform into triangle wave of very large period at L ≫ 1.

Despite of the universality in the large L limit, actual time evolution for the spin occupation can be very sensitive to mesoscopic fluctuations. In order to illustrate this, we plot in Fig. <ref> the same quantity |⟨ψ(t=0) |ψ(t) ⟩| calculated numerically at L=20, g=0.05ω, and Ω = 0.1ω, but for randomly distributed spin energies confined between the same cutoffs, while ϵ_A is taken to be closest to the resonance among all spins. Fig. <ref> (a) and (b) correspond to two typical realizations of disorder. In the first case, the evolution of |⟨ψ(t=0) |ψ(t) ⟩| is quite close to the similar dependence for the idealized equally-spaced distribution, while in the second case spin occupation oscillates with nearly the same period, but it does not reach zero. Such a behavior in the latter case is a direct consequence of mesoscopic fluctuations, i.e., finiteness of spin number L.





We also verify that Eq. (<ref>) is consistent with the normalization condition |⟨ψ(t=0) |ψ(t=0) ⟩|=1. Indeed, it is known that ∑_n=0^∞1/(2n+1)^2=π^2/8 this identity being connected to the famous Basel problem. We would like to stress that the way, this relation appeared in our formalism, seems to be highly nontrivial.


In contrast, occupations in the photon subsystem appear to be generally very small under the initial condition of single spin excited. Namely, the probability amplitude of finding a system in the state with a single photon vanishes in the limit L →∞. This result can be derived from Eq. (<ref>) yielding

    ⟨ 0, ↓…↓ | a |ψ(t) ⟩≃4d/g π^2e^-iϵ_At∑_n=0^∞sin[(2n+1)td/2]/2n+1,

which is a periodical square wave of the amplitude proportional to d/g. It is small provided there are many spin excitation energies within g. Note that two separated roots produce an additional contribution of the order of 1/L. This result is consistent with the fact that the coupling between the dark states and the photon state is very weak, while dark states are predominantly responsible for the system's dynamics in large L limit. Thus, "radiation trapping" effect for homogeneous model <cit.> is naturally recovered. If energy dissipation is present in the system and it is much larger for the cavity compared to spin subsystem, than a coupling to dark states allows to drastically reduce total effective dissipation.

The obtained results are potentially important in the context of quantum information storage and quantum state protection, since they show that dark eigenstates are able to greatly enhance the life time of the single excited spin and also they reveal how stabilization of excitations within spin subsystem induced by these states emerges as the number of spins increases. Perhaps, our findings can be also used to probe properties of mesoscopic ensembles of artificial spins coupled to a cavity. By exciting one of the spins of the ensemble via an additional waveguide and tracing its evolution, it is possible to see whether the ensemble behaves collectively according to the scenario we predict or such a behavior is suppressed due to the disorder and/or decoherence processes.




The problem we study might be also considered as an exactly solvable toy model for the coupled qubit-cavity system in presence of mesoscopic environment of 'parasitic' quantum two-level systems (fluctuators) interacting with the main system via photon degree of freedom and leading to decoherence. This model, however, is unable to account for the fluctuators own environment, so that it is applicable provided fluctuators are more strongly coupled to the qubit-cavity system than to their environment. The model has to be contrasted with the spin-boson model and it corresponds to pure quantum regime when an entanglement between the qubit-cavity system and fluctuators is of importance <cit.>. The existence of the exact solution allows to perform microscopic analysis for such a situation without switching to any simplified and not fully controllable approximation. In this context, we investigated the dynamics of a single qubit-cavity coupled system in presence of a mesoscopic ensemble of fluctuators with strongly detuned excitation frequencies. We found that such an ensemble produces very significant phase drift ('dephasing') for the amplitude of the probability of qubit excitation. We would like to stress that this effect is much stronger than the influence of fluctuators on qubit excited state population – even if Rabi oscillations are almost perfectly reproduced, the phase drift can be significant. Of course, it becomes larger as the detuning between the mean fluctuator frequency and photon mode is decreased. We also expect that presence of individual environments of qubits should lead to the suppression of Zeno-like effect, so that in the regime of strong energy dissipation lifetime of the qubit excited state is limited by the decay into its own environment.




 §.§ Dynamics of the system with single photon state in the initial moment


Let us now apply our general results for Bethe roots in the equally-spaced model to the dynamics of the system starting from a single photon state. For the amplitude of probability to have still a single photon state after some evolution, we obtain from (<ref>)

    ⟨ 0, ↓…↓ | a |ψ(t) ⟩≃d^2/π^2 g^2 e^-i d t/2∑ _α=1^L-1 e^-i ϵ_α t1/1+1/π^2(lnϵ_L-ϵ_α/ϵ_α-ϵ_1)^2+∑ _α=0,L e^-i λ^(α) t1/1+g^2L/(λ^(α)-ϵ_1)(λ^(α)-ϵ_L),


where we separated contributions from dark and bright states. In the regime of weak disorder, Ω≪ g√(L), a maximum of the absolute value of the first term as a function of t scales as ∼Ω^2/g^2L, while the second term starts to represent collective Rabi oscillations of the spin ensemble as a whole with frequency g√(L) and the amplitude approaching 1. In the intermediate regime g√(L)∼Ω, amplitudes of oscillations due to these two contributions are generally of the same order.

As an example, Fig. <ref> shows the evolution of the occupation of a single photon state for L=4 (a), L=12 (b), L=20 (c) and at g=0.05ω, Ω = 0.1ω calculated numerically. This figure evidences that the dynamics for small spin number is essentially chaotic. However, it becomes more ordered as number of spins grows and the system enters an intermediate regime Ω∼ g√(L). Initially, photon is mainly absorbed by a set of dark states, which transform the excitation into the spin subsystem. The photon state occupation in this limit is represented by a superposition of small-amplitude Rabi-like oscillations due to two bright states and quasi-periodical sharp revivals of a period ∼ 1/d due to the set of dark states, which release the excitation from the spin subsystem and then absorb it again. These revivals are only partial. The time delay for the first revival increases with the increase of the number of spins in the ensemble due to the growth of the dark state number. The amplitude of oscillations due to bright states grows as L increases, so that in the regime g√(L)≫Ω they have to be transformed to Rabi oscillations of the whole ensemble with frequency g√(L), while the role of dark states in the dynamics of photon subsystem and for the initial condition considered becomes negligible. Let us stress that, in contrary, they play a dominant role in the dynamics of spin subsystem in this limit for the initial condition considered in the preceding Subsection.

As for the dynamics of spin subsystem, we obtain for the amplitude of probability to have a single spin with energy ϵ_m excited

    ⟨ 0, ↓…↓ | σ_m |ψ(t) ⟩ = ∑ _α=0^L e^-i λ^(α) t1/1+∑_j=1^Lg^2/(λ^(α)-ϵ_j)^2g/λ^(α)-ϵ_m
    ≃4 d/π^2 g e^-i ϵ_m t∑_n=0^∞sin[(2n+1)td/2]/2n+1+∑ _α=0,L e^-i λ^(α) t1/1+g^2L/(λ^(α)-ϵ_1)(λ^(α)-ϵ_L)g/λ^(α)-ϵ_m,

where we again separated contributions from dark and two bright states. This probability is generally small.





 §.§ Dynamics of the system with Bell entangled states encoded in the initial moment into a continuum of spins


Let us now consider the dynamics of the system with constant density of states, when the initial state is one of the Bell states |χ_±⟩=1/√(2)(σ_A^+ ±σ_B^+)|↓↓…↓, 0 ⟩. These states thus can be symmetric and antisymmetric. They represent more sophisticated states which are characterized by a finite quantum entanglement and include more than one spin. The evolution, in this case, is described by Eq. (<ref>). It is known that, in the case of homogeneous model, the antisymmetric Bell state does not decay at all for any L <cit.>. It is of interest to explore an influence of broadening on its stability in the case of a mesoscopic ensemble.

An important quantity in this context is an overlap between the initial state of the system χ_± and its state |ψ(t) ⟩ after the beginning of the evolution

    ⟨χ_± |ψ(t) ⟩=1/2∑ _α=0^L e^-i λ^(α) t1/1+∑_j=1^Lg^2/(λ^(α)-ϵ_j)^2[g/λ^(α)-ϵ_A±g/λ^(α)-ϵ_B]^2.

In the large-L limit, this expression reduces to

    ⟨χ_± |ψ(t) ⟩≃8/π^2e^-i (ϵ_A+ϵ_B) t/2(cos(ϵ_B-ϵ_A/2t) ∑_n=0^∞cos[(2n+1)td/2]/(2n+1)^2±d/(ϵ_A-ϵ_B)sin(ϵ_B-ϵ_A/2t) ∑_n=0^∞sin[(2n+1)td/2]/2n+1).

The correspondent evolution of fidelity defined as |⟨χ_± |ψ(t) ⟩| and calculated numerically at L=20 for three different values of |ϵ_B-ϵ_A|, ϵ_A being in a resonance with ω, is plotted in Fig. <ref>. The agreement between the numerics and explicit result (<ref>) (not plotted in Fig. <ref>) is again spectacular.

We see that instead of the periodic triangle wave fidelity follows much more complicated evolution. The life time of the state χ_+, i.e., the time needed for the fidelity to drop from maximum to the first zero, is of the same order as for the single excited spin, i.e., it is enhanced due to Zeno-like effect. The situation, however, is very different for χ_-. Namely, for the minimum possible separation between energies of two spins A and B, |ϵ_B-ϵ_A|=d, the life time becomes infinite, since the fidelity oscillates periodically but always remains high. Such a behavior in the case of χ_- state is predictable for two isolated spins, but we see that it survives even in presence of an environment of other spins with excitation frequencies close to the resonance, which are expected to strongly interact with a couple of two given spins. In reality, the interaction is strongly suppressed due to a specific structure of a two-spin wave function. The fidelity also does not vanish in the case |ϵ_B-ϵ_A|=2d, but becomes suppressed in its minimum due to the influence of the intermediate spin. Starting from |ϵ_B-ϵ_A|=3d, the minimum fidelity reaches zero, but the life time of χ_- state nevertheless remains longer than the lifetime for χ_+ state. Thus, χ_- turns out to be much more robust than χ_+ especially for small separations between ϵ_A and ϵ_B. However, the effect of difference between ϵ_A and ϵ_B is quite strong and it reflects significant deviations of behavior from the case of homogeneous model. Nevertheless, the evolution of fidelity remains highly sensitive to the plus or minus sign in the Bell state, i.e., quantum interference effects in mesoscopic regime appear as quite robust with respect to broadening.

These results demonstrate that it is possible to stabilize in the spin subsystem collective quantum states involving few spins and not only single spin states. Surprisingly, they can be even more stable than single spin states and are able to support finite entanglement, in principle, for an arbitrary time, which is rather unexpectable, because entanglement is usually considered as a rather fragile entity.



Note that it also follows from Eq. (<ref>) that the first-order correction produced by two separated roots is suppressed in the case of χ_- state due to the minus sign in the sum. This is consistent with the results shown in Fig. <ref>.

Next, let us turn to photon subsystem and evaluate the amplitude of probability of having a single photon state. We readily obtain

    ⟨ 0, ↓…↓ | a |ψ(t) ⟩ = 1/√(2)∑ _α=0^L e^-i λ^(α) t[g/λ^(α)-ϵ_A±g/λ^(α)-ϵ_B] 1/1+∑_j=1^Lg^2/(λ^(α)-ϵ_j)^2.

In leading order in L, it reduces to

    ⟨ 0, ↓…↓ | a |ψ(t) ⟩≃2 √(2)/π^2d/g(e^-i ϵ_A t± e^-i ϵ_B t)
    ∑_n=0^∞sin[(2n+1)td/2]/2n+1

Due to the presence of the term (e^-i ϵ_A t± e^-i ϵ_B t) the occupation of the single photon state in the case of antisymmetric initial condition and at ϵ_A ≈ϵ_B is dramatically reduced compared to the case of a single excited spin considered above, where it is also small. Thus, such a state turns out to be "superdark". This is why χ_- is more robust compared to χ_+. In the case of homogeneous model, χ_- becomes a true Hamiltonian eigenstate totally decoupled from light (subradiant mode <cit.>).

The results of this Subsection might be also used for testing quantum mechanical nature of artificial spin-cavity systems including a manifestation of quantum interference effects, which are crucial for an unambiguous demonstration of "quantumness" of such engineered macroscopic systems. The entangled state between two qubits can be created via additional tunable cavities coupled to these particular qubits. By performing a two-qubit tomography it is then possible to follow the entanglement dynamics and to deeply probe coherent properties of qubits-cavity coupled systems. Particularly, the evolution must be sensitive to the symmetry of Bell state (plus or minus sign).



§ EFFECTS OF STATISTICS


Let us briefly discuss the effects arising due to deviation from equally-spaced distribution towards Gaussian or similar types of statistics characterized by smooth tails. We found numerically that if the typical energy scale associated with the deviation from equally-spaced distribution (such as mean variance for Gaussian distribution) is much larger than g, no qualitative change of behavior for dynamics starting from excitations within spin subsystem is observed at short times ∼ 2π d_0, d_0 being an inverse maximum density of states. Namely, leading order contributions to the quantities found above are qualitatively very similar. This happens because most strong interaction between the photon and spin subsystems is due to spins with energies close to the resonance with ω, while the influence of remaining spins become larger at long times. Nevertheless, subdominant contributions do change since no well defined bright state survives under Gaussian-like distributions with tails; this, for instance, results in absence of Rabi-like oscillations, which are transformed into small amplitude chaotic dependencies similar to the ones visible in Fig. <ref>.

However, rather significant changes do occur provided the scale associated with the deviation from equally-spaced distribution is lowered to ≲ g. In this case, dynamics becomes more irregular. Nevertheless, the typical relaxation time of a single spin excited remains to be long. Moreover, our general conclusion about extreme robustness of antisymmetric Bell state χ_-, which do not decay at all for neighboring spins A and B, remains valid as well.




§ CONCLUSIONS


We studied dynamics of mesoscopic ensemble of qubits coupled to a single mode cavity and concentrated on effects of disorder in qubits excitation frequencies and the regime of a weak excitation. In particular, we analyzed how collective properties of such a coupled system are formed as the number of qubits in the ensemble grows. This is done using a Bethe ansatz solution of Dicke model, which provides a direct access to Hamiltonian exact eigenstates as well as simple and pictorial understanding based on them.

We found that dark states weakly coupled to light gradually emerge and start to play a very important role upon the crossover from few-qubit systems to large ensembles. They are similar to subradiant modes of homogeneous Dicke model, which are totally decoupled from the light. The role of dark states in the free evolution dynamics is more important for initial conditions corresponding to excitations within qubit subsystem.

Our main conclusions are as follows:

(a) Despite of inhomogeneous broadening, single qubit excitation becomes stabilized in the infinite qubit number limit through what can be referred to as Zeno-like effect. It is induced by the gradual formation of dark states.

(b) Surprisingly, antisymmetric entangled Bell state constructed from the individual states of two qubits can be even more stable than single qubit excitation, provided excitation frequencies of these two qubits are not far from each other in the energy space. However, as separation grows, inhomogeneous broadening is able to suppress this effect. Nevertheless, the evolution remains  highly sensitive to the symmetry of the total wave function, which highlights a nontrivial role of quantum interference effects in mesoscopic regime.

(c) The effect of finite number of spins is generally manifested through partial revivals of the initial state, the period of revivals being proportional to the number of qubits. This scenario explains how Rabi-like oscillations in the limit of just few qubits are transformed into collective and highly cooperative behavior in the limit of many qubits.

(d) The role of dark states is less essential for initial condition corresponding to the single photon. They however are still of importance in the mesoscopic regime, when they are able to absorb photon into the qubit subsystem and then to release it quasiperiodically in time. A characteristic period of such revivals grows with the number of qubits, while the amplitude lowers. They coexist with Rabi-like oscillations between the whole ensemble of qubits and photon field so that in the infinite number limit such oscillations are naturally recovered.

Our results provide additional insights to physics of qubits-cavity coupled systems and might be of importance for quantum states engineering, protection, and storage. In particular, we believe that our theoretical predictions can be used to deeply probe both coherent and collective properties of mesoscopic ensembles of artificial spins coupled to cavities. For instance, Zeno-like effect is definitely linked to the formation of collective properties of the ensemble despite of the disorder in excitation frequencies. On the other hand, a sensitivity of the evolution of the entangled Bell state to the symmetry of the wave function is based on quantum interference effects, while experimental demonstration of such effects is crucial for the unambiguous evidence of true "quantumness" of engineered macroscopic artificial systems.






 Useful comments by A. N. Rubtsov, A. A. Elistratov, and S. V. Remizov are acknowledged. This work is supported by Advanced Research Foundation (project no. 7/076/2016-2020). W. V. P. acknowledges a support from RFBR (project no. 15-02-02128). Yu. E. L. acknowledges a support from RFBR (project no. 17-02-01134).








§ TIME DEPENDENCE


Let us consider the dynamics of the system starting from some initial state |ψ(t=0)⟩ corresponding to M=1 this number being conserved during the evolution. We expand the time-dependent wave function over the orthonormal basis |φ_α⟩ as

    |ψ(t) ⟩=∑ _α=0^L C_α(t) |φ_α⟩,

where

    C_α(t=0)=⟨φ_α | ψ(t=0) ⟩.

These coefficients at arbitrary t>0 can be readily found from the Schrödinger equation

    C_α(t)=C_α(t=0) exp (-i λ^(α) t).


Below we consider several initial conditions and general distribution of spin energies. These general results are then used to analyze system's dynamics explicitly for an equally-spaced distribution of spin energies.



 §.§ Single spin excited


Let us obtain the time dependent wave function corresponding to the initial state defined as

    |ψ(t=0)⟩= σ_A^+|↓↓…↓, 0 ⟩,

which describes a single spin excited, while all the remaining spins are in their ground states and there is no photon in the system. Such states have a direct relation to states engineered in Ref. by a spectral hole burning technique. It is of interest to explore its dynamics as the number of spins in the system grows. Note that in the limit L=1 there appear Rabi oscillations between the single spin and photon field. The naive expectation is that an addition of extra spins would just lead to a sort of a chaotization of Rabi oscillations. We will show that this is not the case.

We readily obtain from Eq. (<ref>)

    C_α(t=0)=⟨φ_α | ψ(t=0) ⟩ = 1/√(⟨Φ_α|Φ_α⟩)g/λ^(α)-ϵ_A.

In the explicit form, time dependent wave function reads

    |ψ(t) ⟩=∑ _α=0^L e^-i λ^(α) tg/λ^(α)-ϵ_A1/1+∑_j=1^Lg^2/(λ^(α)-ϵ_j)^2(a^†+∑_j=1^Lg/λ^(α)-ϵ_jσ_j^+)|↓↓…↓, 0 ⟩.





 §.§ Single photon


We may also consider the initial condition of another type, which corresponds to the single photon in a cavity

    |ψ(t=0)⟩ = a^†|↓↓…↓, 0 ⟩.

From Eq. (<ref>), we again find

    C_α(t=0)=⟨φ_α | ψ(t=0) ⟩ = 1.

In the explicit form, time dependent wave function is

    |ψ(t) ⟩=∑ _α=0^L e^-i λ^(α) t1/1+∑_j=1^Lg^2/(λ^(α)-ϵ_j)^2(a^†+∑_j=1^Lg/λ^(α)-ϵ_jσ_j^+)|↓↓…↓, 0 ⟩.





 §.§ Bell state encoded into the spin subsystem


Let us now consider a dynamics of the system with the initial state being one of the two Bell states encoded into two spins within the spin subsystem

    |χ_±⟩=1/√(2)(σ_A^+ ±σ_B^+)|↓↓…↓, 0 ⟩.

Using the developed approach, we arrive at the explicit form of the time dependent wave function

    |ψ(t) ⟩=1/√(2)∑ _α=0^L e^-i λ^(α) t[g/λ^(α)-ϵ_A±g/λ^(α)-ϵ_B] 1/1+∑_j=1^Lg^2/(λ^(α)-ϵ_j)^2(a^†+∑_j=1^Lg/λ^(α)-ϵ_jσ_j^+)|↓↓…↓, 0 ⟩.





§ ROOTS FOR EQUALLY-SPACED DISTRIBUTION


For the confined roots, it is convenient to represent λ^(α) as ϵ_α+δ_α, where δ_α < d. In order to find δ_α, we split the sum in Eq. (<ref>) into two contributions

    λ^(α)-ω/g^2 = ∑_k=0^min(α-2,L-α)[ 1/dk+d+δ_α-1/dk-δ_α] + 
    [∑_k=min(α-2,L-α)+1^α-21/dk+d+δ_α
    - ∑_k=min(α-2,L-α)+1^L-α1/dk-δ_α].

The first contribution is over the energies ϵ_j extending symmetrically from ϵ_α in both sides, while the second one includes contributions of the remaining energies ϵ_j. Since we consider the regime of large L ≫ 1, we may replace upper limits of summation in the first contribution by +∞. This term then reduces to

    ∑_k=0^+∞[ 1/dk+d+δ_α-1/dk-δ_α] = π/d(πδ_α/d).

We may also replace summation by integration in the remaining terms in Eq. (<ref>) as well as λ^(α) by ϵ_α in its left-hand side. We then arrive at the expression for δ_α given by

    δ_α≃d/π^-1( 1/π[d/g^2(ϵ_α - ω) + lnϵ_L-ϵ_α/ϵ_α-ϵ_1] ).

An important special case is a distribution centered around ω, which leads to δ_α≃ d/2 for ϵ_α approaching ω.

Using a similar method, we evaluate ⟨Φ_α|Φ_α⟩ as

    ⟨Φ_α|Φ_α⟩≃ g^2 π^2/d^21/sin^2 (πδ^(α)/d).

The expression of sin^2 (πδ^(α)/d) can be readily obtained from Eq. (<ref>) yielding

    ⟨Φ_α|Φ_α⟩≃ g^2 π^2/d^2[1+1/π^2(lnϵ_L-ϵ_α/ϵ_α-ϵ_1)^2].


Let us now consider two additional roots of Eq. (<ref>), which we denote as λ^(0) and λ^(L), while λ^(0) < ϵ_1 and λ^(L) > ϵ_L. In order to find these roots, we may replace summation by integration in this equation, which is allowed provided ϵ_1-λ^(0)≫ d and λ^(L)-ϵ_L ≫ d. The equations for λ^(0) and λ^(L) are identical and they read as

    λ^(0,L)-ω≃g^2L/Ωlnλ^(0,L)-ϵ_1/λ^(0,L)-ϵ_L.

These are transcendental equations, which can be solved numerically. Let us stress that the dependence of λ^(0,L) on g is nonanalytic. The norms for two corresponding eigenfunctions are

    ⟨Φ_0,L|Φ_0,L⟩≃ 1+ g^2 L/(λ^(0,L)-ϵ_1)(λ^(0,L)-ϵ_L).








1A. Reiserer and G.  Rempe, Rev. Mod. Phys. 87, 1379 (2015).

2B. Hacker,	S. Welte,	G. Rempe, and S. Ritter, Nature 536, 193 (2016).

3H. Bernien,	B. Hensen,	W. Pfaff,	G. Koolstra,	M. S. Blok,	L. Robledo,	T. H. Taminiau,	 M. Markham,	D. J. Twitchen,	L. Childress,	 and R. Hanson, Nature 497, 86 (2013).

MoiseevS. A. Moiseev, S. N. Andrianov, and F. F. Gubaidullin, Phys. Rev. A 82, 022311 (2010).

mem1J. Nunn, K. Reim, K. C. Lee, V. O. Lorenz, B. J. Sussman, I. A. Walmsley, and D. Jaksch, Phys. Rev. Lett. 101, 260502 (2008).

MolmerJ. H. Wesenberg, Z. Kurucz, and K. Mølmer, Phys. Rev. A 83, 023826 (2011).

Molmer1B. Julsgaard and K. Mølmer, Phys. Rev. A 88, 062324 (2013).

ShulgaK. V. Shulga, P. Yang, G. P. Fedorov, M. V. Fistul, M. Weides, and A. V. Ustinov, JETP Lett. 105, 47 (2017).

NorihybZe-Liang Xiang, Sahel Ashhab, J. Q. You, and Franco Nori, Rev. Mod. Phys. 85, 623 (2013).

GirvinR. J. Schoelkopf and S. Girvin, Nature 451, 664 (2008).

MachaP. Macha, G. Oelsner, J.-M. Reiner, M. Marthaler, S. André, G. Schön, U. Hübner, H.-G. Meyer, E. Ilichev, and A. V. Ustinov, Nat. Commun. 5, 5146 (2014).

ZagoskinA. M. Zagoskin, D. Felbacq, and E. Rousseau, EPJ Quantum Tech. 3, 2 (2016).

qubitYu. Makhlin, G. Schön, and A. Shnirman, Rev. Mod. Phys. 73, 357 (2001).

LukinP. L. Stanwix, L. M. Pham, J. R. Maze, D. Le Sage, T. K. Yeung, P. Cappellaro, P. R. Hemmer, A. Yacoby, M. D. Lukin, and R. L. Walsworth, Phys. Rev. B 82, 201201(R) (2010).

Krimer1S. Putz, A. Angerer, D. O. Krimer, R. Glattauer, W. J. Munro, S. Rotter, J. Schmiedmayer, and J. Majer, Nat. Phot. 11, 36 (2017).

Krimer2D. O. Krimer, B. Hartl, and S. Rotter, Phys. Rev. Lett. 115, 033601 (2015).

Krimer3D. O. Krimer, M. Zens, S. Putz, and S. Rotter, Laser  Photon. Rev. 10, 1023 (2016)

eng1 G. Bensky, D. Petrosyan, J. Majer, J. Schmiedmayer,
and G. Kurizki, Phys. Rev. A 86, 012310 (2012).

eng2 J. Cai, F. Jelezko, N. Katz, A. Retzker, and M.B Plenio,
New J. Phys. 14, 093030 (2012).

GaudinM. Gaudin, J. Phys. (Paris) 37, 1087 (1976).

LiebK. Hepp and E. H. Lieb, Ann. Phys. (NY) 76, 453 (1983).

WiegmannA. M. Tsvelick and P. B. Wiegmann, Adv. Physics 32, 453 (1983).

AndreevA. V. Andreev, V. Gurarie, and L. Radzihovsky, Phys. Rev. Lett. 93, 130402 (2004).

YuzbE. A. Yuzbashyan, V. B. Kuznetsov, and B. L. Altshuler, Phys. Rev. B 72, 144524 (2005).

Loss1O. Tsyplyatyev and D. Loss, Phys. Rev. A 80, 023803 (2009); ibid. 82, 024305 (2012).

StronglyC. Sträter, O. Tsyplyatyev, and A. Faribault, Phys. Rev. B 86, 195101 (2012).

ZenoB. Misra and E. C. G. Sudarshan, J. Math. Phys. 18, 756 (1977).

CummingsF. W. Cummings and A. Dorri, Phys. Rev. A 28, 2282 (1983).

GarrawayB. M. Garraway, Phil. Trans. R. Soc. A 369, 1137 (2011).

LozovikN. B. Narozhny, A. M. Fedotov, and Yu. E. Lozovik, Phys. Rev. A 64, 053807 (2001).

ShapiroD. S. Shapiro, A. A. Zhukov, W. V. Pogosov, and Yu. E. Lozovik, Phys. Rev. A 91, 063814 (2015); A. A. Zhukov, D. S. Shapiro, W. V. Pogosov, and Yu. E. Lozovik, Phys. Rev. A 93, 063845 (2016).

TalalaevO. Babelon and D. Talalaev, J. Stat. Mech. 0706, 06013 (2007).

DukJ. Dukelsky, S. Pittel, and G. Sierra, Rev. Mod. Phys. 76, 643 (2004).

YudsonV. I. Rupasov and V. I. Yudson, JETP 60, 927 (1984).

LittlewoodP. R. Eastham and P. B. Littlewood, Phys. Rev. B 64, 235101 (2001).

NucPhys2017W. V. Pogosov, D. S. Shapiro, L. V. Bork, and A. I. Onishchenko, Nucl. Phys. B 919, 218 (2017).

GalperinH. J. Wold, H. Brox, Y. M. Galperin, and J. Bergli, Phys. Rev. B 86, 205404 (2012).






