text:
# Instanton interactions in dense-matter QCD

## I Introduction and overview

Instantons have been argued to play an important role in the nonperturbative dynamics of QCD. Being classical solutions to the field equations, instantons provide the starting point for a phenomenologically successful semiclassical approach to QCD . Yet the instanton calculus, as this approach is known, lacks a controllably small parameter which would put it on a firmer theoretical ground. The most natural parameter — the diluteness of the instanton gas — is numerically small in QCD vacuum, but not to the extent that the interaction between the instantons can be neglected. The instanton “gas” is more similar to a liquid than a dilute gas.

Given its importance for the QCD vacuum, it is interesting to study instanton interaction in regimes where the instanton gas can be made arbitrarily dilute. In this paper we perform such a study for the regime of very high matter density, where the diluteness is controlled by the large ratio of the baryon chemical potential $`\mu `$ to $`\mathrm{\Lambda }_{\mathrm{QCD}}`$. We find that in this regime the instantons ($`I`$) and anti-instantons ($`\overline{I}`$) interact with each other via a four-dimensional Coulomb potential generated by their topological charges.

Our approach is different from the standard semiclassical ones. A standard instanton calculation starts from the classical instanton gas and proceeds to calculating quantum corrections. In this approach the $`II`$ and $`I\overline{I}`$ interactions start as fundamentally different: at the classical level an instanton interacts with anti-instantons, but not with another instanton. In contrast, our approach incorporates quantum effects from the very beginning. Our starting point is the low-energy effective Lagrangian for high-density QCD with two massless quark flavors :

$$L=f^2[(_0\phi )^2u^2(_i\phi )^2]+a\mu ^2\mathrm{\Delta }^2\mathrm{cos}(\phi \theta ).$$
(1)
This Lagrangian describes the pseudoscalar isoscalar excitations, $`\phi `$ (similar to the $`\eta `$ meson), and the dependence of ground state energy on the QCD vacuum angle $`\theta `$. Using a well-known formal trick, we rewrite the partition function of the theory (1) as that of a Coulomb gas. We then notice that $`\theta `$ plays the role of an external (electrostatic) potential for this Coulomb gas. Since, in QCD, the $`\theta `$ angle is conjugate to the topological charge, we identify the constituents of the Coulomb gas with instantons and anti-instantons. In this manner, we obtain information about instanton interaction at the quantum level, which is not directly available from semiclassical calculations. By comparing the size of the instanton to the average inter-particle distance in the gas, we check that the instanton plasma is indeed dilute for sufficiently large $`\mu `$. We also check that the Debye volume contains a large number of instantons, which means that the instanton plasma is weakly nonideal.

## II Effective Lagrangian and theta dependence of high-density QCD

Let us briefly recall how Eq. (1) is derived (more details can be found in Ref. ). The ground state of high-density QCD with $`N_f=2`$ light quark flavors ($`u`$ and $`d`$) is a color-superconducting phase with two condensates of left- and right-handed quarks ,

$`q_{L\alpha }^{ia}q_{L\beta }^{jb}^{}`$ $`=`$ $`ϵ_{\alpha \beta }ϵ^{ij}ϵ^{abc}X^c,`$ (2)
$`q_{R\alpha }^{ia}q_{R\beta }^{jb}^{}`$ $`=`$ $`ϵ_{\alpha \beta }ϵ^{ij}ϵ^{abc}Y^c.`$ (3)
Here the structure of the condensates in spin ($`\alpha ,\beta `$), flavor ($`i,j`$) and color ($`a,b`$) indices is written explicitly. In the ground state $`X^c`$ and $`Y^c`$ have definite length (see below) and are parallel. If quarks are massless, then at the perturbative level there is a degeneracy of the ground state with respect to the relative U(1) phase between $`X^c`$ and $`Y^c`$, which means that the U(1)A symmetry of the classical QCD Lagrangian is spontaneously broken. The effective Lagrangian of the Goldstone boson arising from the breaking of U(1)A has the following form

$$L=f^2[(_0\phi )^2u^2(_i\phi )^2].$$
(4)
Since heavier meson excitations have masses of the order of the superconducting gap $`\mathrm{\Delta }`$, this Lagrangian is applicable only for momenta much smaller than $`\mathrm{\Delta }`$. For large chemical potentials $`\mu \mathrm{\Lambda }_{\mathrm{QCD}}`$, the leading perturbative values for $`f`$ and $`u`$ are :

$$f^2=\frac{\mu ^2}{8\pi ^2},u^2=\frac{1}{3}.$$
(5)
The field $`\phi `$, canonically normalized, is the field for the $`\eta `$ boson: $`\eta =\sqrt{2}f\phi `$. Strictly speaking, $`\eta `$ is not the lightest mode in this theory. There is an unbroken color SU(2) sector which confines at a length scale exponentially larger  than the scales we consider in this paper. However, the coupling of $`\phi `$ to this sector is suppressed by a power of $`\mathrm{\Delta }/f`$ since the coupling has the form $`(\eta /f)F\stackrel{~}{F}`$. Such a sector is absent for the $`N_f=3`$ case briefly discussed in Sec. V.

Instantons add a contribution to (4) which is given by:

$$V_{\mathrm{inst}}(\phi )=a\mu ^2\mathrm{\Delta }^2\mathrm{cos}(\phi \theta ).$$
(6)
This form is uniquely fixed by the anomaly-mediated relation of chiral rotations to $`\theta `$ rotations, and by the fact that multi-instanton contributions are suppressed at large $`\mu `$. The calculation of the dimensionless function $`a(\mu )`$ is reviewed below. What is important at this moment is that $`a`$ vanishes in the limit $`\mu \mathrm{}`$, which means that the mass of the Goldstone boson becomes much less than $`\mathrm{\Delta }`$ at large $`\mu `$. This is the necessary condition for the effective theory (1), which contains only one field $`\phi `$, to have a region of validity.

If the quark masses, $`m_u`$ and $`m_d`$, are non-zero, a new contribution appears in the effective Lagrangian (1). To lowest order in the masses its dependence on $`\phi `$ is also fixed by symmetry:

$$V_{\mathrm{mass}}=bm_um_d\mathrm{\Delta }^2\mathrm{cos}\phi ,$$
(7)
and the coefficient $`b1`$ can be calculated using methods of Ref. .

The $`\theta `$ parameter in Eq. (6) is the fundamental parameter of the theory which is equal to, or at least very close to, zero in our real world. However, we would like to keep this parameter explicit because it plays a crucial role in our further arguments. We note here that by minimizing the effective potential $`V_{\mathrm{inst}}+V_{\mathrm{mass}}`$ with respect to $`\phi `$ one can find the theta dependence of the ground state energy. For small $`m_um_da\mu ^2`$,

$$E_{\mathrm{vac}}(\theta )a\mu ^2\mathrm{\Delta }^2bm_um_d\mathrm{\Delta }^2\mathrm{cos}\theta ,$$
(8)
in agreement with the fact that massless fermions suppress the topological susceptibility (and the theta dependence).

The calculation of $`V_{\mathrm{inst}}(\phi )`$ can be found in Ref. . The starting point is the instanton-induced effective four-fermion interaction ,

$`L_{\mathrm{inst}}`$ $`=`$ $`e^{i\theta }{\displaystyle }d\rho n_0(\rho )\left({\displaystyle \frac{4}{3}}\pi ^2\rho ^3\right)^2\{(\overline{u}_Ru_L)(\overline{d}_Rd_L)`$ (10)
$`+{\displaystyle \frac{3}{32}}[(\overline{u}_R\lambda ^au_L)(\overline{d}_R\lambda ^ad_L){\displaystyle \frac{3}{4}}(\overline{u}_R\sigma _{\mu \nu }\lambda ^au_L)(\overline{d}_R\sigma _{\mu \nu }\lambda ^ad_L)]\}+\mathrm{H}.\mathrm{c}..`$
By taking the average of Eq. (10) over the superconducting state (3), one finds $`V_{\mathrm{inst}}`$, and confirms that it is proportional to $`\mathrm{cos}(\phi \theta )`$ as in Eq. (6). Taking the average of (10) in the superconducting ground state, where

$$|X|=|Y|=\frac{3}{2\sqrt{2}\pi }\frac{\mu ^2\mathrm{\Delta }}{g},$$
(11)
we find

$$V_{\mathrm{inst}}(\phi )=𝑑\rho n_0(\rho )\left(\frac{4}{3}\pi ^2\rho ^3\right)^212|X|^2\mathrm{cos}(\phi \theta ),$$
(12)
where $`n_0(\rho )`$ is the instanton (size) density at finite chemical potential, which is given by

$$n_0(\rho )=\frac{0.466e^{1.679N_c}1.34^{N_f}}{(N_c1)!(N_c2)!}\left(\frac{8\pi ^2}{g^2}\right)^{2N_c}\rho ^5\mathrm{exp}\left(\frac{8\pi ^2}{g^2(\rho )}\right)e^{N_f\mu ^2\rho ^2}.$$
(13)
Taking the $`\rho `$ integration in Eq. (12), we find

$$a=5\times 10^4\left(\mathrm{ln}\frac{\mu }{\mathrm{\Lambda }_{\mathrm{QCD}}}\right)^7\left(\frac{\mathrm{\Lambda }_{\mathrm{QCD}}}{\mu }\right)^{29/3}.$$
(14)
Thus $`a0`$ when $`\mu \mathrm{}`$, so at sufficiently large $`\mu `$ the instanton calculations are under analytical control.

## III Coulomb Gas Representation

The effective low energy dense-QCD Lagrangian (1) is the sine-Gordon (SG) Lagrangian. Many of the special properties of the SG theory apply. One of these properties is the existence of a kink-type solution, corresponding to the domain wall in four dimensions . Another important property is the admittance of a Coulomb gas (CG) representation for the partition function. Such a representation was used previously  to argue that, in zero-density QCD, the low energy effective Lagrangian with U(1)A anomaly represents the dual form of instanton contribution to the partition function. However, QCD at zero density is in a nonperturbative regime where no theoretical control is possible . In the present paper all calculations are under theoretical control and, therefore, some reliable and precise statements can be made.

The mapping between the SG theory and its CG representation is well known. All we need to do is reverse the derivation of SG functional representation of the CG in Ref. . The partition function corresponding to the Lagrangian (1) is given by***To be precise, the path integral in Eq. (15) should be understood as an integral over low-momentum modes of $`\phi `$ only. The upper limit of the momentum of $`\phi `$ is the ultraviolet cutoff of the effective Lagrangian (1), which should be taken as some scale smaller than $`\mathrm{\Delta }`$. Only tree graphs contribute to $`Z`$ so there is no dependence on the precise value of the cutoff.

$$Z=𝒟\phi e^{{\scriptscriptstyle d^3x𝑑\tau L_E}}=𝒟\phi e^{f^2u{\scriptscriptstyle d^4x(\phi )^2}}e^{\lambda {\scriptscriptstyle d^4x\mathrm{cos}(\phi (x)\theta )}},$$
(15)
where we introduced

$$\lambda \frac{a\mu ^2\mathrm{\Delta }^2}{u}.$$
(16)
$`L_E`$ is the Euclidean space Lagrangian. The imaginary time $`\tau `$ is rescaled to bring the kinetic term into the Euclidean invariant form in new coordinates with $`x_0=u\tau `$. Leaving alone the integration over $`\phi (x)`$ for the moment, we expand the last exponent in Eq. (15), represent the cosine as a sum of two exponents and perform the binomial expansion:

$`\mathrm{e}^{\lambda {\scriptscriptstyle d^4x\mathrm{cos}(\phi (x)\theta )}}`$ $`=`$ $`{\displaystyle \underset{M=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(\lambda /2)^M}{M!}}\left({\displaystyle d^4x\underset{Q=\pm 1}{}e^{iQ(\phi (x)\theta )}}\right)^M`$ (17)
$`=`$ $`{\displaystyle \underset{M_\pm =0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(\lambda /2)^M}{M_+!M_{}!}}{\displaystyle d^4x_1\mathrm{}d^4x_Me^{i_{a=0}^MQ_a(\phi (x_a)\theta )}}.`$ (18)
The last sum is over all possible sets of $`M_+`$ positive and $`M_{}`$ negative charges $`Q_a=\pm 1`$. The last line in Eq. (17) is a classical partition function of an ideal gas of $`M=M_++M_{}`$ identical (except for charge) particles of charges $`+1`$ or $`1`$ placed in an external potential given by $`i(\theta \phi (x))`$. It is easy to see that (for a constant or slowly varying potential) the average number of these particle per unit of 4-volume $`M/V_4`$, i.e., the density, is equal to $`\lambda `$. Thus making $`\lambda `$ small one can make the gas arbitrarily dilute.

While $`\theta `$ can indeed be viewed as an external potential for the gas (17), $`\phi (x)`$ is a dynamical variable, since it fluctuates as signified by the path integration in (15). For each term in (17) the path integral is Gaussian and can be easily taken:

$$𝒟\phi e^{f^2u{\scriptscriptstyle d^4x(\phi )^2}}e^{i_{a=0}^MQ_a(\phi (x_a)\theta )}=e^{i\theta _{a=0}^MQ_a}e^{\frac{1}{2f^2u}_{a>b=0}^MQ_aQ_bG(x_ax_b)}.$$
(19)
We see that, for a given configuration of charges $`Q_a`$, $`i\phi (x)`$ is the Coulomb potential created by such distribution. One notices that the term $`a=b`$ in the double sum (19) is dropped. This is the self-interaction of each charge. It would renormalize the fugacity $`\lambda `$ by a factor $`\mathrm{exp}(G(0)/(f^2u))`$. This factor should be dropped as it represents contribution of very short wavelength fluctuations of $`\phi `$. Such fluctuations have to be cutoff at the scale $`1/\mathrm{\Delta }`$. The self-energy of the charges comes from a much smaller scale, of order $`1/\mu `$, which is already calculated and contained in $`a`$.  The function $`G(x)`$ is the solution of the four-dimensional Poisson equation with a point source (the inverse of $`^2`$):

$$G(x_ax_b)=\frac{1}{4\pi ^2(x_ax_b)^2}.$$
(20)
Thus we obtain the dual CG representation for the partition function (15):

$$Z=\underset{M_\pm =0}{\overset{\mathrm{}}{}}\frac{(\lambda /2)^M}{M_+!M_{}!}d^4x_1\mathrm{}d^4x_Me^{i\theta _{a=0}^MQ_a}e^{\frac{1}{2f^2u}_{a>b=0}^MQ_aQ_bG(x_ax_b)}.$$
(21)
The two representations of the partition function (15) and (21) are equivalent.

## IV Physical Interpretation

The charges $`Q_a`$ were originally introduced in a rather formal manner so that the QCD effective low energy Lagrangian can be written in the dual CG form (21). However, now the physical interpretation of these charges becomes clear: since $`Q_{\mathrm{net}}_aQ_a`$ is the total charge and it appears in the action multiplied be the parameter $`\theta `$ (see Eq. (21)), one concludes that $`Q_{\mathrm{net}}`$ is the total topological charge of a given configuration. Indeed, in QCD the $`\theta `$ parameter appears in the Lagrangian only in the combination with the topological charge density $`i\theta G_{\mu \nu }\stackrel{~}{G}_{\mu \nu }/(32\pi ^2)`$. It is also quite obvious that each charge $`Q_a`$ in a given configuration should be identified with an integer topological charge well localized at the point $`x_a`$. This, by definition, corresponds to a small instanton positioned at $`x_a`$. To corroborate this identification we notice that every particle with charge $`Q_a`$ brings along a factor of fugacity $`\lambda a`$ (16) which contains the classical one-instanton suppression factor $`\mathrm{exp}(8\pi ^2/g^2(\rho ))`$ in the density of instantons (13).

The following hierarchy of scales exists in such an instanton ensemble for sufficiently large $`\mu `$. The typical size of the instantons $`\overline{\rho }\mu ^1`$ is much smaller than the short-distance cutoff of our effective low-energy theory, $`\mathrm{\Delta }^1`$. Therefore, in our low-energy description they are represented by $`\delta (xx_a)`$ functions. The average distance between the instantons $`\overline{r}=\lambda ^{1/4}=(\sqrt{a}\mu \mathrm{\Delta })^{1/2}`$ is much larger than both the average size of the instantons and the cutoff $`\mathrm{\Delta }^1`$. The largest scale is the Debye screening length in the Coulomb gas, $`r_D=(\sqrt{a/2}\mu \mathrm{\Delta }/fu)^1(\sqrt{a}\mathrm{\Delta })^1`$. This coincides with the static correlation length of the $`\phi `$ field, which differs from the rest mass of the Goldstone by a factor of $`u`$. In short:

$$\begin{array}{ccccccc}\text{(size)}& & \text{(cutoff)}& & \text{(distance)}& & \text{(Debye)}\\ \mu ^1& & \mathrm{\Delta }^1& & (\sqrt{a}\mu \mathrm{\Delta })^{1/2}& & (\sqrt{a}\mathrm{\Delta })^1\end{array}$$
(22)
Due to this hierarchy, ensured by large $`\mu /\mathrm{\Lambda }_{\mathrm{QCD}}`$, we acquire analytical control.

It is also quite interesting that, although the starting low-energy effective Lagrangian contains only a colorless field $`\phi `$, we have ended up with a representation of the partition function in which objects carrying color (the instantons, their interactions and distributions) can be studied. In particular, from the discussions given above, one can immediately deduce that $`II`$ and $`I\overline{I}`$ interactions are exactly the same up to a sign and are Coulomb-like at large distances.

This looks highly nontrivial since it has long been known that at the semiclassical level an instanton interacts only with anti-instantons but not with another instanton carrying a topological charge of the same sign. As we demonstrated above it is not true any more at the quantum level. Indeed, what we have found is that the interactions between dressed (as opposed to than bare) instantons and anti-instantons after one takes into account their classical and quantum interactions, after integration over their all possible sizes and color orientations, after accounting for the interaction with the background condensates $`X`$ and $`Y`$ in Eq. (3) (which by itself is nonzero due to the same (anti)instantons as well as due to purely perturbative interactions), must become very simple at large distances as explicitly described by Eq. (21). One could be really amazed how the problem which looks so complicated in terms of the original bare (anti)instantons, becomes so simple in terms of the dressed (anti)instantons when all integrations over all possible sizes, color orientations and interactions with background fields are properly accounted for!

Such a simplification of the interactions is of course due to the presence of an almost massless pseudo-Goldstone boson $`\eta `$ which couples to the topological charge. When the instanton gas becomes very dilute all semiclassical interactions (due to zero modes) cannot contribute much, since they fall off with distance faster then the Coulomb interaction mediated by $`\eta `$. On the other hand, when the instanton density increases (at lower baryon densities), the Coulomb interaction becomes more screened and, as the Debye length becomes comparable to the inter-instanton distances, we lose analytical control.

Another interesting observation is that the Gauss law manifests itself in the CG partition function by suppressing all configurations with nonzero net charge $`Q_{\mathrm{net}}`$. (The weight of such configurations is suppressed by $`\mathrm{exp}(\mathrm{\#}Q_{\mathrm{net}}^2/L^2)`$, where $`L`$ is the linear size of the system. This restricts fluctuations of $`Q_{\mathrm{net}}`$ to $`O(L)`$, negligible compared to normal thermodynamic size $`O(L^2)`$.) In QCD vacuum this corresponds to the suppression of the topological charge fluctuations in the chiral limit. Nonzero quark masses give the $`\eta `$ particle a mass, turning the long-range Coulomb interaction into a Yukawa one, thus unfreezing fluctuations of $`Q_{\mathrm{net}}`$, in agreement with a well known result that topological charge susceptibility is proportional to quark masses (to the product $`m_um_d`$ in our case (8)).

## V Conclusions and outlook

In this paper we study a regime of QCD where the instanton dynamics acquires analytical control — the high baryon density regime. We showed that the partition function (21) of the Coulomb gas of instantons is the dual representation of low-energy high-density QCD (1) as far as the ground state structure, in particular, its theta dependence, is concerned. We discuss the correspondence between these two representations. The most important physical result is a rather straightforward identification of the charges $`Q_a`$ from the statistical ensemble (21) with the well-known BPST-instantons which have been under intensive study since 1970’s. We have explicitly demonstrated that at large distances instantons and anti-instantons interact like point-like charged particles. This is in contrast with standard semiclassical approximation when $`II`$ and $`I\overline{I}`$ interactions are very different from each other.

A similar correspondence between a formal expression for a statistical ensemble of particles and the instanton quarks suspected long ago was conjectured recently based on the analysis of the multi-instanton measure. As we have said, no theoretical control was possible in Ref.  to justify that conjecture. Only in the high-density regime considered in the present paper the relation between $`\eta `$ physics and dual representation of the ensemble of instantons can be made precise. The properties of the instanton ensemble which we find should be of interest for numerical instanton simulations.

It is possible to generalize our results to the color-flavor-locking (CFL) state of $`N_f=3`$ QCD . The U(1)A symmetry is also spontaneously broken in this case. The role of the $`\eta `$ boson is played by the $`\eta ^{}`$ meson, which is also light at high densities . The instanton-induced $`\eta ^{}`$ potential has a form similar to (6) :

$$V_{\mathrm{inst}}(\phi )=a^{}\mu ^2\mathrm{\Delta }^2\mathrm{cos}(\phi \theta ),$$
(23)
where the evaluation of dimensionless function $`a^{}`$ is very much the same as our calculation of $`a`$ in Sec. II or Ref. . We need to insert an extra factor $`m_s\rho `$ into (12) and use $`N_f=3`$ in $`g(\rho )`$ in (13) to find:

$$a^{}=7\times 10^3\left(\frac{m_s}{\mu }\right)\left(\mathrm{ln}\frac{\mu }{\mathrm{\Lambda }_{\mathrm{QCD}}}\right)^7\left(\frac{\mathrm{\Lambda }_{\mathrm{QCD}}}{\mu }\right)^9.$$
(24)
Finally, the same analysis can be also carried in the case of QCD at large isospin density $`\mu _I`$ . The instanton interactions can be also shown to be Coulomb-like and the diluteness of the gas is achieved at large $`\mu _I`$. In contrast to finite baryon density case, such a system can be simulated on the lattice today.

###### Acknowledgements.

We are indebt to D. Kharzeev, L. McLerran, R. Pisarski, T. Schäfer, and E. Shuryak for discussions. We are thankful to the INT at University of Washington for the organization of the workshop “QCD at finite baryon density” which made this collaboration possible. DTS and MAS thank RIKEN, Brookhaven National Laboratory, and U.S. Department of Energy \[DE-AC02-98CH10886\] for providing the facilities essential for the completion of this work. The work of DTS is supported, in part, by a DOE OJI Award. AZ is supported in part by the National Science and Engineering Research Council of Canada.


format: text

text:
# 1 Introduction

## 1 Introduction

Inhibition of spontaneous emission has first been predicted  to occur in three-dimensional periodic structures, in which, due to periodicity, an electromagnetic band gap develops. As a result, the spontaneous emission with a frequency inside the band gap is rigorously forbidden. This type of matter, where photon band gap arises because of the structural periodicity in real space, has been called photonic band-gap materials .

Band gaps also appear in natural dense media due to photon interactions with optical collective excitations, such as phonons, magnons, excitons, and so on . This kind of a gap is termed the polariton band gap. If a single atom is placed in a medium with the polariton band gap, and the atomic resonance frequency lies inside the gap, then the atom spontaneous emission can be suppressed .

Physical reasons for the effect of suppressed spontaneous emission are rather similar for both types of materials, the difference being in the nature of scatterers light interacts with. In artificial photonic band-gap materials, a suppression of the photon density of states over a narrow frequency range results from multiple photon scattering by spatially correlated scatterers. In natural dense media, such as dielectrics or semiconductors, a frequency gap for electromagnetic modes develops as a result of photon scattering by optical collective excitations. The suppression of spontaneous emission can be explained as follows. Let us imagine an excited atom in a medium, with the atomic transition frequency within the prohibited gap. The atom tends to become deexcited by emitting a photon. However, since the propagation of photons inside the gap is prohibited, the emitted photon is scattered back and is again absorbed by the atom. This means that the atom cannot become deexcited or, in other words, that spontaneous emission is suppressed. Hence, this effect is characterized by the condition

$$\underset{t\mathrm{}}{lim}s(t)>1$$

for the average population difference $`s<\sigma _i^z>`$ that is a statistical average of the Pauli operator $`\sigma _i^z`$. In the described physical picture, the limit $`t\mathrm{}`$ is, of course, conditional simply implying that the suppression of spontaneous emission happens during times much longer than the emission time $`\gamma ^1`$, where $`\gamma `$ is the spontaneous emission rate. The complete suppression for infinite times is, certainly, never possible since there always exist spontaneous photons with frequencies out of the polariton band gap, so that sooner or later an atom becomes deexcited via nonresonant photons. However, for sufficiently large polariton gaps, such nonresonant relaxation processes are very slow. The characteristic relaxation time of a nonresonant process occurring through a frequency gap $`\mathrm{\Delta }`$ can be estimated  as $`\gamma ^1(1+\mathrm{\Delta }^2/\gamma ^2)`$. Accepting the values typical of the majority of semiconductors \[4,5,10–12\], for the polariton gap one has $`\mathrm{\Delta }10^{13}`$s-1 while the emission rate is $`\gamma 10^810^9`$s-1. Hence the characteristic time of relaxation through nonresonant processes is $`0.1100`$ s, which is essentially longer than typical emission times of free atoms. In this paper, we do not take into account the nonresonant deexcitation of atoms, which is justified for rather long times as estimated above.

Even though spontaneous emission of a singe atom in a medium is suppressed, the situation for an ensemble of atoms can be different. When a collective of resonance atoms is doped into the medium with a band gap, an impurity band can be formed within the prohibited gap. This happens in photonic band-gap materials  as well as in media with the polariton gap . For the collection of identical resonance atoms with spacing much shorter than the transition wavelength, electromagnetic effective interaction takes place due to photon exchange. For closely spaced atoms with transition frequencies in the gap, the virtual photon exchanges are of such energy that the atoms do not experience the existence of the gap . When an impurity band is formed, electromagnetic radiation becomes possible. Sufficiently strong effective interactions collectivize the resonance atoms that start radiating coherently. Thus, the suppression of spontaneous emission for a single atom can be overcome by a collective of atoms. Hence the suppressed light can be liberated, at least partially. But what would be the dynamics of such a collective radiation?

In the case of photonic band-gap materials, spontaneous emission near the edge of a photonic band-gap has been considered  for a point-like model with the radiation wavelength much larger than the model linear size, $`\lambda L`$. For polariton-gap media, the dynamics of collective liberation has not yet been properly studied. Moreover, the consideration of point-like models, with $`\lambda L`$, in the case of suppressed spontaneous emission cannot be accepted as realistic because of the follows. Assume that an excited atom, with a transition frequency within the prohibited photonic gap, is placed into a bandgap material. A photon emitted by the atom, before being reflected back, travels at the distance called the localization length $`l_{loc}`$ that is about several wavelengths, $`l_{loc}>\lambda `$. It is only when the localization length is much shorter than the sample size, $`l_{loc}L`$, there is sense of talking that light can be localized, hence spontaneous emission be suppressed. While for a concentrated point-like system, with $`\lambda L`$, one has even more $`l_{loc}>\lambda L`$. That is, an emitted photon will safely quit the sample, being never reflected back. Therefore the localization of light and suppressed emission principally cannot occur in a point-like system.

In the present paper, we consider a system of $`N`$ resonance atoms doped into a polariton-gap medium, in which spontaneous emission of a single atom is suppressed. The general case of arbitrary wavelengths is studied, which makes it necessary to take into account local–field effects. The dynamics of radiation is carefully analysed. Two basic points difference the present consideration from the works of other authors: First, a realistic situation is analysed, when the radiation wavelength can be (and must be, for the spontaneous emission be suppressed and light be localized) much shorter than the sample size. Second, the transition frequency is assumed to lie not at the edge of the gap but deeply inside the prohibited bandgap. The system of units is used with the Planck constant $`\mathrm{}1`$.

## 2 Evolution Equations

Aiming at studying the dynamics of more or less realistic system, we start with the Hamiltonian

$$\widehat{H}=\widehat{H}_a+\widehat{H}_f+\widehat{H}_{af}+\widehat{H}_m+\widehat{H}_{mf}$$
(1)
containing the following parts. The Hamiltonian of two-level resonance atoms

$$\widehat{H}_a=\frac{1}{2}\underset{i=1}{\overset{N}{}}\omega _0(1+\sigma _i^z),$$
(2)
with the transition frequency $`\omega _0`$. The radiation-field Hamiltonian

$$\widehat{H}_f=\frac{1}{8\pi }\left(𝐄^2+𝐇^2\right)𝑑𝐫,$$
(3)
with electric field $`𝐄`$ and magnetic field $`𝐇=\times 𝐀`$. The vector potential $`𝐀`$ is assumed to satisfy the Coulomb gauge calibration $`𝐀=0`$ and the commutation relation

$$[E^\alpha (𝐫,t),A^\beta (𝐫^{},t)]=4\pi ic\delta _{\alpha \beta }\delta (𝐫𝐫^{}),$$

where $`c`$ is the light velocity. The atom-field interaction is given by the standard dipole Hamiltonian

$$\widehat{H}_{af}=\underset{i=1}{\overset{N}{}}\left(\frac{1}{c}𝐣_i𝐀_i+𝐝_i𝐄_{0i}\right),$$
(4)
in which $`𝐀_i𝐀(𝐫_i,t)`$; the transition-current and the transition-dipole operators, respectively, are

$$𝐣_i=i\omega _0\left(𝐝\sigma _i^+𝐝^{}\sigma _i^{}\right),𝐝_i=𝐝\sigma _i^++𝐝^{}\sigma _i^{},$$
(5)
where $`𝐝`$ is the transition dipole and $`\sigma _i^\pm `$ are the ladder operators; and $`𝐄_{0i}`$ is an external field, if any. The Hamiltonian of matter, $`\widehat{H}_m`$, can be modelled by a collection of oscillators . And matter-field interactions are given by the Hamiltonian

$$\widehat{H}_{mf}=\frac{1}{c}\underset{i=1}{\overset{N_0}{}}𝐣_{mi}𝐀_i,$$
(6)
where $`𝐣_{mi}`$ is a current produced by matter at the point $`𝐫_i`$; $`N_0`$ being the number of oscillators modelling the matter.

When considering stationary properties of electromagnetic field, one usually passes to the momentum-energy representation by expanding operators in Fourier series. However this representation is not as convenient for studying dynamical properties of a system of atoms, especially when their transition wavelength is much shorter than the system size. Therefore we shall employ here the space-time representation whose brief description can be found in book  and detailed analysis in Refs. .

The Heisenberg equations of motion, with the Hamiltonian (1), yield for the atomic variables the equations

$$\frac{d\sigma _i^{}}{dt}=i\omega _0\sigma _i^{}+\left(k_0𝐝𝐀_ii𝐝𝐄_{0i}\right)\sigma _i^z,$$

$$\frac{d\sigma _i^z}{dt}=2\left(k_0𝐝𝐀_ii𝐝𝐄_{0i}\right)\sigma _i^+2\left(k_0𝐝^{}𝐀_i+i𝐝^{}𝐄_{0i}\right)\sigma _i^{},$$
(7)
where $`k_0\omega _0/c`$. The Heisenberg equations for the electric field and vector potential lead to the Maxwell operator equations, combining which with the Coulomb gauge calibration, one gets

$$\left(^2\frac{1}{c^2}\frac{^2}{t^2}\right)𝐀=\frac{4\pi }{c}𝐉,$$
(8)
with the density of current

$$𝐉(𝐫,t)=\underset{i=1}{\overset{N}{}}𝐣_i(t)\delta (𝐫𝐫_i)+\underset{i=1}{\overset{N_0}{}}𝐣_{mi}(t)\delta (𝐫𝐫_i).$$
(9)
The formal operator solution to Eq. (8) reads

$$𝐀(𝐫,t)=𝐀_{vac}(𝐫,t)+\frac{1}{c}𝐉(𝐫^{},t\frac{1}{c}|𝐫𝐫^{}|)\frac{d𝐫^{}}{|𝐫𝐫^{}|},$$
(10)
where $`𝐀_{vac}`$ is a solution of the related uniform equation. Substituting the density of current (9) into Eq. (10) results in the vector potential

$$𝐀=𝐀_{vac}+𝐀_{rad}+𝐀_{mat}$$
(11)
consisting of the potential $`𝐀_{vac}`$ of vacuum fluctuations, potential $`𝐀_{rad}`$ produced by radiating atoms, and the potential $`𝐀_{mat}`$ induced by the medium,

$$𝐀_{rad}(𝐫_i,t)=\underset{j}{\overset{N}{}}\frac{1}{cr_{ij}}𝐣_j\left(t\frac{r_{ij}}{c}\right),$$

$$𝐀_{mat}(𝐫_i,t)=\underset{j}{\overset{N_0}{}}\frac{1}{cr_{ij}}𝐣_{mj}\left(t\frac{r_{ij}}{c}\right),$$
(12)
where the summation excludes the terms with $`j=i`$, and the notation $`r_{ij}|𝐫_{ij}|,𝐫_{ij}𝐫_i𝐫_j`$, is used. It is worth emphasizing that the relation (10) is an exact expression immediately following from the evolution equation (8) and identical to the latter. No far-zone approximation is involved here. This approximation could be done  by expanding $`|𝐫𝐫^{}|^1`$ in the integrand of Eq. (10) in powers of $`r^{}/r`$, assuming that $`r|𝐫|r^{}|𝐫^{}|`$. But we do not invoke here such an expansion. Substituting Eqs. (10) to (12) into the Hamiltonian (4) and the evolution equations (7), it is easy to notice that all these equations contain an effective dipole-dipole coupling of atoms responsible for near-zone local-field effects . Expression (10) is not a final solution for the observable field but it is an operator relation having sense only in the combination with the Heisenberg equations (7). The system of Eqs. (7) and (10) follows directly from the initial Hamiltonian (1) and is often employed in optics .

Let us introduce the quantities

$$u_i(t)<\sigma _i^{}(t)>,s_i(t)<\sigma _i^z(t)>,$$
(13)
in which the angle brackets imply the statistical averaging over the atomic degrees of freedom only, not touching the vacuum and matter degrees of freedom. For the double correlators, we shall employ the mean-field-type decoupling

$$<\sigma _i^\alpha \sigma _j^\beta >=<\sigma _i^\alpha ><\sigma _j^\beta >(ij).$$
(14)
It is worth emphasizing that, since only the atomic degrees of freedom have been involved in the averaging, the decoupling (14) does not kill quantum effects related to the vacuum and the matter degrees of freedom, which have not yet been averaged out. The variables corresponding to these degrees of freedom will enter the equations of motion for the quantities (13) as operator variables or they can be modelled by stochastic variables . This is why the decoupling (14) may be called the stochastic mean-field approximation.

Another problem concerns the way of treating the retardation appearing in Eqs. (12). We cannot use the Markov approximation just ignoring this retardation since this would eliminate all local-field effects that are so important for realistic many-atom systems with wavelength shorter than the system size . But we may employ the Born approximation

$$<\sigma _j^{}\left(t\frac{r_{ij}}{c}\right)>=u_j(t)\mathrm{exp}(ik_0r_{ij})$$
(15)
retaining the account of local fields. Recall that considering resonance atoms, one always keeps in mind that the intensities of interactions of atoms with all fields are assumed to be much smaller than the transition frequency $`\omega _0`$. If this were not so, then the notion of resonance atoms as such would loose its sense, because strong fields, essentially shifting the atomic energy levels, would completely change the classification of these levels and the values of transition frequencies. Under the assumption that the transition frequency $`\omega _0`$ is the largest energy scale, from Eqs. (7) it immediately follows that the approximation (15) is completely justified \[20–22\]. This can also be called the quasirelativistic approximation since in the nonrelativistic limit, when $`c\mathrm{}`$, Eq. (15) becomes an identity.

The equations for the quantities (13) follow from Eqs. (7) where we introduce the notation for the effective field

$$f_i(t)k_0<𝐝𝐀_{rad}(𝐫_i,t)>i𝐝𝐄_{0i}+\xi _i(t),$$
(16)
in which the last term is the random field

$$\xi _i(t)k_0𝐝\left[𝐀_{vac}(𝐫_i,t)+𝐀_{mat}(𝐫_i,t)\right].$$
(17)
As usual, we also need to take into account the level width $`\gamma _1`$ and the line width $`\gamma _2`$. The introduction of these widths is done in the standard way  keeping in mind that the values of $`\gamma _1`$ and $`\gamma _2`$, in general, differ from the emission rate $`\gamma `$ of a free atom. This is because the longitudinal relaxation rate $`\gamma _1`$ is assumed to effectively include the influence of matter which atoms are dopped in, for instance, the influence of lattice vibrations. The transverse relaxation rate $`\gamma _2`$ differs from $`\gamma `$ because of the same reason of incorporating the influence of matter, for example, by taking account of dynamic and, partially, of inhomogeneous broadening. Thus, $`\gamma _1`$ and $`\gamma _2`$ are to be treated as phenomenological parameters introduced in the commonly accepted way . In this way, we obtain the equations for the average transition amplitude

$$\frac{du_i}{dt}=(i\omega _0+\gamma _2)u_i+f_is_i,$$
(18)
and the average population difference

$$\frac{ds_i}{dt}=2(u_i^{}f_i+f_i^{}u_i)\gamma _1(s_is_0).$$
(19)
The last term in Eq. (19) allows for the suppression of spontaneous emission in the case of a single atom. Really, to return to the latter case, we have to put $`f_i=0`$. Then we get the equation for a single atom, with the solution $`s_i(t)=s_0`$, as is discussed in the Introduction. But, in the case of an ensemble of atoms, the presence of nonlinear terms related to the effective field (16) can make the solutions to the evolution equations quite nontrivial. These equations are to be complimented by an equation either for $`u_i^{}`$ or for $`|u_i|^2u_i^{}u_i`$. For the latter equation, we find

$$\frac{d|u_i|^2}{dt}=2\gamma _2|u_i|^2+s_i(u_i^{}f_i+f_i^{}u_i).$$
(20)
Equations (18) to (20) are the main evolution equations we shall consider. These are the stochastic differential equations since they contain the random field (17) including the vacuum and matter degrees of freedom. Electromagnetic fluctuations coupled with collective excitations of a medium are called polaritons . Therefore Eq. (17) represents the polariton filed acting on atoms.

Strictly speaking, one should distinguish two different physical cases, depending on the relation between the longitudinal relaxation time $`T_1\gamma _1^1`$ and the characteristic time of nonresonant processes $`T_n=\gamma ^1(1+\mathrm{\Delta }^2/\gamma ^2)`$ discussed in the Introduction. As is explained in the latter, it makes sense of talking about the suppressed spontaneous emission of an atom only for times $`tT_n`$. This suppression can be either dynamic, when $`T_1T_n`$, or static, when $`T_1T_n`$. In the first case, one has to retain the parameter $`\gamma _1`$ in the evolution equations, while in the second case, one has to set $`\gamma _1=0`$. For generality, we shall consider both these situations starting with the more general case of finite $`\gamma _1`$, from which it is easy to pass to the particular case of $`\gamma _1=0`$.

## 3 Scale Separation

In what follows, we consider the case when an external field is used only for exciting atoms at the initial time, but after this the field is switched off, $`𝐄_{0i}=0`$. It is worth saying a few words about the possibility of exciting resonance atoms whose transition frequency lies in the polariton band gap. There can be several ways of doing this. The first way is by exciting the atoms with a short pulse of a resonance external field. The amplitude of the latter is to be sufficiently large since a weak electromagnetic wave with a frequency within the polariton gap cannot propagate through the medium because of total reflection. However, when the incident field is strong enough, a monocromatic wave can penetrate into the medium with a polariton gap due to nonlinear effects . Another possibility for exciting the atoms in a polariton gap could be through a third level with the related transition frequency lying outside the gap. It could also be possible to get the population inversion at the desired level through a set of levels outside the gap, by using nonresonance pumping.

It is convenient to define the arithmetic averages

$$u(t)\frac{1}{N}\underset{i=1}{\overset{N}{}}u_i(t),s(t)\frac{1}{N}\underset{i=1}{\overset{N}{}}s_i(t)$$
(21)
of the variables (13) and also the average polariton field

$$\xi (t)\frac{1}{N}\underset{i=1}{\overset{N}{}}\xi _i(t).$$
(22)
Let us introduce the collective frequency $`\mathrm{\Omega }`$ and the collective width $`\mathrm{\Gamma }`$, respectively,

$$\mathrm{\Omega }\omega _0+\gamma _2g^{}s,\mathrm{\Gamma }\gamma _2(1gs),$$
(23)
where the effective atomic coupling parameters are

$$g\frac{3\gamma }{4\gamma _2N}\underset{ij}{\overset{N}{}}\frac{\mathrm{sin}k_0r_{ij}}{k_0r_{ij}},g^{}\frac{3\gamma }{4\gamma _2N}\underset{ij}{\overset{N}{}}\frac{\mathrm{cos}k_0r_{ij}}{k_0r_{ij}},$$
(24)
and the notation

$$\gamma \frac{4}{3}k_0^3d_0^2,k_0\frac{\omega _0}{c},𝐝d_0𝐞_d$$

is employed. Then from Eqs. (18) to (20) in the mean-field approximation, we obtain

$$\frac{du}{dt}=\left(i\mathrm{\Omega }+\mathrm{\Gamma }\right)u+s\xi +\gamma _2(g+ig^{})su^{}𝐞_d^2,$$
(25)
$$\frac{ds}{dt}=4\gamma _2g|u|^22(u^{}\xi +\xi ^{}u)\gamma _1(ss_0)$$

$$2\gamma _2\left[(g+ig^{})(u^{}𝐞_d)^2+(gig^{})(𝐞_d^{}u)^2\right],$$
(26)
$$\frac{d|u|^2}{dt}=2\mathrm{\Gamma }|u|^2+s(u^{}\xi +\xi ^{}u)+\gamma _2s\left[(g+ig^{})(u^{}𝐞_d)^2+(gig^{})(𝐞_d^{}u)^2\right].$$
(27)
Equations (25) to (27) form a system of nonlinear stochastic differential equations. It is worth stressing that the rotating wave approximation (RWA), well-known in resonance optics , cannot be blindly applied for simplifying this system of equations. Really, the standard usage of RWA is as follows. One assumes that the polarization variable $`u`$ behaves approximately as $`\mathrm{exp}(i\mathrm{\Omega }t)`$, that is, rotates with the frequency $`\mathrm{\Omega }`$, while the population difference $`s`$ does not rotate. With this assumption, one neglects in the right-hand sides of the evolution equations the terms that behave qualitatively differently to the assumed behaviour of the left-hand sides. Thus, in Eq. (25) one should drop the counter-rotating term containing $`u^{}`$ as well as the non-rotating term $`s\xi `$. Similarly, in Eqs. (26) and (27) for non-rotating quantities, one should neglect the rotating terms with $`u^2`$ and $`(u^{})^2`$ as well as with $`u\xi ^{}`$ and $`\xi ^{}u`$. But, as is evident, such a procedure would completely loose information of polariton degrees of freedom $`\xi `$. Hence, RWA is not applicable to such stochastic differential equations. Although the RWA reasoning is not justified for these equations, one may try to simplify the latter by formally omitting only those terms that would be dropped if the polariton field $`\xi `$ were absent, but leaving untouched the terms containing $`\xi `$. Then one would come to the system of equations

$$\frac{du}{dt}=(i\mathrm{\Omega }+\mathrm{\Gamma })u+s\xi ,$$

$$\frac{ds}{dt}=4\gamma _2g|u|^22(u^{}\xi +\xi ^{}u)\gamma _1(ss_0),$$

$$\frac{d|u|^2}{dt}=2\mathrm{\Gamma }|u|^2+s(u^{}\xi +\xi ^{}u).$$

This is yet a system of nonlinear, because of Eqs. (23), stochastic differential equations.

To solve this system of equations (25) to (27), we shall use the scale separation approach \[30–32\] which is a variant of the Krylov-Bogolubov averaging method  generalized to the case of stochastic differential equations. To this end, we need, first of all, to classify the functional variables of the system (25) to (27) onto fast and slow. This can be done by accepting the usual inequalities

$$\frac{\gamma _1}{\mathrm{\Omega }}1,\frac{\gamma _2}{\mathrm{\Omega }}1,\left|\frac{\mathrm{\Gamma }}{\mathrm{\Omega }}\right|1$$
(28)
and assuming that the average energy of atomic interactions with the polariton field is also much less than $`\mathrm{\Omega }`$. Then from Eqs. (25) to (27) it immediately follows that the variable $`u`$ is to be classified as fast as compared to the slow variables $`s`$ and $`|u|^2`$. The next step is to solve Eq. (25) for the fast variable treating the slow variables as quasi-invariants, which yields

$$u(t)=\left[u_0+s_0^te^{(i\mathrm{\Omega }+\mathrm{\Gamma })t^{}}\xi (t^{})𝑑t^{}\right]e^{(i\mathrm{\Omega }+\mathrm{\Gamma })t},$$
(29)
where $`u_0u(0)`$. This solution is to be substituted into the right-hand sides of the equations (26) and (27) for the slow variables, with averaging these right-hand sides over polariton degrees of freedom and over explicitly entering time. In accomplishing this procedure, we define the effective atom-polariton coupling

$$\alpha \frac{\mathrm{Re}}{s\mathrm{\Gamma }}\underset{\tau \mathrm{}}{lim}\frac{1}{\tau }_0^\tau u^{}(t)\xi (t)dt,$$
(30)
where $`\mathrm{Re}`$ means the real part and $`\mathrm{}`$ implies the averaging over polariton degrees of freedom. With the form (29) and the condition $`\xi (t)=0`$, this gives

$$\alpha =\frac{\mathrm{Re}}{\mathrm{\Gamma }}\underset{\tau \mathrm{}}{lim}\frac{1}{\tau }_0^\tau 𝑑t_0^te^{(i\mathrm{\Omega }+\mathrm{\Gamma })(tt^{})}\xi ^{}(t)\xi (t^{})dt^{}.$$

By assumption, the atom-polariton coupling is weak,

$$|\alpha |1.$$
(31)
In order to understand the structure of the coupling (30), let us model the polariton field by an ensemble of $`N_0`$ oscillators, so that

$$\xi (t)=\frac{1}{\sqrt{N_0}}\underset{k}{}\gamma _k\left(b_ke^{i\omega _kt}+b_k^{}e^{i\omega _kt}\right),$$
(32)
where $`b_k`$ and $`b_k^{}`$ are Bose operators with the properties

$$b_k^{}b_k^{}=\delta _{kk^{}}n_k,b_kb_k^{}=0,$$

$$b_kb_k^{}^{}=\delta _{kk^{}}(1+n_k),\frac{1}{N_0}\underset{k}{}n_k=1,$$

and the effective width $`\gamma _k`$ is of order $`\gamma `$ outside the polariton gap but is zero inside this gap,

$$\gamma _k=0,\omega _1<\omega _k<\omega _2,\mathrm{\Delta }_p\omega _2\omega _1.$$

For the polariton field (32), the coupling (30) becomes

$$\alpha =\frac{1}{N_0}\underset{k}{}|\gamma _k|^2\left[\frac{n_k}{(\omega _k\mathrm{\Omega })^2+\mathrm{\Gamma }^2}+\frac{1+n_k}{(\omega _k+\mathrm{\Omega })^2+\mathrm{\Gamma }^2}\right].$$
(33)
We assume that the collective atomic frequency $`\mathrm{\Omega }`$ lies deeply inside the polariton gap of width $`\mathrm{\Delta }_p`$. If this width is much smaller than the frequencies at the edges of the gap, $`\mathrm{\Delta }_p\omega _1`$, then Eq. (33) can be approximated by the form

$$\alpha \frac{4\gamma ^2}{\mathrm{\Delta }_p^2+4\mathrm{\Gamma }^2}.$$

After substituting Eq. (29) into the right-hand sides of Eqs. (26) and (27), with accomplishing the described averaging and introducing the new function

$$w|u|^2\alpha s^2,$$
(34)
we obtain the equations

$$\frac{ds}{dt}=4g\gamma _2w\gamma _1(ss_0),\frac{dw}{dt}=2\gamma _2(1gs)w$$
(35)
for the slow variables.

Let us note that the system of three equations (25) to (27) has been derived from the set of $`3N`$ equations (18) to (20), where $`i=1,2,\mathrm{},N`$, by employing a mean-field type approximation for the averages (21). As is obvious, it is impossible to deal with a system of $`3N`$ equations for realistically large $`N\mathrm{}`$, because of which one always has to invoke some approximation permitting one to reduce the system of untreatably large number of equations to a treatable finite dimensional dynamical system. The most common such a way of reduction is by using the uniform approximation which assumes that the variables $`u_i`$ and $`s_i`$ do not depend on the index $`i`$. This approximation is known  to give a good description when there are no electromagnetic spatial structures and the considered sample is sufficiently large for boundary effects being neglected. The uniform approximation reduces Eqs. (18) to (20) to the same system of equations (25) to (27). Therefore, the usage of the averages (21), together with the mean-field approximation, is mathematically identical to the uniform approximation. However, from the physical point of view, dealing with the averages is more preferable since an averaging description serves as a reasonable first approximation even when the sample is nonuniform .

## 4 Liberation of Light

Suppose that at the initial time the atoms are excited so that their average population difference $`s_0s(0)>0`$. For a single atom, when $`g=0`$, from Eqs. (35) one gets $`s(t)=s_0`$, that is, the emission is suppressed. For a collective of atoms, when $`g0`$, the system of equations (35) becomes nonlinear, and its solutions essentially depend on the value of the atom-atom coupling parameter $`g`$, defined in Eq. (24).

According to the notion of suppressed emission, discussed in the Introduction, we need, first, to consider the stationary solutions of the evolution equations. For Eqs. (35), there are two stationary solutions. One of them,

$$s_1^{}=s_0,w_1^{}=0,$$
(36)
tells that emission is suppressed. And another one,

$$s_2^{}=\frac{1}{g},w_2^{}=\frac{\gamma _1(gs_01)}{4\gamma _2g^2},$$
(37)
shows that light can, at least partially, be liberated. The asymptotic stability of the stationary solutions (36) and (37) can be studied involving the Lyapunov analysis. This is done by calculating the Jacobian matrix associated with Eqs. (35) and finding its eigenvalues that are

$$\lambda ^\pm =\frac{1}{2}[\gamma _1+\gamma _2(1gs)]\frac{1}{2}\left\{\left[\gamma _12\gamma _2(1gs)\right]^232\gamma _2^2g^2w\right\}^{1/2}.$$

These eigenvalues, evaluated at the corresponding fixed points (36) and (37), define the characteristic exponents

$$\lambda _1^+=\gamma _1,\lambda _1^{}=2\gamma _2(1gs_0),$$

$$\lambda _2^\pm =\frac{\gamma _1}{2}\left\{1\pm \left[1+8\frac{\gamma _2}{\gamma _1}(1gs_0)\right]^{1/2}\right\}.$$
(38)
The real parts of the expressions in Eq. (38) are the Lyapunov exponents whose signs characterize the stability of the corresponding fixed points.

The Lyapunov analysis shows that if $`gs_0<1`$, the stationary solution (36) is a stable node, while the fixed point (37) is a saddle point. When $`gs_0=1`$, both fixed points, (36) and (37), merge together becoming neutral. In the interval

$$1<gs_01+\frac{\gamma _1}{8\gamma _2},$$

solution (36) is a saddle point, while the fixed point (37) is a stable node. When

$$gs_0>1+\frac{\gamma _1}{8\gamma _2},$$

the fixed point (36) remains a saddle point, and solution (37) becomes a stable focus with the characteristic exponents

$$\lambda _2^\pm =\frac{\gamma _1}{2}i\omega _{osc},$$

where the oscillation frequency is

$$\omega _{osc}\frac{\gamma _1}{2}\left[8\frac{\gamma _2}{\gamma _1}(gs_01)1\right]^{1/2}.$$

This frequency defines the asymptotic oscillation period

$$T_{osc}\frac{2\pi }{\omega _{osc}}=\frac{4\pi }{\sqrt{8(gs_01)\gamma _1\gamma _2\gamma _1^2}}.$$
(39)
In the case when $`\gamma _1\gamma _2,s_01`$, and $`g1`$, the latter simplifies to

$$T_{osc}2\pi \sqrt{\frac{T_1T_2}{2g}},$$

where $`\gamma _1T_11`$ and $`\gamma _2T_21`$, so that $`T_1`$ and $`T_2`$ are the longitudinal and transverse relaxation times, respectively.

In this way, if the coherent interactions between atoms are weak, and the atom-atom coupling parameter $`g`$ is small, so that $`gs_0<1`$, then the system tends to the stationary solution (36), that is, there is no liberation of light. When the atomic interactions become stronger, so that $`gs_0>1`$, the system tends to another stable stationary solution (37), and a partial liberation of light occurs since

$$s_2^{}=\frac{1}{g}<s_0(gs_0>1).$$

The portion of excitation that remains in the atomic ensemble decreases, with increasing $`g`$, as $`s_2^{}0`$. The qualitative change of the asymptotic behaviour of solutions to Eqs. (35) happens when $`gs_0=1`$. This equality defines the bifurcation point for Eqs. (35), when the dynamical system is structurally unstable . This bifurcation point separates the regions where emission remains suppressed ($`gs_0<1`$) and where light becomes partially liberated $`(gs_0>1)`$. For sufficiently strong atomic interactions, such that $`gs_0>1+\gamma _1/8\gamma _2`$, the liberation of light occurs by means of a series of pulses that, at asymptotically large times, are separated one from another by the period (39).

To describe the coherent pulse, occurring when $`\gamma _1=0`$, we may consider Eqs. (35) omitting there the relaxation term with $`\gamma _1`$. Then this system of nonlinear equations can be solved exactly resulting in the solution

$$s=\frac{\gamma _0}{g\gamma _2}\mathrm{tanh}\left(\frac{tt_0}{\tau _0}\right)+\frac{1}{g},w=\frac{\gamma _0^2}{4g^2\gamma _2^2}\mathrm{sech}^2\left(\frac{tt_0}{\tau _0}\right),$$
(40)
where the radiation width $`\gamma _0`$ is given by the relation

$$\gamma _0^2=\mathrm{\Gamma }_0^2+4g^2\gamma _2^2\left(|u_0|^2\alpha _0s_0^2\right)$$
(41)
in which $`u_0u(0)`$, $`\alpha _0\alpha (0)`$, and

$$\mathrm{\Gamma }_0\gamma _2(1gs_0),\gamma _0\frac{1}{\tau _0},$$

and where the delay time is

$$t_0=\frac{\tau _0}{2}\mathrm{ln}\left|\frac{\gamma _0\mathrm{\Gamma }_0}{\gamma _0+\mathrm{\Gamma }_0}\right|.$$
(42)
Introducing the critical atom-polariton coupling

$$\alpha _c\frac{(gs_01)^2}{4g^2s_0^2}+\frac{|u_0|^2}{s_0^2},$$
(43)
the radiation width can be written as

$$\gamma _0=2g|s_0|\gamma _2\sqrt{\alpha _c\alpha _0}.$$
(44)
The value (43) is termed critical since the coupling $`\alpha _0`$ cannot exceed $`\alpha _c`$ for the solution (40) to remain finite. The restriction $`\alpha _0\alpha _c`$ specifies condition (31) assumed for the validity of the averaging method. This restriction is not as severe as far as $`\alpha 1`$ while, for large $`g`$, $`\alpha _c1/4`$.

At the time $`t=t_0`$, the solutions (40) and (34) yield

$$s(t_0)=\frac{1}{g},w(t_0)=s_0^2(\alpha _c\alpha _0),|u(t_0)|^2=s_0^2(\alpha _c\alpha _0)+\frac{\alpha }{g^2}.$$
(45)
Then the system of atoms achieves the maximum of coherence.

Consider more in detail the case when the system of atoms is initially completely inverted, $`s_0=1`$, there is no triggering pulse, $`u_0=0`$, and the atom-atom coupling is strong, $`g>1`$. Then the critical value (43) becomes

$$\alpha _c=\left(\frac{g1}{2g}\right)^2.$$
(46)
If $`\alpha \alpha _c`$, the radiation width (44) gives

$$\gamma _0=(g1)\gamma _2\left[1\frac{2g^2\alpha _0}{(g1)^2}\right].$$
(47)
The corresponding radiation time is

$$\tau _0=\frac{T_2}{g1}\left[1+\frac{2g^2\alpha _0}{(g1)^2}\right].$$
(48)
And for the delay time (42), we have

$$t_0=\frac{T_2}{2(g1)}\mathrm{ln}\left|\frac{(g1)^2g^2\alpha _0}{g^2\alpha _0}\right|.$$
(49)
When the atomic coupling is strong, then

$$\alpha _c\frac{1}{4}(u_0=0,g1).$$
(50)
The radiation time (48) for $`g1`$ becomes

$$\tau _0\frac{T_2}{g}(1+2\alpha _0).$$
(51)
Hence, $`\tau _0`$ can be much smaller than $`T_2`$. Since $`gN`$, we have $`\tau _0N^1`$, which is typical of superradiance . Under condition (50), the delay time (49) is

$$t_0\frac{T_2}{2g}|\mathrm{ln}\alpha _0|.$$
(52)
Note that if $`\alpha _00`$, then $`t_0\mathrm{}`$, which means that emission would be suppressed for very long time. The superradiant burst develops at the delay time $`t_0`$ and lasts during the radiation time $`\tau _0`$.

In the case when $`\alpha _c\alpha _01`$, the radiation time increases, as compared to the opposite case when $`\alpha _0\alpha _c`$, being

$$\tau _0=\frac{T_2}{2g|s_0|\sqrt{\alpha _c\alpha _0}},$$
(53)
while the delay time shortens becoming

$$t_0=\frac{2g|s_0|T_2}{(gs_01)^2}\sqrt{\alpha _c\alpha _0}.$$
(54)
In this case, the radiation can be coherent if $`\tau _0T_2`$, which requires $`g\sqrt{\alpha _c\alpha _0}1`$.

The superradiant character of emission is connected with the radiation time being inversely proportional to the number of radiators, $`\tau _0N^1`$. Another characteristic describing the level of coherence in the radiating system is the radiation intensity. The total intensity of radiation, averaged over fast oscillations,

$$I(t)=I_{inc}(t)+I_{coh}(t),$$
(55)
consists of two terms, the intensity of incoherent radiation

$$I_{inc}=\frac{1}{2}\omega _0\gamma N(1+s)$$
(56)
and the intensity of coherent radiation

$$I_{coh}=\omega _0\gamma \phi _sN^2|u|^2,$$
(57)
where $`\phi _s`$ is the shape factor  given by the integral over spherical angles,

$$\phi _s=\frac{3}{8\pi }\left|𝐧\times 𝐞_d\right|^2F(k_0𝐧)𝑑\mathrm{\Omega }(𝐧),F(𝐤)\left|\frac{1}{N}\underset{i=1}{\overset{N}{}}e^{i𝐤𝐫_i}\right|^2.$$
(58)
For the point-like system, for which $`k_00`$, we have $`\phi _s1`$, and the intensity of coherent radiation is proportional to the number of radiators squared, $`I_{coh}N^2`$, which is typical of the Dicke model. For a finite-size system, the shape factor essentially depends on the relation between the radiation wavelength $`\lambda `$ and the characteristic sizes of the sample . Thus, for a cylindrical sample of radius $`R`$ and length $`L`$, we have

$`\phi _s\{\begin{array}{ccc}\frac{3\lambda }{8L},& \frac{\lambda }{2\pi L}1,& \frac{R}{L}1,\\ \frac{3}{8}\left(\frac{\lambda }{\pi R}\right)^2,& \frac{\lambda }{2\pi R}1,& \frac{L}{R}1\end{array}`$
for pencil-like or disk-like shapes, respectively. Hence for the first case $`I_{coh}N^{5/3}`$, while for the second case $`I_{coh}N^{4/3}`$.

To compare the intensities of coherent and incoherent radiation, it is convenient to introduce the coherence coefficient  defined as the ratio

$$C_{coh}(t)\frac{I_{coh}(t)}{I_{inc}(t)}.$$
(61)
For the coherent radiation intensity (57), using the relation (34), we get

$$I_{coh}(t)=\omega _0\gamma \phi _sN^2(w+\alpha s^2).$$

Therefore, the coherence coefficient (59) is

$$C_{coh}=2\phi _sN\frac{w+\alpha s^2}{1+s}.$$
(62)
At the moment of the maximal coherence of the superradiant burst, according to Eqs. (45), we find

$$C_{coh}(t_0)=2g\phi _ss_0^2N\frac{\alpha _c\alpha _0}{1+g}.$$

If $`C_{coh}>1`$, the radiation is predominantly coherent. For the case when $`u_0=0,s_0=1`$, and $`g1`$, the latter expression gives

$$C_{coh}(t_0)\frac{1}{2}\phi _sN.$$

If the sample has a pencil-like or disk-lake shape, then

$`C_{coh}(t_0)\{\begin{array}{cc}\frac{3\pi }{16}\left(\frac{R}{\lambda }\right)^2\rho \lambda ^3,& pencil\\ \frac{3}{16\pi }\left(\frac{L}{\lambda }\right)\rho \lambda ^3,& disk,\end{array}`$
where $`\rho `$ is the density of resonance atoms. When this density is sufficiently high, and $`\lambda R`$ or $`\lambda L`$, the coherence coefficient can be very large, which would mean that the radiation is almost purely coherent.

The general condition for light to be, at least partially, liberated is

$$s(\mathrm{})<s_0.$$
(65)
From equations (40), we have

$$s(\mathrm{})=\frac{\gamma _0}{g\gamma _2}+\frac{1}{g},w(\mathrm{})=0,$$

which, with the radiation width (44), gives

$$s(\mathrm{})=2|s_0|\sqrt{\alpha _c\alpha _0}+\frac{1}{g}.$$
(66)
Therefore the liberation condition (61) becomes

$$g(s_0+2|s_0|\sqrt{\alpha _c\alpha _0})>1.$$
(67)
In the case, when $`u_0=0`$, and taking into account that $`\alpha _01`$, equation (62) reduces to

$`s(\mathrm{})=\{\begin{array}{cc}s_0,& gs_0<1\\ s_0+2/g,& gs_0>1.\end{array}`$ (70)
The condition (63) simplifies to

$$gs_0>1.$$
(71)
It is interesting that, although the limit values $`s(\mathrm{})`$ are different for the cases when $`\gamma _1`$ is finite or when it is zero, but the liberation condition (65) remains the same. Under this condition, a system of atoms can radiate, though spontaneous emission of a single atom is suppressed. The case of a single atom can be recovered by setting $`g0`$. Then, as is obvious, condition (65) can never hold true. It is only when the density of doped atoms is sufficiently high, so that the atom-atom coupling $`g`$ becomes sufficiently large, satisfying condition (65), the radiation of atoms is possible being due to collective effects.

## 5 Discussion

We considered a system of resonance atoms doped into a medium with a polariton band gap. Spontaneous emission of a single atom with the transition frequency inside the polariton gap is suppressed. However an ensemble of atoms with their transition frequencies in the gap can radiate due to effective coherent interactions between the atoms. If this intertaction is sufficiently strong, light is partially liberated. The collective liberation of light occurs through one or a series of superradiant pulses.

The dynamics of light liberation has been analysed employing the scale separation approach \[30–32\], which is a generalization of the averaging method  to stochastic differential equations. This approach provides an efficient tool for treating complicated systems of nonlinear evolution equations, as has earlier been demonstrated for the problems of superradiant spin relaxation in nonequilibrium magnets \[30–32\] and of nonlinear dynamics of atoms in magnetic traps \[37–39\], where nonlinear phenomena are of crucial importance. Collective liberation of light is also a principally nonlinear phenomenon, with nonlinearity caused by coherent atomic interactions.

In conclusion, let us give some estimates for the characteristic parameters considered in the text. The polariton effect is well developed in many dielectrics and semiconductors . For instance, it is intensively studied in such semiconductors as CuCl, CuBr, CdSe, ZnSe, GaAs, GaSb, InAs, AlAs, and SiC \[10–12,40\]. The polariton gap in such materials develops around the frequency $`10^{14}`$ s-1, with the gap width $`\mathrm{\Delta }_p10^{13}`$ s-1. By assumption, the atomic transition frequency is inside the polariton gap, i.e. $`\omega _010^{14}`$ s-1. Hence the radiation wavelength is $`\lambda 10^3`$ cm. For the line width, one may take $`\gamma _210^9`$ s-1.

For the initially inverted atoms, with $`s_0=1`$, the critical atom-atom coupling parameter, above which the collective radiation becomes possible, is $`g_c=1`$. For the parameters $`g>g_c`$, the collective liberation of light occurs. If $`\gamma \gamma _2`$, then the atomic coupling parameters, defined in Eq. (24), are $`g\rho \lambda ^3`$ and $`g^{}\rho \lambda ^3`$, where $`\rho `$ is the density of atoms. Consequently, the critical density of atoms, providing $`g_c=1`$, is $`\rho _c\lambda ^3`$, which gives $`\rho _c10^9`$ cm-3. When $`\rho >\rho _c`$, radiation becomes possible because of the formation of an impurity band inside the polariton band. The width of the impurity band is of the order of the collective line width $`\mathrm{\Gamma }\gamma _2g\gamma \rho \lambda ^3`$. The latter becomes larger than the polariton gap, when the density $`\rho >\mathrm{\Delta }_p/\gamma _2\lambda ^3`$, that is, $`\rho >10^{13}`$ cm-3. For such densities of resonance atoms, the polariton gap can be overlapped by the impurity band. At the density $`\rho 10^{13}`$ cm-3, the atom-atom coupling is $`g10^4`$, and the collective line width is $`\mathrm{\Gamma }10^{13}`$ s-1.

With the densities $`\rho <10^{13}`$ cm-3, the atom-atom coupling $`g<10^4`$, while the atom-polariton coupling is $`\alpha \gamma ^2/\mathrm{\Delta }_p^2`$, that is $`\alpha 10^8`$. The radiation time $`\tau _0T_2/g`$ is larger than $`10^{13}`$ s but can be much smaller than $`T_210^9`$ s-1. And the delay time $`t_0\tau _0|\mathrm{ln}\alpha |`$ is an order longer than the radiation time. In this way, for the density of doped atoms $`\rho 10^910^{13}`$ cm-3, the effective atom-atom coupling is $`g110^4`$. Then the delay time is $`t_010^{12}10^8`$ s and the radiation time $`\tau _010^{13}10^9`$ s.

I am grateful to M.R. Singh for useful discussions and hospitality during my stay at the Centre for Interdisciplinary Studies in Chemical Physics at the University of Western Ontario, London, Canada, where this paper was started. I appreciate the Senior Fellowship of the University of Western Ontario. This work has also been supported by the Bogolubov-Infeld Grant of the State Agency for Atomic Energy, Poland, and by a Grant of the São Paulo State Research Foundation, Brazil.


text:
# Fidelity of Single Qubit Maps

thanks: to whom correspondence should be addressed at the Clarendon Laboratory (March 15, 2001)

## Abstract

We give a simple way of characterising the average fidelity between a unitary and a general operation on a single qubit which only involves calculating the fidelities for a few pure input states.

Consider transforming a quantum system by a desired operator, $`U`$. In practice, we may not be able to implement $`U`$ exactly but may actually apply a map $`𝒮`$ instead. In general, the map $`𝒮`$ is a superoperator. To gauge how closely $`𝒮`$ approximates $`U`$ (or *vice versa*), we require a measure of how close the output states are to each other, *i.e.* given identical input states, $`\rho _{in}`$, how does $`U[\rho _{in}]`$ compare with $`𝒮[\rho _{in}]`$. For the most general case where both output states are given by density operators, we can define their fidelity as uhlmann76 ,

$$F(\rho _1,\rho _2)=\left(\text{Tr}\left(\sqrt{\sqrt{\rho _1}\rho _2\sqrt{\rho _1}}\right)\right)^2,$$
(1)
which simplifies to

$$F(|\psi \psi |,\rho )=\text{Tr}\left(|\psi \psi |\rho \right),$$
(2)
in the case where at least one of the states is pure bruss98 . This situation arises naturally when $`U`$ is unitary and the input state, $`\rho _{in}`$, is pure, as in this case $`U[\rho _{in}]`$ will also be pure.

In the case where our system is a single qubit, the allowed operations are restricted to affine contractions of the Bloch Ball. We will consider the case where $`U`$ is a desired unitary transformation (rotation) of the Bloch ball, and $`𝒮`$ is a superoperator. It is assumed that $`𝒮`$ is a linear, trace-preserving map on the space of single qubit density operators ruskai01 . We will define $`\overline{F}`$ as the average fidelity over all pure input states,

$$\overline{F}=\frac{1}{4\pi }F_{|\psi \psi |}\text{d}\mathrm{\Omega },$$
(3)
where

$$F_{|\psi \psi |}=\text{Tr}\left(U|\psi \psi |U^{}𝒮\left[|\psi \psi |\right]\right).$$
(4)
The state $`|\psi \psi |`$ can be expressed in the basis of the Pauli spin matrices

$$\begin{array}{cc}\hfill |\psi \psi |& =\frac{1}{2}\left(\mathrm{𝟏}+\left(\begin{array}{cc}\mathrm{cos}\theta & \mathrm{sin}\theta e^{i\varphi }\\ \mathrm{sin}\theta e^{i\varphi }& \mathrm{cos}\theta \end{array}\right)\right)\hfill \\ & =\frac{\sigma _0}{2}+\mathrm{sin}\theta \mathrm{cos}\varphi \frac{\sigma _x}{2}+\mathrm{sin}\theta \mathrm{sin}\varphi \frac{\sigma _y}{2}+\mathrm{cos}\theta \frac{\sigma _z}{2}\hfill \\ & =\underset{j=0,x,y,z}{}c_j(\theta ,\varphi )\frac{\sigma _j}{2}.\hfill \end{array}$$
(5)
Hence Eq. (3) can be expressed as

$$\begin{array}{cc}\hfill \overline{F}& =\frac{1}{4\pi }_{\theta =0}^\pi _{\varphi =0}^{2\pi }\text{Tr}\left(U\underset{j}{}c_j(\theta ,\varphi )\frac{\sigma _j}{2}U^{}𝒮\left[\underset{k}{}c_k(\theta ,\varphi )\frac{\sigma _k}{2}\right]\right)\mathrm{sin}\theta \text{d}\varphi \text{d}\theta \hfill \\ & =\underset{jk}{}\left(\frac{1}{4\pi }_\theta _\varphi c_jc_k\mathrm{sin}\theta \text{d}\theta \text{d}\varphi \right)\text{Tr}\left(U\frac{\sigma _j}{2}U^{}𝒮\left[\frac{\sigma _k}{2}\right]\right)\hfill \end{array}$$
(6)
where we have used the linearity of $`U`$ and $`𝒮`$. The integrals of the coefficients, $`c_{jk}`$, are easily evaluated by symmetry: the cross terms vanish bowdrey01  leaving the simpler expression

$$\begin{array}{cc}\hfill \overline{F}& =\underset{jk}{}\left(\frac{2\delta _{j0}\delta _{k0}+\delta _{jk}}{3}\right)\text{Tr}\left(U\frac{\sigma _j}{2}U^{}𝒮\left[\frac{\sigma _k}{2}\right]\right)\hfill \\ & =\text{Tr}\left(U\frac{\sigma _0}{2}U^{}𝒮\left[\frac{\sigma _0}{2}\right]\right)+\frac{1}{3}\underset{j=x,y,z}{}\text{Tr}\left(U\frac{\sigma _j}{2}U^{}𝒮\left[\frac{\sigma _j}{2}\right]\right)\hfill \\ & =\frac{1}{2}+\frac{1}{3}\underset{j=x,y,z}{}\text{Tr}\left(U\frac{\sigma _j}{2}U^{}𝒮\left[\frac{\sigma _j}{2}\right]\right)\hfill \end{array}$$
(7)
where we have used the fact that $`𝒮[\sigma _0/2]`$ is a density matrix and thus has unit trace. In order to express the average fidelity in terms of states, and to give a more intuitive picture of the above expression, we use the substitutions,

$$\begin{array}{cc}\hfill \frac{\sigma _j}{2}& =\frac{\sigma _0+\sigma _j}{2}\frac{\sigma _0}{2}=\rho _j\rho _0\hfill \\ & =\frac{\sigma _0}{2}\frac{\sigma _0\sigma _j}{2}=\rho _0\rho _j,\hfill \end{array}$$
(8)
where $`\rho _{\pm j}`$ represents a pure state in the $`\pm j`$-direction and $`\rho _0`$ is the maximally mixed state. This gives the two equivalent expressions,

$$\overline{F}=\frac{1}{2}+\frac{1}{3}\underset{j=x,y,z}{}\text{Tr}\left(U\rho _jU^{}𝒮\left[\rho _j\right]U\rho _jU^{}𝒮\left[\rho _0\right]\right)$$
(9)
$$\overline{F}=\frac{1}{2}+\frac{1}{3}\underset{j=x,y,z}{}\text{Tr}\left(U\rho _jU^{}𝒮\left[\rho _j\right]U\rho _jU^{}𝒮\left[\rho _0\right]\right)$$
(10)
and taking their average yields,

$$\begin{array}{cc}\hfill \overline{F}& =\frac{1}{2}+\frac{1}{6}\underset{j=x,y,z}{}\text{Tr}\left(U\rho _jU^{}𝒮\left[\rho _j\right]+U\rho _jU^{}𝒮\left[\rho _j\right]U(\rho _j+\rho _j)U^{}𝒮\left[\rho _0\right]\right)\hfill \\ & =\frac{1}{2}+\frac{1}{6}\underset{j=x,y,z}{}\text{Tr}\left(U\rho _jU^{}𝒮\left[\rho _j\right]+U\rho _jU^{}𝒮\left[\rho _j\right]2U\rho _0U^{}𝒮\left[\rho _0\right]\right)\hfill \\ & =\frac{1}{2}+\frac{1}{6}\underset{j=x,y,z}{}\left(\text{Tr}\left(U\rho _jU^{}𝒮\left[\rho _j\right]\right)+\text{Tr}\left(U\rho _jU^{}𝒮\left[\rho _j\right]\right)1\right)\hfill \\ & =\frac{1}{6}\underset{j=\pm x,\pm y,\pm z}{}\left(\text{Tr}\left(U\rho _jU^{}𝒮\left[\rho _j\right]\right)\right).\hfill \end{array}$$
(11)
Hence, the fidelity of the superoperator $`𝒮`$ with the unitary operator $`U`$ can be calculated by simply averaging the fidelities of the six axial pure states on the Bloch sphere, $`\{\rho _{+x}`$, $`\rho _x`$, $`\rho _{+y}`$, $`\rho _y`$, $`\rho _{+z}`$, $`\rho _z\}`$. In the case where $`𝒮`$ is unital ($`𝒮[\rho _0]=\rho _0`$) it can be seen from Eq. 9 that this reduces to an average over only three states, $`\{\rho _{+x}`$, $`\rho _{+y}`$, $`\rho _{+z}\}`$; similarly Eq. 10 shows that the average can be taken over $`\{\rho _x`$, $`\rho _y`$, $`\rho _z\}`$.

###### Acknowledgements.

We thank E. Galvão and L. Hardy for helpful conversations. M.D.B. and A.J.S. thank EPSRC (UK) for a research fellowship. D.K.L.O thanks CESG (UK) for financial support. J.A.J. is a Royal Society University Research Fellow. This work is in part a contribution from the Oxford Centre for Molecular Sciences, which is supported by the UK EPSRC, BBSRC, and MRC.