# Quantum analysis of Rydberg atom cavity detector for dark matter axion search

## I Introduction

Search for the so-called “invisible” axions  as non-baryonic dark matter is one of the most challenging issues in particle physics and cosmology . The mass range $`m_a10^6\mathrm{eV}10^3\mathrm{eV}`$ is still open for the cosmic axions . (The light velocity is set to be $`c=1`$ throughout the present article.) As originally proposed by Sikivie , the basic idea for the dark matter axion search is to convert axions into microwave photons in a resonant cavity under a strong magnetic field via the Primakoff process. It is, however, extremely difficult to detect the cosmic axions due to their unusually weak interactions with ordinary matters. Pioneering tries with amplification-heterodyne-method were published already . An advanced experiment by the US group is currently continued, and some results have been reported, where the KSVZ axion of mass $`2.9\times 10^6\mathrm{eV}`$ to $`3.3\times 10^6\mathrm{eV}`$ is excluded at the 90% confidence level as the dark matter in the halo of our galaxy .

In this paper, we describe a quite efficient scheme for the dark matter axion search, where Rydberg atoms are utilized to detect the axion-converted microwave photons . An experimental apparatus called CARRACK I (Cosmic Axion Research with Rydberg Atoms in resonant Cavities in Kyoto) is now running to search for the dark matter axions over a 10 % mass range around $`10^5\mathrm{eV}`$. Based on the performance of CARRACK I, a new large-scale apparatus CARRACK II has been constructed recently to search for the axions over a wide range of mass . Clearly, in order to derive quantitative and rigorous results from the axion search with this type of Rydberg atom cavity detectors, it is essential to develop the quantum theoretical formulations and calculations on the interaction of cosmic axions with Rydberg atoms via the Primakoff process in the resonant cavity . Here, we present the details of these quantum calculations. This theoretical analysis actually provides important guides for the detailed design of Rydberg atom cavity detector, which will be examined separately in forthcoming papers. It is also fascinating that these quantum theoretical investigations on the axion-photon-atom interaction will be useful in viewpoint of the applications of the cavity quantum electrodynamics to fundamental researches.

This paper is organized as follows. In Sec. II, we describe the Rydberg atom cavity detector for the dark matter axion search. In the quantum theoretical point of view, this Rydberg atom cavity detector can be treated as the dynamical system of interacting oscillators with dissipation. In Sec. III, these quantum oscillators describing the photons, axions and atoms are introduced appropriately. Then, in Sec. IV, the interaction Hamiltonians are provided in terms of these quantum oscillators. In Sec. V, the quantum dynamics of interacting oscillators with dissipation is generally formulated in the Liouville picture. A series of master equations is derived to determine the time evolution of the quatum averages of the occupation numbers and higher-order correlations of these oscillators. By applying this formulation to the axion-photon-atom system in the resonant cavity, we can, in particular, calculate the number of the atoms in the upper state which are excited by absorbing the axion-converted and thermal background photons. In Sec. VI, we examine some characteristic properties of the axion-photon-atom interaction by solving analytically the master equation for the simple case with time-independent atom-photon coupling. Then, in Sec. VII, in order to make precise estimates on the sensitivity of the Rydberg atom cavity detector, we elaborate the calculations by taking into account the motion and uniform distribution of the Rydberg atoms in the incident beam and also the spatial variation of the electric field in the cavity. In Sec. VIII, the dependence of the detection sensitivity on the relevant experimental parameters is discussed. Detailed numerical calculations are performed in Sec. IX, and the quantum evolution of the axion-photon-atom system in the resonant cavity is determined precisely. The counting rates of signal and noise are calculated with the steady solutions of the master equation. The sensitivity of the Rydberg atom cavity detector is then estimated from these calculations. Finally, we summarize the present quantum analysis in Sec. X. Appendices are devoted to some supplementary issues.

## II Rydberg atom cavity detector for dark matter axion search

A schematic diagram of the Rydberg atom cavity detector (CARRACK I) is shown in Fig. 1. The axions are first converted into photons under a strong magnetic field in the conversion cavity. Then, the photons are transferred to the detection cavity via a coupling hole, and they are absorbed there by Rydberg atoms. The detection cavity is set to be free from magnetic field to avoid the complexity of the atomic energy levels due to the Zeeman splitting. The Rydberg atoms, the transition frequency of which is tuned approximately to the cavity resonant frequency $`1\mathrm{G}\mathrm{H}\mathrm{z}`$, are prepared initially in a lower state with principal quantum number $`n100`$ by exciting alkaline atoms in the ground state with laser excitation just in front of the detection cavity. It should, however, be noted that atoms in the upper state with $`n^{}(>n)`$ are not generated at this stage. The atoms prepared in this way are excited to the upper state by absorbing the microwave photons in the detection cavity, and they are detected quite efficiently with the selective field ionization method  just out of the cavity. By cooling the resonant cavity system down to about 10 mK with a dilution refrigerator in high vacuum, the thermal background photons can be reduced sufficiently to obtain a significant signal-to-noise ratio. Hence the Rydberg atom cavity detector, which is free from the amplifier noise by itself, is expected to be quite efficient for the dark matter axion search.

In the quantum theoretical point of view, the Rydberg atom cavity detector can be treated as the dynamical system of the interacting oscillators with dissipation which appropriately describe the axions, photons and Rydberg atoms. In the following sections we develop detailed quantum theoretical formulations and calculations for this dynamical system. Specifically, the analysis is made by taking into account the actual experimental situation such as the motion and uniform distribution of the Rydberg atoms in the incident beam and also the spatial variation of the electric field in the cavity. This quantum treatment provides proper estimates on the sensitivity of the Rydberg atom cavity detector for dark matter axion search.

## III Axion-photon-atom system in the cavity

We first identify the quantum oscillators which appropriately describe the photons, axions and Rydberg atoms in the resonant cavity. This provides the basis for investigating various properties of the axion-photon-atom system.

### A Resonant mode of photons

The electric field operator in the cavity $`𝒱`$ is given by

$$𝑬(𝐱,t)=(\mathrm{}\omega _c/2ϵ_0)^{1/2}[𝜶(𝐱)c(t)+𝜶^{}(𝐱)c^{}(t)]$$
(1)   
for the radiation mode with a resonant frequency $`\omega _c`$, where $`ϵ_0`$ is the dielectric constant. The creation and annihilation operators satisfy the usual commutation relation

$$[c,c^{}]=1.$$
(2)   
The mode vector field $`𝜶(𝐱)`$ is normalized by the condition

$$_𝒱|𝜶(𝐱)|^2d^3x=1.$$
(3)   
The whole cavity $`𝒱`$ may be viewed as a combination of two subcavities, the conversion cavity $`𝒱_1`$ with volume $`V_1`$ and the detection cavity $`𝒱_2`$ with volume $`V_2`$, which are coupled together:

$$𝒱=𝒱_1𝒱_2.$$
(4)   
The axion-photon conversion takes place in $`𝒱_1`$ under the strong mangetic field, while the Rydberg atoms are excited by absorbing the photons in $`𝒱_2`$. It is then suitable to divide the mode vector as

$$𝜶(𝐱)=𝜶_1(𝐱)+𝜶_2(𝐱),$$
(5)   
where $`𝜶_1(𝐱)=\mathrm{𝟎}`$ for $`𝐱𝒱_2`$ and $`𝜶_2(𝐱)=\mathrm{𝟎}`$ for $`𝐱𝒱_1`$, respectively. The normalization condition of $`𝜶(𝐱)`$ is rewritten as

$$_{𝒱_1}|𝜶_1(𝐱)|^2d^3x+_{𝒱_2}|𝜶_2(𝐱)|^2d^3x=1.$$
(6)   
The actual cavity is designed so that neglecting the small joint region the subcavities $`𝒱_1`$ and $`𝒱_2`$ admit the mode vectors $`𝜶_1^0(𝐱)`$ and $`𝜶_2^0(𝐱)`$ (up to the normalization and complex phase), respectively, whose frequencies are tuned to be almost equal to some common value $`\omega _c^0`$. In this situation, as confirmed by numerical calculations and experimental observations, two nearby eigenmodes with the frequencies $`\omega _c,\omega _c^{}\omega _c^0`$ are obtained for the whole cavity $`𝒱`$. Then, the mode vector $`𝜶(𝐱)`$ is constructed approximately of $`𝜶_1(𝐱)𝜶_1^0(𝐱)`$ and $`𝜶_2(𝐱)𝜶_2^0(𝐱)`$ with significant magnitudes in both $`𝒱_1`$ and $`𝒱_2`$. The conversion of the cosmic axions takes place predominantly to the radiation mode which is resonant with the axions satisfying the condition $`|\omega _cm_a/\mathrm{}|\gamma _a`$ (axion width) $``$ small fraction of $`\gamma `$ (cavity damping rate), as will be discussed later. The cavity can be designed so as to give a sufficient separation of $`|\omega _c\omega _c^{}|>\mathrm{several}\gamma `$ for the nearby modes with strong coupling between $`𝒱_1`$ and $`𝒱_2`$. Therefore, in the search for the signal from the cosmic axions, the one resonant mode can be extracted solely for the electric field in a good approximation, as given in Eq. (1), whose frequency $`\omega _c`$ is supposed to be close enough to the axion frequency $`\omega _a=m_a/\mathrm{}`$.

### B Coherent mode of cosmic axions

The axion field operator is expanded as usual in terms of the continuous modes:

$`\varphi (𝐱,t)`$ $`=`$ $`\mathrm{}^{1/2}{\displaystyle \frac{d^3k}{(2\pi )^32\omega _k}\left[a_𝐤(t)\mathrm{e}^{i𝐤𝐱}+a_𝐤^{}(t)\mathrm{e}^{i𝐤𝐱}\right]}`$ (7)   
$``$ $`\varphi ^+(𝐱,t)+\varphi ^{}(𝐱,t),`$ (9)   
where $`\varphi ^+`$ and $`\varphi ^{}`$ represent the positive and negative frequency parts, respectively. The cosmic axions form a coherent state with momentum distribution $`\eta _a(𝐤)`$,

$$|\eta _a=N_\eta ^{1/2}\mathrm{exp}\left(\frac{d^3k}{(2\pi )^32\omega _k}\eta _a(𝐤)a_𝐤^{}\right)|0,$$
(10)   
where $`N_\eta `$ is the normalization constant to ensure the condition $`\eta _a|\eta _a=1`$. The velocity dispersion of the galactic axions is expected to be very small as $`\beta _a10^3`$ so that $`\eta _a(𝐤)`$ takes significant values only in a small region $`_a`$ where $`|𝐤|\beta _am_a/\mathrm{}`$ .

The axion field operator (9) may be divided into the low-momentum part $`\varphi __a`$ with $`𝐤_a`$ and the residual part $`\varphi _{\mathrm{res}}`$:

$$\varphi (𝐱,t)=\varphi __a(𝐱,t)+\varphi _{\mathrm{res}}(𝐱,t).$$
(11)   
Since the coherent region $`_a`$ is chosen so as to give $`\eta _a|\varphi |\eta _a\eta _a|\varphi __a|\eta _a`$, the residual part $`\varphi _{\mathrm{res}}`$ does not provide significant contributions in the following calculations, i.e., $`\eta _a|\varphi _{\mathrm{res}}|\eta _a0`$. It is further noticed that the spatial variation of $`\varphi __a(𝐱,t)`$ is negligible in the cavity region, i.e., $`\mathrm{e}^{i𝐤𝐱}1`$ for $`|𝐤|\beta _am_a/\mathrm{}`$ and $`𝐱𝒱_1`$. This is because the de Broglie wavelength of the axions, $`\lambda _ah/(\beta _am_a)100\mathrm{m}`$ typically for $`m_a10^5\mathrm{eV}`$ and $`\beta _a10^3`$, is much longer than the microwave cavity scale $`0.1\mathrm{m}`$. In this situation to ensure approximately $`\varphi __a(𝐱,t)\varphi __a(\mathrm{𝟎},t)`$ in the cavity, the coherent axion mode can be identified as

$`a(t)`$ $`=`$ $`(\mathrm{}\mathrm{\Sigma }_a)^{1/2}\varphi __a^+(\mathrm{𝟎},t)`$ (12)   
$`=`$ $`\mathrm{\Sigma }_a^{1/2}{\displaystyle __a}{\displaystyle \frac{d^3k}{(2\pi )^32\omega _k}}a_𝐤(t).`$ (13)   
Here, the normalization factor is given by

$$\mathrm{\Sigma }_a=__a\frac{d^3k}{(2\pi )^32\omega _k}\frac{1}{2m_a}\left(\frac{\beta _am_a}{2\pi \mathrm{}}\right)^3,$$
(14)   
so that the coherent mode operator satisfies the canonical commutaion relation,

$$[a,a^{}]=1.$$
(15)   
The normalization factor $`\mathrm{\Sigma }_a`$ is expressed in terms of the mean velocity $`\beta _a`$ of galactic axions with $`\omega _km_a/\mathrm{}`$ and $`\beta _a1`$. It should, however, be realized that this parameter $`\beta _a`$ is introduced just for the normalization of the coherent axion mode. The actual mean velocity of axions is rather determined by the axion spectrum, which may slightly be different from the parameter $`\beta _a`$ given in Eq. (14). The normalization factor $`\mathrm{\Sigma }_a`$ is canceled out anyway in calculating the axion signal rate, as explicitly seen later.

The coherent cosmic axions can be described effectively in terms of a single mode oscillator, as shown in the above. Then, the energy spread of the galactic axions can be taken into account as the damping rate of the coherent axion mode,

$$\gamma _a\beta _a^2m_a/\mathrm{}.$$
(16)   
The expectation value of the number operator $`a^{}a`$ of the coherent axion mode is determined by the energy density $`\rho _a`$ of the cosmic axions as follows. The number density operator for the axion field is given by $`\widehat{n}_a(𝐱,t)=(2\mathrm{})^1(i\varphi ^{}_0\varphi ^+i_0\varphi ^{}\varphi ^+)`$, and its expectation value is calculated for $`𝐱𝒱_1`$ and $`\beta _a1`$ as

$`\eta _a|\widehat{n}_a(𝐱,0)|\eta _a`$ $``$ $`\mathrm{}^2m_a\eta _a|\varphi __a^{}(\mathrm{𝟎},0)\varphi __a^+(\mathrm{𝟎},0)|\eta _a`$ (17)   
$``$ $`\rho _a/m_a.`$ (18)   
Then, by considering Eq.(13) to relate the relevant part of axion field operator $`\varphi __a`$ to the coherent mode operators $`a`$ and $`a^{}`$, we can calculate the occupation number of the coherent axion mode as

$$\overline{n}_a=\eta _a|a^{}(0)a(0)|\eta _a\left(\frac{2\pi \mathrm{}}{\beta _am_a}\right)^3\left(\frac{\rho _a}{m_a}\right).$$
(19)   
This expression implies that the coherent axion mode is normalized suitably in a box with a volume of $`(\text{de Broglie wavelength})^3`$.

### C Oscillator for Rydberg atoms

The Rydberg atom is treated well as a two-level system in the resonant cavity. The relevant lower and upper atomic states are represented by $`|n`$ and $`|n^{}`$, respectively. They are connected by the electric dipole transition with frequency $`\omega _b=(E_n^{}E_n)/\mathrm{}`$, which is actually fine-tuned to be almost equal to the cavity frequency $`\omega _c`$ by utilizing the small Stark shift. This two-level atomic system is described in terms of spin 1/2 like operators 

$`D^+`$ $`=`$ $`|n^{}n|,D^{}=|nn^{}|,`$ (20)   
$`D^3`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(|n^{}n^{}||nn|\right).`$ (21)   
These operators satisfy the usual SU(2) Lie algebra,

$$[D^+,D^{}]=2D^3,[D^3,D^\pm ]=\pm D^\pm .$$
(22)   
The free Hamiltonian for the two-level system is then given by

$$H_{\mathrm{atom}}=\mathrm{}\omega _bD^3.$$
(23)   
The Rydberg atoms are initially prepared in the lower state in the present detectoin scheme for cosmic axions. It is also expected that the probability of the excitation of an atom to the upper state by absorbing a photon is small enough. This is because the average number of photons in the cavity is much smaller than one at sufficiently low temperatures ($`10\mathrm{m}\mathrm{K}`$). Therefore, the atoms almost remain in the lower state, providing a good approximation

$$D^3n|D^3|n=\frac{1}{2}$$
(24)   
in the first commutation relation of Eq.(22). Then, the two-level atomic system operators may be substituted by those of an oscillator as 

$$bD^{},b^{}D^+,$$
(25)   
which satisfy the commutation relation

$$[b,b^{}]=1.$$
(26)   
The second commutation relation of Eq.(22) is also reproduced with $`D^3=\frac{1}{2}+bb^{}`$ up to the higher order terms, and accordingly the free atomic Hamiltonian is expressed effectively as

$$H_{\mathrm{atom}}=\mathrm{}\omega _bb^{}b,$$
(27)   
where the irrelevant constant $`\mathrm{}\omega _b/2`$ is discarded. This treatment of the Rydberg atom in terms of the quantum oscillator is valid if the probability in the upper state (or the lower state) is very small, i.e., $`b^{}b1`$, as is the case in the present scheme.

### D Characteristics of the photons, axions and atoms

We have seen so far that the resonant photons, coherent cosmic axions and Rydberg atoms in the resonant cavity are described in terms of the appropriate quantum oscillators with dissipation due to the couplings to the relevant reservoirs. The characteristic properties of these quantum oscillators are summarized as follows.

The thermal photon number $`\overline{n}_c`$ of the resonant mode is determined by the cavity temperature $`T_c`$ as

$$\overline{n}_c=\left(\mathrm{e}^{\mathrm{}\omega _c/k_\mathrm{B}T_c}1\right)^1.$$
(28)   
The damping rate $`\gamma _c`$ of photons is described in terms of the quality factor $`Q`$ of the cavity as

$$\gamma _c\gamma =5\times 10^{10}\mathrm{eV}\mathrm{}^1\left(\frac{\mathrm{}\omega _c}{10^5\mathrm{eV}}\right)\left(\frac{2\times 10^4}{Q}\right).$$
(29)   
The coherent axion mode is normalized in a box of the de Broglie wavelength, as seen in Eq. (19). The axion number is then estimated as

$$\overline{n}_a=5.7\times 10^{25}\left(\frac{\rho _a}{0.3\mathrm{GeVcm}^3}\right)\left(\frac{10^3}{\beta _a}\right)^3\left(\frac{10^5\mathrm{eV}}{m_a}\right)^4,$$
(30)   
where the enegy density of the cosmic axions $`\rho _a`$ is taken to be equal to that of the galactic dark halo $`\rho _{\mathrm{halo}}0.3\mathrm{GeVcm}^3`$. The dissipation of the coherent axions is characterized by their energy spread as

$$\gamma _a\beta _a^2m_a/\mathrm{}=0.02\gamma \left(\frac{\beta _a^2Q}{10^6\times 2\times 10^4}\right),$$
(31)   
where the resonant condition $`m_a\mathrm{}\omega _c`$ is considered. All the atoms are prepared in the lower state so that

$$\overline{n}_b=0.$$
(32)   
The atomic dissipation is determined by its lifetime $`\tau _b`$ as $`\gamma _b=\tau _b^1`$. The lifetime of the Rydberg state is typically $`\tau _b10^3\mathrm{s}`$ for $`n100`$ in the vacuum . Since the transitions to the off-resonant states are highly suppressed in the resonant cavity, the atomic lifetime may even be longer. Hence, the damping rate of Rydberg atoms is estimated to be at most

$$\gamma _b=6.6\times 10^{13}\mathrm{eV}\mathrm{}^1\left(\frac{10^3\mathrm{s}}{\tau _b}\right),$$
(33)   
which is actually much smaller than the photon damping rate $`\gamma `$.

## IV Interactions in the resonant cavity

In this section, we provide the interaction Hamiltonians in terms of the quantum oscillators. They determine the quantum evolution of the axion-photon-atom system in the resonant cavity, which will be investigated in detail in the following sections.

### A Axion-photon interaction

The axion-photon interaction under a strong static magnetic field with flux density $`𝑩_0`$ is described by the Lagrangian density

$$_a=\mathrm{}^{1/2}ϵ_0g_{a\gamma \gamma }\varphi 𝑬𝑩_0,$$
(34)   
where $`\mathrm{}^{1/2}`$ and $`ϵ_0`$ (dielectric constant) are explicitly factored out so that the Lagrangian density has the right dimension $`_a\mathrm{}\mathrm{s}^1\mathrm{m}^3\mathrm{eVm}^3`$ with the axion-photon-photon coupling constant $`g_{a\gamma \gamma }\mathrm{eV}^1`$, as given below. The axion-photon-photon coupling constant is calculated  as

$$g_{a\gamma \gamma }=c_{a\gamma \gamma }\frac{\alpha }{2\pi ^2}\frac{m_a}{f_\pi m_\pi }\frac{(1+Z)}{\sqrt{Z}},$$
(35)   
where $`Z=m_u/m_d`$, and

$$c_{a\gamma \gamma }=\frac{E}{C}\frac{2(4+Z)}{3(1+Z)}$$
(36)   
with

$$E=\mathrm{Tr}Q_{\mathrm{PQ}}Q_{\mathrm{em}}^2,C\delta _{ab}=\mathrm{Tr}Q_{\mathrm{PQ}}\lambda _a\lambda _b.$$
(37)   
The parameter $`c_{a\gamma \gamma }`$ represents the variation of the axion-photon-photon coupling depending on the respective Peccei-Quinn models such as the so-called KSVZ  and DFSZ .

The original Lagrangian density for the axion-photon-photon coupling in Eq. (34) provides the effective interaction Hamiltonian between the coherent axion mode $`a`$ and the resonant radiation mode $`c`$,

$$H_{ac}=\mathrm{}\kappa (a^{}c+ac^{}).$$
(38)   
The axion-photon conversion in the cavity $`𝒱_1`$ is well described with this interaction Hamiltonian. The coupling constant $`\kappa `$ is determined for $`\omega _cm_a/\mathrm{}`$ by considering the relations $`|\eta _a|\varphi ^\pm |\eta _a|\mathrm{}(\rho _a/m_a^2)^{1/2}`$ from Eq. (18) and $`|\eta _a|a|\eta _a|=\overline{n}_a^{1/2}`$ from Eq. (19) in calculating $`\eta _a|_{𝒱_1}(_a)d^3x|\eta _a\eta _a|H_{ac}|\eta _a`$ :

$`\kappa `$ $`=`$ $`\mathrm{}^{1/2}g_{a\gamma \gamma }ϵ_0^{1/2}B_{\mathrm{eff}}\left[\left({\displaystyle \frac{\beta _am_a}{2\pi \mathrm{}}}\right)^3{\displaystyle \frac{V_1}{2}}\right]^{1/2}`$ (39)   
$`=`$ $`4\times 10^{26}\mathrm{eV}\mathrm{}^1\left({\displaystyle \frac{g_{a\gamma \gamma }}{1.4\times 10^{15}\mathrm{GeV}^1}}\right)\left({\displaystyle \frac{B_{\mathrm{eff}}}{4\mathrm{T}}}\right)`$ (40)   
$`\times `$ $`\left({\displaystyle \frac{\beta _am_a}{10^3\times 10^5\mathrm{eV}}}\right)^{3/2}\left({\displaystyle \frac{V_1}{5000\mathrm{c}\mathrm{m}^3}}\right)^{1/2},`$ (41)   
where

$$B_{\mathrm{eff}}=\zeta _1GB_0,$$
(42)   
and $`B_0`$ is the maximal density of the external magnetic flux. The axion-photon-photon coupling constant $`g_{a\gamma \gamma }`$ is taken here to be the value expected from the DFSZ axion model  at $`m_a=10^5\mathrm{eV}`$.

This effective axion-photon coupling $`\kappa `$ apparently involves the $`\beta _a`$ dependence coming from the normalization given in Eq. (14). This $`\beta _a`$ dependence is, however, canceled with that of $`\overline{n}_a`$ given in Eq. (19) in calculating the signal rate which is actually proportional to $`(\kappa /\gamma )^2\overline{n}_a`$. The form factor for the magnetic field is given by

$$G=\zeta _1^1V_1^{1/2}\left|_{𝒱_1}d^3x𝜶_1(𝐱)[𝑩_0(𝐱)/B_0]\right|$$
(43)   
with

$$\zeta _1=\left[_{𝒱_1}d^3x|𝜶_1(𝐱)|^2\right]^{1/2}.$$
(44)   
This additional factor $`\zeta _1`$ ($`<1`$ as seen from Eq.(6)) represents the effective reduction of the axion-photon conversion which is due to the fact that the magnetic field is applied only in the conversion cavity $`𝒱_1`$. We may obtain the effective magnetic field strength $`B_{\mathrm{eff}}4\mathrm{T}`$, as taken in Eq.(41), by using typically a magnet of $`B_07\mathrm{T}`$ and the cavity system with $`G=\sqrt{0.7}`$ of $`\mathrm{TM}_{010}`$ mode and $`\zeta _10.7`$ for the conversion cavity.

It should here be remarked that the possible interaction term $`\mathrm{}\kappa (ac+a^{}c^{})`$ is not included in Eq. (38). This is justified as follows. The characteristic time scale for the evolution of the system considered in the present scheme is placed by the lifetime of photons $`\tau _\gamma =\gamma ^1`$ in the cavity. Hence, although the interaction term $`\mathrm{}\kappa (ac+a^{}c^{})`$ representing the processes with energy change of $`\mathrm{}\omega _c+m_a`$ is also obtained from the original axion-photon-photon coupling $`_a`$, its effects are suppressed sufficiently for $`Q1`$ by the energy conservation in this time scale $`\tau _\gamma `$.

### B Atom-photon interaction

We next consider the interaction between the Rydberg atoms and the resonant photons. The Rydberg atoms are utilized for counting photons in the cavity. They are excited by absorbing the photons through the electric dipole transition with frequency $`\omega _b`$. The emission and absorption of photon by the two-level atomic system is described by the interaction Hamiltonian

$`H_\mathrm{\Omega }`$ $`=`$ $`\mathrm{}\mathrm{\Omega }[D^+c+D^{}c^{}]`$ (45)   
$``$ $`\mathrm{}\mathrm{\Omega }[b^{}c+bc^{}].`$ (46)   
Here, we may assume for simplicity, though not essential, that the mode vector field is described in terms of a uniform complex polarization vector $`\mathit{ϵ}_2`$ ($`|\mathit{ϵ}_2|=1`$) and a real profile function $`f_2(𝐱)`$ as

$$𝜶_2(𝐱)=\mathit{ϵ}_2\overline{V}_2^{1/2}f_2(𝐱)$$
(47)   
with $`\overline{V}_2|𝜶_2(𝐱)|_{\mathrm{max}}^2`$. The profile function is normalized for convenience by the condition $`|f_2(𝐱)|_{\mathrm{max}}=1`$ at the antinodal positions. Then, the intrinsic atom-photon coupling constant $`\mathrm{\Omega }`$ is evaluated at the antinodal position in terms of the electric dipole transition matrix element $`d=|\mathit{ϵ}_2𝒅|`$ projected to the direction of polarization:

$$\mathrm{\Omega }=\frac{d}{\mathrm{}}\left(\frac{\mathrm{}\omega _c}{2ϵ_0\overline{V}_2}\right)^{1/2},$$
(48)   
where $`\overline{V}_2V_2`$ is expected in accordance with the electric field normalization condition (6). The atom-photon coupling at an arbitrary atomic position in the cavity is also given with the profile function by

$$\mathrm{\Omega }(𝐱)=\mathrm{\Omega }f_2(𝐱).$$
(49)   
In the actual experimental system, the Rydberg atoms are injected as a uniform beam. Then, certain number $`N`$ of Rydberg atoms are constantly present in the cavity. Such an ensemble of identical atoms behaves as a collective system in the interaction with the resonant radiation mode. The novelty of Rydberg atom physics is to provide situations where this collective behavior appears for a relatively small number of atoms, typically $`N10^310^6`$ . Suppose for simplicity that the $`N`$ atoms are at the antinodal position ($`\mathrm{\Omega }(𝐱)=\mathrm{\Omega }`$) in the detection cavity. Then, the effective coupling between the radiation mode and the $`N`$ atoms is actually given by

$$H_{bc}=\underset{i=1}{\overset{N}{}}H_\mathrm{\Omega }(b_i,c)=\mathrm{}\mathrm{\Omega }_N(b^{}c+bc^{}),$$
(50)   
where the collective atomic mode operator is defined by

$$b\frac{1}{\sqrt{N}}\underset{i=1}{\overset{N}{}}b_i$$
(51)   
satisfying the commutation relation $`[b,b^{}]=1`$. (Hereafter we may use without confusion the operators $`b`$ and $`b^{}`$ for the collective atomic mode rather than the single atomic mode.) The collective atom-photon coupling $`\mathrm{\Omega }_N=\mathrm{\Omega }\sqrt{N}`$  is estimated typically as

$$\mathrm{\Omega }_N=1\times 10^{10}\mathrm{eV}\mathrm{}^1\left(\frac{\mathrm{\Omega }}{5\times 10^3\mathrm{s}^1}\right)\left(\frac{N}{10^3}\right)^{1/2}.$$
(52)   
Hence, if the number of Rydberg atoms is as large as $`N10^3`$, this collective coupling can be comparable to the cavity damping rate $`\gamma `$ as given in Eq. (29). When the atomic motion and distribution in the cavity are taken into account, the atom-photon coupling should be modified suitably. This point will be treated in Sec. VII.

## V Quantum evolution of the system

We here formulate the quantum dynamics of interacting oscillators with dissipation. A systematic procedure is developed in the Liouville picture to calculate the quantum averages of the particle numbers and higher-order correlations. In practice, a series of master equations is derived for such quantities from the Fokker-Planck equation based on the Liouville picture. (The quantum evolution of the system can also be described in the Langevin picture, which is considered in Appendix B.) By applying this formulation, the quatum evolution of the axion-photon-atom system in the resonant cavity is determined with the effective interaction Hamiltonians presented in the preceeding section. Then, we can calculate, in particular, the number of the atoms in the upper state which are excited by absorbing the axion-converted and thermal photons.

### A Evolution in Liouville picture

The reduced density matrix of the damped oscillators $`q_i`$ (such as $`a,b,c`$) interacting each other obeys the Liouville equation ,

$$\frac{d\rho }{dt}=\frac{1}{i\mathrm{}}[H,\rho ]+\mathrm{\Lambda }\rho .$$
(53)   
The Liouvillian relaxations are represented by the operator $`\mathrm{\Lambda }\rho `$, which may explicitly be given by

$`\mathrm{\Lambda }\rho `$ $`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle \frac{\gamma _i}{2}}\left[2q_i\rho q_i^{}q_i^{}q_i\rho \rho q_i^{}q_i\right]`$ (54)   
$`+`$ $`{\displaystyle \underset{i}{}}\gamma _i\overline{n}_i\left[q_i^{}\rho q_i+q_i\rho q_i^{}q_i^{}q_i\rho \rho q_iq_i^{}\right],`$ (55)   
where $`\gamma _i`$ and $`\overline{n}_i`$ are the damping rates and equilibrium occupation numbers, respectively. The total Hamiltonian is given by

$$H=\underset{i}{}\mathrm{}\omega _iq_i^{}q_i+\underset{ij}{}\mathrm{}\mathrm{\Omega }_{ij}(t)q_i^{}q_j,$$
(56)   
where $`\mathrm{\Omega }_{ij}(t)=\mathrm{\Omega }_{ji}^{}(t)`$ represent the interaction terms, which may be time-dependent, as is the case for the interaction between photons and moving atoms in the cavity (see Eq. (64)).

The solution of the Liouville equation is generally obtained in terms of the creation and annihilation operators as

$$\rho =\rho (q_i,q_i^{},t).$$
(57)   
The quantum average of a physical quantity represented by a relevant operator $`𝒪`$ is evaluated with this density matrix by the formula

$$\overline{𝒪}(t)=𝒪\mathrm{Tr}[𝒪(q_i,q_i^{},t)\rho (q_i,q_i^{},t)].$$
(58)   
In practical calculations, it may be useful to take the coherent state basis . Then, the density operator $`\rho `$ is represented by the “classical” time-dependent distribution function $`P(\alpha _i,\alpha _i^{},t)`$, and the quantum average is calculated by

$$\overline{𝒪}(t)=\underset{i}{}d^2\alpha _i𝒪_\mathrm{n}(\alpha _i,\alpha _i^{},t)P(\alpha _i,\alpha _i^{},t).$$
(59)   
Here the subscript “n” means the the expression of the operator $`𝒪`$ when it is written in the normal ordered form. For example, $`𝒪_\mathrm{n}(\alpha ,\alpha ^{})=\alpha ^{}\alpha +1`$ for $`𝒪(q,q^{})=qq^{}`$.

The Liouville equation is expressed in the coherent state basis as the Fokker-Planck equation to determine the distribution function $`P`$:

$`{\displaystyle \frac{P}{t}}`$ $`=`$ $`i_{ij}(t){\displaystyle \frac{}{\alpha _i}}(\alpha _jP)i_{ij}^{}(t){\displaystyle \frac{}{\alpha _i^{}}}(\alpha _j^{}P)`$ (60)   
$`+`$ $`𝒟_{ij}{\displaystyle \frac{^2P}{\alpha _i\alpha _j^{}}}.`$ (61)   
Here, the effective Hamiltonian and diffusion term for the damped oscillators are represented by the matrices 

$`_{ij}(t)`$ $`=`$ $`\left(\omega _i{\displaystyle \frac{i}{2}}\gamma _i\right)\delta _{ij}+\mathrm{\Omega }_{ij}(t)(1\delta _{ij}),`$ (62)   
$`𝒟_{ij}`$ $`=`$ $`\gamma _i\overline{n}_i\delta _{ij}.`$ (63)   
They are given explicitly for the axion-photon-atom system under consideration by

$$(t)=\left(\begin{array}{ccc}\omega _b\frac{i}{2}\gamma _b& \mathrm{\Omega }(t)& 0\\ \mathrm{\Omega }(t)& \omega _c\frac{i}{2}\gamma _c& \kappa \\ 0& \kappa & \omega _a\frac{i}{2}\gamma _a\end{array}\right)$$
(64)   
with $`\omega _a=m_a/\mathrm{}`$, and

$$𝒟=\left(\begin{array}{ccc}\gamma _b\overline{n}_b& 0& 0\\ 0& \gamma _c\overline{n}_c& 0\\ 0& 0& \gamma _a\overline{n}_a\end{array}\right).$$
(65)   
The collective atom-photon coupling becomes time-dependent through the atomic motion along the $`x`$ axis with velocity $`v`$ as

$$\mathrm{\Omega }(t)=\mathrm{\Omega }_Nf(vt).$$
(66)   
Here, it is assumed for simplicity that the $`N`$ atoms are injected together as a pulsed beam. This time-dependence is determined by the atomic velocity $`v`$ and the electric field profile $`f(x)f_2(x,y_0,z_0)`$ of the detection cavity along the atomic beam which is injected in the $`x`$ direction from the point $`(0,y_0,z_0)`$ in the $`y`$-$`z`$ plane.

### B Master equations

We may be interested in the self-adjoint multiple moments of the oscillators, 

$`𝒩_{i_1\mathrm{}i_pj_1\mathrm{}j_p}^{(2p)}(t)q_{i_1}^{}\mathrm{}q_{i_p}^{}q_{j_1}\mathrm{}q_{j_p}`$ (67)   
$`={\displaystyle \underset{i}{}d^2\alpha _i\alpha _{i_1}^{}\mathrm{}\alpha _{i_p}^{}\alpha _{j_1}\mathrm{}\alpha _{j_p}P(\alpha _i,\alpha _i^{},t)},`$ (68)   
where $`2p`$ denotes the number of involved operators. A series of master equations is then obtained for these moments by considering the Fokker-Planck equation for the distribution function:

$$\frac{d𝒩^{(2p)}}{dt}=i^{(2p)}(t)𝒩^{(2p)}+𝒟^{(2p2)}𝒩^{(2p2)}.$$
(69)   
Here the linear operators $`^{(2p)}`$ and $`𝒟^{(2p2)}`$ acting on the moments are defined by 

$`(^{(2p)}(t)𝒩^{(2p)})_{IJ}`$ $``$ $`{\displaystyle \underset{k}{}}𝒩_{i_1\mathrm{}i_pj_1\mathrm{}j_k^{}\mathrm{}j_p}^{(2p)}_{j_k^{}j_k}^\mathrm{T}(t)`$ (70)   
$``$ $`{\displaystyle \underset{k}{}}_{i_ki_k^{}}^{}(t)𝒩_{i_1\mathrm{}i_k^{}\mathrm{}i_pj_1\mathrm{}j_p}^{(2p)},`$ (71)   
$`(𝒟^{(2p2)}𝒩^{(2p2)})_{IJ}`$ $``$ $`{\displaystyle \underset{k,l}{}}𝒟_{i_kj_l}𝒩_{I_k^{}J_l^{}}^{(2p2)}`$ (72)   
with abbreviation of the indices $`Ii_1\mathrm{}i_p`$, $`Jj_1\mathrm{}j_p`$, $`I_k^{}i_1\mathrm{}i_{k1}i_{k+1}\mathrm{}i_p`$ and $`J_l^{}j_1\mathrm{}j_{l1}j_{l+1}\mathrm{}j_p`$.

We examine, in particular, the second-order moment $`𝒩_{ij}(t)𝒩_{ij}^{(2)}(t)`$. The number of the atoms which are excited by absorbing the photons is given by the diagonal component of $`𝒩(t)`$ as

$$n_b(t)=b^{}b=𝒩_{bb}(t).$$
(73)   
The master equation for $`𝒩(t)`$ reads

$$\frac{d𝒩}{dt}=i𝒩^\mathrm{T}(t)+i^{}(t)𝒩+𝒟.$$
(74)   
This linear equation with an inhomogeneous term may be solved formally as follows. We first introduce a new matrix $`𝒩^{}(t)`$ instead of $`𝒩(t)`$ by

$$𝒩^{}(t)=𝒰^1(t)𝒩(t)𝒰^{1\mathrm{T}}(t).$$
(75)   
The linear transformation $`𝒰(t)`$ representing the time evolution due to $`(t)`$ is given by

$$𝒰(t)=\mathrm{P}\left[i\mathrm{exp}_0^t(\tau )𝑑\tau \right]$$
(76)   
satisfying the equation $`d𝒰/dt=i(t)𝒰`$ with $`𝒰(0)=\mathrm{𝟏}`$, where “P” denotes the chronological product. Then, the master equation (74) is reduced as

$$\frac{d𝒩^{}}{dt}=𝒰^1(t)𝒟𝒰^{1\mathrm{T}}(t).$$
(77)   
By considering the initial condition $`𝒩^{}(0)=𝒩(0)`$ with $`𝒰(0)=\mathrm{𝟏}`$ and Eq. (63) for $`𝒟`$, we can solve the above equation as

$$𝒩_{ij}^{}(t)=_{ij}^k(t)\overline{n}_k+𝒩_{ij}(0),$$
(78)   
where 

$$_{ij}^k(t)=_0^t\gamma _k𝒰_{ik}^1(\tau )𝒰_{kj}^{1\mathrm{T}}(\tau )𝑑\tau .$$
(79)   
Then, we obtain the solution as

$$𝒩_{ij}(t)=_{ij}^k(t)\overline{n}_k+𝒰_{ik}^{}(t)𝒩_{kl}(0)𝒰_{lj}^\mathrm{T}(t),$$
(80)   
where

$$_{ij}^k(t)=\left[𝒰^{}(t)^k(t)𝒰^\mathrm{T}(t)\right]_{ij}.$$
(81)   
The same result is also obtained in the Langevin picture (see Appendix B). This solution approaches to certain asymptotic value after a long enough time as

$$𝒩_{ij}(t)_{ij}^k(t)\overline{n}_k(t\mathrm{max}[\gamma _i^1]).$$
(82)   
In practical calculations, the master equation (74) will be solved numerically with the time-dependent effective Hamiltonian $`(t)`$ as given in Eq. (64).

## VI Aspects of axion-photon-atom interaction

We can see some characteristic properties of the axion-photon-atom interaction in the resonant cavity by examining the simple case with the constant atom-photon coupling $`\mathrm{\Omega }(t)=\mathrm{\Omega }_N`$ for the effective Hamiltonian $`(t)=`$ in Eq. (64). In this case, as derived in the Appendix A, the analytic solution is obtained for the particle numbers (with the condition $`\overline{n}_b=0`$) as

$$n_i(t)=q_i^{}q_i=r_{ic}(t)\overline{n}_c+r_{ia}(t)\overline{n}_a,$$
(83)   
where

$$r_{ij}(t)=\underset{m,n}{}g_{ij}^mg_{ij}^n\left[\left(1\frac{\gamma _j}{\mathrm{\Lambda }_{mn}}\right)\mathrm{e}^{\mathrm{\Lambda }_{mn}t}+\frac{\gamma _j}{\mathrm{\Lambda }_{mn}}\right]$$
(84)   
with 

$`g_{ij}^k`$ $`=`$ $`\underset{si\lambda _k}{lim}(s+i\lambda _k)(s\mathrm{𝟏}+i)_{ij}^1,`$ (85)   
$`\mathrm{\Lambda }_{mn}`$ $`=`$ $`i(\lambda _m^{}\lambda _n).`$ (86)   
Here, the atomic damping rate $`\gamma _b`$ may be neglected for simplicity, since it is sufficiently smaller than $`\gamma _a`$ and $`\mathrm{\Omega }_N`$ ($`\gamma _b0.001\gamma `$ with $`\tau _b10^3\mathrm{s}`$ for $`\omega _c10^5\mathrm{eV}`$ and $`Q10^4`$). The condition $`\omega _b=\omega _c`$ may also be taken for definiteness, since the atomic transition frequency should be tuned almost equal to the cavity frequency. Then, the eigenvalues of the Hamiltonian $``$ are given by

$`\lambda _1`$ $`=`$ $`\omega _c{\displaystyle \frac{i}{4}}\gamma +i{\displaystyle \frac{(\gamma ^216\mathrm{\Omega }_N^2)^{1/2}}{4}},`$ (87)   
$`\lambda _2`$ $`=`$ $`\omega _c{\displaystyle \frac{i}{4}}\gamma i{\displaystyle \frac{(\gamma ^216\mathrm{\Omega }_N^2)^{1/2}}{4}},`$ (88)   
$`\lambda _3`$ $`=`$ $`\omega _a{\displaystyle \frac{i}{2}}\gamma _a,`$ (89)   
where

$$(\gamma ^216\mathrm{\Omega }_N^2)^{1/2}=\{\begin{array}{cc}\sqrt{\gamma ^216\mathrm{\Omega }_N^2}\hfill & (\mathrm{\Omega }_N/\gamma 1/4)\hfill \\ i\sqrt{16\mathrm{\Omega }_N^2\gamma ^2}\hfill & (\mathrm{\Omega }_N/\gamma >1/4)\hfill \end{array}.$$
(90)   
If the number of atoms $`N`$ (or the atomic beam intensity $`I_{\mathrm{Ryd}}`$) is not so large giving $`\mathrm{\Omega }_N/\gamma <1/4`$, the damping rate of the eigenmode of $`\lambda _1`$ is smaller than $`\gamma /2`$, and that of $`\lambda _2`$ lies between $`\gamma /2`$ and $`\gamma `$. On the other hand, in the strong coupling region of $`\mathrm{\Omega }_N/\gamma >1/4`$ for the collective atom-photon interaction the eignemodes of $`\lambda _1`$ and $`\lambda _2`$ form a doublet around the frequency $`\omega _c`$ with the same damping rate $`\gamma /2`$ (Rabi splitting).

Among the damping rates $`\mathrm{Re}[\mathrm{\Lambda }_{mn}]`$ of the respective terms in the factors $`r_{ij}(t)`$ representing the contributions of thermal photons and axions, $`\mathrm{Re}[\mathrm{\Lambda }_{33}]=\gamma _a`$ ($`0.01\gamma `$ for $`\beta _a10^3`$ and $`Q10^4`$, typically) is the smallest one. The rate $`\mathrm{Re}[\mathrm{\Lambda }_{11}]4\gamma (\mathrm{\Omega }_N/\gamma )^2`$ for $`\mathrm{\Omega }_N/\gamma 0.1`$ may also be comparable to the smallest rate $`\mathrm{Re}[\mathrm{\Lambda }_{33}]`$. The atomic transit time through the cavity is, on the other hand, given by

$$t_{\mathrm{tr}}=L/v$$
(91)   
with the detection cavity length $`L`$ and the atomic velocity $`v`$. This transit time provides the effective cut-off for the axion-photon-atom interaction in the cavity. (Here we assume for simplicity that the atoms have the uniform velocity.) It is typically $`t_{\mathrm{tr}}400\tau _\gamma `$ with $`L=0.2\mathrm{m}`$ and $`v=350\mathrm{m}\mathrm{s}^1`$ for $`m_a=10^5\mathrm{eV}`$ and $`Q=2\times 10^4`$ in the case of the detection apparatus such as CARRACK I. Hence, the transit time can be regareded to be long enough compared to $`(\mathrm{Re}[\mathrm{\Lambda }_{mn}])^1`$, i.e.,

$$t_{\mathrm{tr}}>\mathrm{several}\gamma _a^1,$$
(92)   
so that the respective particle numbers will almost reach the asymptotic values as

$$r_{ij}(t_{\mathrm{tr}})r_{ij}(\mathrm{})=\underset{m,n}{}g_{ij}^mg_{ij}^n\frac{\gamma _j}{\mathrm{\Lambda }_{mn}}.$$
(93)   
By using Eqs. (87) – (89) and the explicit matrix form (A21) for $`(s\mathrm{𝟏}+i)^1`$ given in Appendix A, we can calculate the coefficients $`g_{ij}^k`$ in Eq. (85). Then, we can show the relations

$$r_{bc}(\mathrm{})=r_{cc}(\mathrm{})=1.$$
(94)   
This implies that if the axion-photon interaction is turned off, the numbers of the photons and excited atoms reach the same asymptotic value $`\overline{n}_c`$ of the thermal photon number at $`T_c`$ :

$$n_b[cb]n_c[cc]\overline{n}_c.$$
(95)   
The number of axion-converted photons is, on the other hand, given approximately by

$$n_c[ac]r_{ca}(\mathrm{})\overline{n}_a.$$
(96)   
The number of excited atoms due to the axion-converted photons is also given by

$$n_b[acb]r_{ba}(\mathrm{})\overline{n}_a.$$
(97)   
These contributions from the axions are essentially determined by the factors

$$g_{ca}^m,g_{ba}^m\frac{\kappa }{(\lambda _m\lambda _k)(\lambda _m\lambda _l)},$$
(98)   
where $`mk,l`$ and $`kl`$. The detuning of axion mass from the cavity frequency is given by

$$\mathrm{\Delta }\omega _a\omega _a\omega _c.$$
(99)   
Then, the above factors are enhanced if the axion detuning $`\mathrm{\Delta }\omega _a`$ lies in the range where the resonant conditions $`\lambda _3\lambda _1`$ and/or $`\lambda _3\lambda _2`$ are satisfied.

In the leading order of the axion-photon coupling $`\kappa `$ with Eq. (98), the factors for the axion contributions are given with $`r_{ca}(\mathrm{}),r_{ba}(\mathrm{})(\kappa /\gamma )^2`$. Then, by noting the relation $`(\kappa /\gamma )^2\overline{n}_a(\rho _a/m_a)V_1`$ from Eqs. (19) and (41) it is found that the axion-converted photons $`n_c[ac]`$ and the number of atoms $`n_b[acb]`$ excited by such photons are both proportional to the number of axions contained in the conversion cavity. It is here relevant to define the form factors for the axion contributions with respect to the axion detuning ($`t_{\mathrm{tr}}>\mathrm{several}\gamma _a^1\gamma ^1`$) by

$$\sigma _{ia}(\mathrm{\Delta }\omega _a)=\frac{r_{ia}(t_{\mathrm{tr}})}{[\kappa /(\gamma /2)]^2}\frac{r_{ia}(\mathrm{})}{[\kappa /(\gamma /2)]^2}(i=b,c).$$
(100)   
These form factors $`\sigma _{ba}(\mathrm{\Delta }\omega _a)`$ for $`acb`$ (solid lines) and $`\sigma _{ca}(\mathrm{\Delta }\omega _a)`$ for $`ac`$ (dotted lines) are plotted together in Fig. 2 for some typical values of the atom-photon coupling, $`\mathrm{\Omega }_N/\gamma =0.1,0.5,1.0`$.

The behavior of $`\sigma _{ba}(\mathrm{\Delta }\omega _a)`$ is, in particular, understood by noting its dominant contribution which is in fact given by the term of $`|g_{ba}^{m=3}|^2`$ ($`\gamma _a\gamma `$) as

$$\sigma _{ba}(\mathrm{\Delta }\omega _a)\frac{4\mathrm{\Omega }_N^2\gamma ^2}{4\mathrm{\Delta }\omega _a^2\gamma (\gamma 4\gamma _a)+(4\mathrm{\Delta }\omega _a^24\mathrm{\Omega }_N^2+\gamma _a\gamma )^2}.$$
(101)   
The peak of this form factor appears at

$$\mathrm{\Delta }\omega _a(\mathrm{peak})=\{\begin{array}{cc}0\hfill & (\mathrm{\Omega }_N\overline{\mathrm{\Omega }})\hfill \\ \pm \sqrt{\mathrm{\Omega }_N^2\overline{\mathrm{\Omega }}^2}\hfill & (\mathrm{\Omega }_N>\overline{\mathrm{\Omega }})\hfill \end{array},$$
(102)   
where $`\overline{\mathrm{\Omega }}=\gamma /\sqrt{8}+O(\gamma _a)`$. (Although $`\sigma _{ba}(0)`$ is apparently divergent for $`\mathrm{\Omega }_N=\frac{1}{2}(\gamma _a\gamma )^{1/2}`$ in the approximate formula of Eq. (101), the other contributions also become significant around $`\mathrm{\Delta }\omega _a=0`$ so as to give a finite value of $`\sigma _{ba}(0)`$.) It is found that for $`\mathrm{\Omega }_N\gamma /4`$ the peak value of the axion signal decreases due to the Rabi splitting, which approaches to $`\sigma _{ba}(\mathrm{\Delta }\omega _a(\mathrm{peak}))1`$ for $`(\mathrm{\Omega }_N/\gamma )^21`$. On the other hand, for the appropriate atom-photon coupling such as

$$\mathrm{\Omega }_N(\gamma _a\gamma )^{1/2}$$
(103)   
the signal is enhanced significantly, as observed in Fig. 2 ($`\mathrm{\Omega }_N=0.1\gamma `$ and $`\gamma _a=0.02\gamma `$), by virtue of the narrow width (long coherence) of the galactic axions as

$$\sigma _{ba}(\mathrm{\Delta }\omega _a)(\gamma _a/\gamma )^1,$$
(104)   
when the axion detuning becomes small enough to satisfy the condition

$$|\mathrm{\Delta }\omega _a|\gamma _a.$$
(105)   
Here, the energy uncertainty $`\mathrm{}/t_{\mathrm{tr}}`$ due to the atomic motion is assumed to be smaller than the axion width $`\mathrm{}\gamma _a`$.

The form factor $`\sigma _{ca}(\mathrm{\Delta }\omega _a)`$ determines the equilibrium number $`r_{ca}(\mathrm{})\overline{n}_a=4\sigma _{ca}(\mathrm{\Delta }\omega _a)(\kappa /\gamma )^2\overline{n}_a`$ of the axion-converted photons in the cavity. As seen in Fig. 2, its peak value is $`\sigma _{ca}1`$ almost independently of the atom-photon coupling $`\mathrm{\Omega }_N`$. It is here, in particular, interesting to observe a narrow dip in $`\sigma _{ca}(\mathrm{\Delta }\omega _a)`$ around $`\mathrm{\Delta }\omega _a=0`$ for $`\mathrm{\Omega }_N/\gamma =0.1`$. This indicates that the axion-converted photons are efficiently absorbed by the atoms for $`\mathrm{\Omega }_N/\gamma (\gamma _a/\gamma )^{1/2}`$. For larger $`\mathrm{\Omega }_N`$, two separate peaks appear in $`\sigma _{ca}(\mathrm{\Delta }\omega _a)`$ due to the Rabi splitting.

Some characteristic features concerning the axion-photon-atom interaction in the resonant cavity have been discussed so far by making the analytic calculations for the simple case with the constant atom-photon coupling $`\mathrm{\Omega }_N`$. They will indeed be confirmed in Sec. IX by performing detailed numerical calculations for the realistic situation with the continuous atomic beam passing through the spatially varying electric field.

## VII Detection sensitivity with continuous atomic beam

In order to make precise estimates for the sensitivity of the Rydberg atom cavity detector, we have to take into accout (i) the motion and (ii) the almost uniform distribution of the atoms in the incident beam as well as (iii) the spatial variation of the electric field in the cavity. We here elaborate the calculations presented in the preceeding sections by treating these points appropriately.

The electric field felt by the atoms varies with time through the atomic motion. Accordingly, the effect of atomic motion can be incorporated by introducing the relevant time-dependence for the atom-photon coupling, which is determined by the profile of the electric field in the cavity, as given in Eq. (66). On the other hand, in order to treat the spatial distribution of the atoms, we divide the atoms in the cavity into $`K`$ bunches with a fixed beam intensity

$$I_{\mathrm{Ryd}}=N/t_{\mathrm{tr}}.$$
(106)   
Here, $`N`$ is the total number of Rydberg atoms in the cavity. Then, the coutinuous atomic beam will actually be realized for $`K1`$.

The collective atomic mode of each bunch is denoted by $`b_i`$ ($`i=1,2,\mathrm{},K`$), and the effective Hamiltonian is given by a $`(K+2)\times (K+2)`$ matrix

$$(t)=\left(\begin{array}{cc}(\omega _b\frac{i}{2}\gamma _b)\mathrm{𝟏}& \begin{array}{cc}\mathrm{\Omega }_1(t)\hfill & 0\\ \mathrm{}\hfill & \mathrm{}\\ \mathrm{\Omega }_K(t)\hfill & 0\end{array}\\ \begin{array}{ccc}\mathrm{\Omega }_1& \mathrm{}& \mathrm{\Omega }_K(t)\\ 0& \mathrm{}& 0\end{array}& \begin{array}{cc}\omega _c\frac{i}{2}\gamma _c& \kappa \\ \kappa & \omega _a\frac{i}{2}\gamma _a\end{array}\end{array}\right).$$
(107)   
Now suppose that the $`i`$-th atomic bunch locates around

$$x_i=(i1)\delta x$$
(108)   
for a time interval

$$M\delta tt<(M+1)\delta t,$$
(109)   
where $`M=0,1,2,\mathrm{}`$, and

$$\delta x=L/K,\delta t=t_{\mathrm{tr}}/K.$$
(110)   
Then, the collective atom-photon coupling for the $`i`$-th bunch containing $`N/K`$ atoms is given with the electric field profile by

$$\mathrm{\Omega }_i(t)=(\mathrm{\Omega }_N/\sqrt{K})f(x_i+v\overline{t}),$$
(111)   
where

$$\overline{t}tM\delta t(0\overline{t}<\delta t).$$
(112)   
After each time interval of $`\delta t`$, the $`K`$-th atomic bunch leaves the cavity, and a new one comes in. Then, the collective atomic modes of the respective bunches should be replaced as

$$b_ib_{i+1},$$
(113)   
which is represented by a $`(K+2)\times (K+2)`$ matrix

$$𝒫_{}=\left(\begin{array}{cccccccc}0& 0& 0& \mathrm{}& 0& 0& 0& 0\\ 1& 0& 0& \mathrm{}& 0& 0& 0& 0\\ 0& 1& 0& \mathrm{}& 0& 0& 0& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& 0& \mathrm{}& 0& 0& 0& 0\\ 0& 0& 0& \mathrm{}& 1& 0& 0& 0\\ 0& 0& 0& \mathrm{}& 0& 0& 1& 0\\ 0& 0& 0& \mathrm{}& 0& 0& 0& 1\end{array}\right).$$
(114)   
In accordance with this replacement of the collective atomic modes at each period of $`\delta t`$, the time evolution of the particle number matrix $`𝒩_{ij}(t)`$ should be determined by dividing it into the corresponding parts in time as

$$\overline{𝒩}(\overline{t},M)𝒩(t)(M\delta tt<(M+1)\delta t).$$
(115)   
These parts are connected together at $`t=M\delta t`$ with a suitable matching condition

$$\overline{𝒩}(0,M+1)=𝒫_{}\overline{𝒩}(\delta t,M)𝒫_{}^\mathrm{T}.$$
(116)   
This condition is written down explicitly as

$$\begin{array}{cc}\overline{𝒩}_{ij}(0,M+1)=\overline{𝒩}_{i1j1}(\delta t,M)\hfill & [2i,jK],\hfill \\ \overline{𝒩}_{1j}(0,M+1)=0\hfill & [1jK+2],\hfill \\ \overline{𝒩}_{i1}(0,M+1)=0\hfill & [1iK+2],\hfill \\ \overline{𝒩}_{ij}(0,M+1)=\overline{𝒩}_{i1j}(\delta t,M)\hfill & [2iK,j=a,c],\hfill \\ \overline{𝒩}_{ij}(0,M+1)=\overline{𝒩}_{ij1}(\delta t,M)\hfill & [i=a,c,2jK],\hfill \\ \overline{𝒩}_{ij}(0,M+1)=\overline{𝒩}_{ij}(\delta t,M)\hfill & [i,j=a,c].\hfill \end{array}$$
(117)   
Here, the second and third lines imply that the incoming atomic mode ($`b_1`$) does not have any correlation initially with the other quantum modes ($`b_2,\mathrm{},b_K,c,a`$) already interacting in the cavity.

In this way, by solving the master equation (74) for many time intervals, we obtain the steady solution

$$\overline{𝒩}(\overline{t})\overline{𝒩}(\overline{t},M1).$$
(118)   
It in fact appears independently of the choice of the initial value $`𝒩(0)`$, which is due to the dissipation of the quantum modes for the axions, photons and atoms. The number of excited atoms contained in each atomic bunch is then given by

$$n_{b_i}(\overline{t})=\overline{𝒩}_{ii}(\overline{t})=n_{b_i}^a(\overline{t})+n_{b_i}^c(\overline{t}).$$
(119)   
Here, the contributions of the axions and thermal photons are proporional to the respective particle numbers as

$$n_{b_i}^a(\overline{t})=r_{b_ia}(\overline{t})\overline{n}_a,n_{b_i}^c(\overline{t})=r_{b_ic}(\overline{t})\overline{n}_c.$$
(120)   
The distribution of excited atoms in the cavity is determined with this steady solution as

$$\overline{\rho }_b(x)=\frac{n_{b_i}(\overline{t}=x/v)}{\delta x}((i1)\delta xx<i\delta x).$$
(121)   
The $`K`$-th atomic bunch exits the cavity at $`\overline{t}=\delta t`$, and the excited atoms contained there are detected. Accordingly, the counting rates for the contributions of the axions and thermal photons are calculated, respectively, for large enough $`K`$ and $`M`$ by

$$R_s=\frac{n_{b_K}^a(\delta t)}{\delta t},$$
(122)   
$$R_n=\frac{n_{b_K}^c(\delta t)}{\delta t}.$$
(123)   
By using these counting rates, the measurement time required to search for the axion signal at the confidence level of $`m\sigma `$ is estimated as

$$\mathrm{\Delta }t=\frac{m^2(1+R_n/R_s)}{R_s}.$$
(124)   
In the search for the axions with unknown mass, the cavity frequency ($`\omega _b=\omega _c`$ for definiteness) should be changed with appropriate scanning step $`\mathrm{\Delta }\omega _c`$. The total scanning time over a 10% frequency range is then given by

$$t_{\mathrm{tot}}=\frac{0.1\omega _c}{\mathrm{\Delta }\omega _c}\mathrm{\Delta }t.$$
(125)   
The frequency step $`\mathrm{\Delta }\omega _c`$ is determined by considering the resonant condition for the absorption of the axion-converted photons by the Rydberg atoms. Specifically, as seen in Eqs. (103), (104) and (105), the signal rate is enhanced significantly for the axion detuning $`|\mathrm{\Delta }\omega _a|\gamma _a`$ with $`\mathrm{\Omega }_N/\gamma (\gamma _a/\gamma )^{1/2}`$. Hence, the scanning frequency step should be taken as

$$\mathrm{\Delta }\omega _c\gamma _a.$$
(126)   
Then, if the axion really exists in a mass range searched with this frequency step, the resonant condition for the axion signal can be satisfied at a certain scanning step.

## VIII Dependence on the relevant parameters

We here consider how the counting rates of the axion signal and thermal photon noise depend on the relevant experimental parameters, before presenting detailed numerical calculations in the next section. Although the notation for the case of $`K=1`$ with constant atom-photon coupling $`\mathrm{\Omega }_N`$ is used for simplicity, the essential features argued below are indeed valid for the realistic case with continuous atomic beam.

The energy scales involved in the effective Hamiltonian are as follows:

$`\begin{array}{ccc}\mathrm{𝐝𝐞𝐭𝐮𝐧𝐢𝐧𝐠𝐬}\hfill & :& \mathrm{\Delta }\omega _a\omega _a\omega _c,\mathrm{\Delta }\omega _b\omega _b\omega _c,\hfill \\ \mathrm{𝐜𝐨𝐮𝐩𝐥𝐢𝐧𝐠𝐬}\hfill & :& \kappa ,\mathrm{\Omega }_N,\hfill \\ \mathrm{𝐝𝐚𝐦𝐩𝐢𝐧𝐠𝐬}\hfill & :& \gamma _a,\gamma _b,\gamma _c\gamma .\hfill \end{array}`$

The atomic transit time is also important to determine the time evolution of the system, which is controlled by the atomic velocity $`v`$ with a fixed cavity length $`L`$ as

$$v=L/t_{\mathrm{tr}}.$$
(127)   
The effects of the dissipation of axions and atoms are rather small compared to that of the photons in the actual situation with $`\gamma _a,\gamma _b\gamma `$. The atomic transit time is in fact much longer than the other characteristic time scales except for $`\kappa ^1`$, e.g., $`t_{\mathrm{tr}}400\tau _\gamma `$ for $`v=350\mathrm{m}\mathrm{s}^1`$ with $`L=0.2\mathrm{m}`$, $`Q=2\times 10^4`$ and $`m_a=10^5\mathrm{eV}`$. In this situation, the quantum averaged particle numbers approach to the asymptotic values with the factors $`r_{ij}(t_{\mathrm{tr}})r_{ij}(\mathrm{})`$, as given in Eq. (93). Then, the counting rates of the signal and noise, which are calculated by $`R_s=r_{ba}(t_{\mathrm{tr}})\overline{n}_a/t_{\mathrm{tr}}`$ and $`R_n=r_{bc}(t_{\mathrm{tr}})\overline{n}_c/t_{\mathrm{tr}}`$, are found to be roughly proportional to the atomic velocity $`v`$.

In the weak atom-photon coupling region of $`\mathrm{\Omega }_N/\gamma 0.1`$, the counting rates $`R_s`$ and $`R_n`$ are roughly proportional to $`\mathrm{\Omega }_N^2NI_{\mathrm{Ryd}}`$, as seen in Eq. (101) for the axion signal. On the other hand, when the atom-photon coupling becomes significant with sufficient atomic beam intensity $`I_{\mathrm{Ryd}}`$, the noise rate saturates to certain value $`R_n\overline{n}_c/t_{\mathrm{tr}}`$, and the signal rate $`R_s`$ is maximized, as seen in Eqs. (103) and (104), for the atom-photon coupling

$$\mathrm{\Omega }_N/\gamma (\gamma _a/\gamma )^{1/2}0.1\left(\frac{\beta _a}{10^3}\right)\left(\frac{Q}{10^4}\right)^{1/2}.$$
(128)   
It is possible to attain this optimal coupling in the actual detection apparatus by preparing the appropriate number of Rydberg atoms with a suitable laser system. The beam intensity of Rydberg atoms $`I_{\mathrm{Ryd}}=N/t_{\mathrm{tr}}`$ providing this optimal value (128) for the collective atom-photon coupling $`\mathrm{\Omega }_N=\mathrm{\Omega }\sqrt{N}`$ is chosen depending on the relevant parameters as

$$I_{\mathrm{Ryd}}v\beta _a^2m_a^2Q^1\mathrm{\Omega }^2.$$
(129)   
It is typically estimated as $`I_{\mathrm{Ryd}}7\times 10^5\mathrm{s}^1`$ with $`N400`$ atoms for $`m_a=10^5\mathrm{eV}`$, $`\beta _a=10^3`$, $`Q=2\times 10^4`$, $`\mathrm{\Omega }=5\times 10^3\mathrm{s}^1`$, $`v=350\mathrm{m}\mathrm{s}^1`$ and $`L=0.2\mathrm{m}`$ ($`t_{\mathrm{tr}}400\tau _\gamma `$). On the other hand, if the atom-photon coupling becomes too strong, the signal rate decreases due to the appearance of Rabi splitting, as seen in Fig. 2.

We have observed that the atom-photon interaction can not be treated perturbatively for $`\mathrm{\Omega }_N/\gamma 0.1`$ . The conversion between atoms and photons appears to be reversible in this case. On the other hand, since the axion-photon coupling is extremely small, the conversion of axions to photons can well be treated perturbatively as an irreversible process. The signal rate $`R_s`$ is then found to be proportional to $`(\kappa /\gamma )^2\overline{n}_a`$ in a very good approximation. Then, by considering the relations (19), (29) and (41) with $`\omega _cm_a/\mathrm{}`$, the dependence of the signal rate on the relevant parameters is specified as

$$R_svc_{a\gamma \gamma }^2B_{\mathrm{eff}}^2Q^2V_1(\rho _a/m_a)\overline{\sigma }_{ba}.$$
(130)   
Here, the form factor $`\overline{\sigma }_{ba}=\sigma _{ba}(\pm \mathrm{\Delta }\omega _c/2)`$ for $`\mathrm{\Delta }\omega _a=\pm \mathrm{\Delta }\omega _c/2`$ ($`\mathrm{\Delta }\omega _b=0`$) is also indicated. This form factor actually depends on the choices of the atom-photon coupling $`\mathrm{\Omega }_N`$ and the scanning frequency step $`\mathrm{\Delta }\omega _c`$, as seen in Eqs. (103), (104) and (105). It is noticed in Eq. (130) that the signal rate is proportional to the number of axions contained in the conversion cavity $`(\rho _a/m_a)V_1`$. The dependence of the signal rate on the axion mass $`m_a`$ is further specified by considering the relation for the conversion cavity volume

$$V_1m_a^2.$$
(131)   
This relation is due to the fact that the diameter of the cavity is proportional to $`m_a^1`$ while its length is fixed.

To summarize these arguments, the optimal situation, as given in Eqs. (103), (104) and (105), can be set up experimentally for the dark matter axion search with Rydberg atom cavity detector. Then, the counting rates are expected to behave with respect to the changes of the relevant parameters as

$`R_s`$ $``$ $`vc_{a\gamma \gamma }^2B_{\mathrm{eff}}^2Q\beta _a^2m_a^3\rho _a,`$ (132)   
$`R_n`$ $``$ $`v\overline{n}_c.`$ (133)   
Here, the noise rate is proportional to the thermal photon number $`\overline{n}_c[m_a/T_c]`$ which is determined as a function of the ratio $`m_a/T_c`$. The resonant value $`\overline{\sigma }_{ba}(\gamma _a/\gamma )^1\beta _a^2Q^1`$ is obtained for the axion-photon-atom interaction with the small enough axion detuning $`|\mathrm{\Delta }\omega _a|\mathrm{\Delta }\omega _c/2\gamma _a`$. (This is valid as long as the quantum uncertainty of energy $`\mathrm{}/t_{\mathrm{tr}}`$ due to the finite atomic transit time is smaller than the axion width $`\mathrm{}\gamma _a`$, as ensured in the actual detection system.) It should also be remarked that the small detuning $`\mathrm{\Delta }\omega _b`$ of the atomic transition frequency only shifts slightly the location of the peak of the axion signal from $`\mathrm{\Delta }\omega _a(\mathrm{peak})=0`$ to $`\mathrm{\Delta }\omega _a(\mathrm{peak})\mathrm{\Delta }\omega _b`$. By taking these signal and noise rates in the optimal case, the measurement time $`\mathrm{\Delta }t`$ at each scanning step is estimated with respect to the relevant parameters as

$$\mathrm{\Delta }t\{\begin{array}{cc}v^1c_{a\gamma \gamma }^4B_{\mathrm{eff}}^4Q^2\beta _a^4m_a^6\rho _a^2\overline{n}_c\hfill & (R_s/R_n<1)\hfill \\ c_{a\gamma \gamma }^2B_{\mathrm{eff}}^2Q^2\beta _a^2m_a^3\rho _a^1\hfill & (R_s/R_n1)\hfill \end{array}.$$
(134)   
The total scanning time is then estimated with the appropriate scanning step $`\mathrm{\Delta }\omega _c\gamma _a`$ as

$$t_{\mathrm{tot}}\mathrm{\Delta }t/\beta _a^2.$$
(135)   
Here the negative power of $`\beta _a`$ appearing in the right-hand side is indeed compensated by the positive power of $`\beta _a`$ contained in Eq. (134) for $`\mathrm{\Delta }t`$. The sensitivity seems to be improved apparently for the high atomic velocity in Eqs. (134) and (135). It should, however, be remarked that the condition $`\mathrm{}/t_{\mathrm{tr}}<\mathrm{}\gamma _a`$ for the resonant axion-photon-atom interaction is no longer ensured if the atomic velocity beccomes too high. The high atomic beam intensity is even required in such a case in accordance with Eq. (129). Hence, we will in fact see in the next section that the preferable atomic velocity is $`v100\mathrm{m}\mathrm{s}^11000\mathrm{m}\mathrm{s}^1`$ for the actual experimental apparatus.

Some essential features have been discussed so far concerning the detection sensitivity of the Rydberg atom cavity detector. In the next section, they will indeed be confirmed by detailed numerical calculations for the realistic situation with continuous atomic beam.

## IX Numerical analysis

Numerical calculations have been performed to determine precisely the quantum evolution of the axion-photon-atom system in the resonant cavity, where some practical values are taken for the experimental parameters such as

$`\begin{array}{c}m_a=3\times 10^6\mathrm{eV}3\times 10^5\mathrm{eV},\hfill \\ \beta _a=3\times 10^43\times 10^3,\hfill \\ \rho _a=\rho _{\mathrm{halo}}=0.3\mathrm{GeVcm}^3,\hfill \\ T_c=10\mathrm{m}\mathrm{K}15\mathrm{m}\mathrm{K},\hfill \\ Q=1\times 10^47\times 10^4,B_{\mathrm{eff}}=4\mathrm{T},\hfill \\ V_1=5000\mathrm{c}\mathrm{m}^3\left(m_a/10^5\mathrm{eV}\right)^2,L=0.2\mathrm{m},\hfill \\ I_{\mathrm{Ryd}}=10^3\mathrm{s}^110^7\mathrm{s}^1,\hfill \\ v=350\mathrm{m}\mathrm{s}^110000\mathrm{m}\mathrm{s}^1,\hfill \\ \mathrm{\Omega }=5\times 10^3\mathrm{s}^1,\tau _b=10^3\mathrm{s}.\hfill \end{array}`$

The steady state is realized by solving the master equation (74) for a long enough time interval $`M\delta t`$ ($`M1`$). It in practice appears independently of the choice of initial value $`𝒩(0)`$ due to the finite damping rates of the axions, photons and atoms. The spatial distribution of the excited atoms in the cavity is calculated by Eq. (121) with this steady solution. The contributions of the galactic axions $`\overline{\rho }_b^{[a]}(x)`$ and the thermal photons $`\overline{\rho }_b^{[\gamma ]}(x)`$ are depicted together in Figs. 3 and 4 for the case of DFSZ axion with the tanh-type and sine-type electric field profiles, respectively, where the relevant parameters are taken as $`m_a=10^5\mathrm{eV}`$ ($`\omega _c=\omega _a`$), $`T_c=12\mathrm{m}\mathrm{K}`$, $`Q=2\times 10^4`$, $`\mathrm{\Omega }_N/\gamma =0.1`$ and $`v=350\mathrm{m}\mathrm{s}^1`$. These results obtained with $`K=10`$ and $`M=20`$ in fact appear to be smooth enough indicating that the steady state in the case with continuous atomic beam is realized in a good approximation. It is here remarkable that the form of electric field profile is reflected in the distribution of the excited atoms. In the following calculations, we take the tanh-type electric field profile for definiteness. The essential results concerning the interaction between axions and atoms mediated by photons in the resonant cavity, however, hold even in the cases with more realistic electric profiles.

We can compute the signal and noise rates with the steady solutions obtained in this way. (The following calculations are made for $`K=5`$ and $`M=10`$ with the tanh-type electric field profile. The results are changed less than 10% even if larger $`K`$ and $`M`$ are taken.) It is readily checked by these calculations that the signal rate $`R_s`$ is indeed proportional to $`c_{a\gamma \gamma }^2B_{\mathrm{eff}}^2V_1(\rho _a/m_a)`$, as shown in Eq. (130). Hence, the signal rate for the KSVZ axion is larger than that for the DFSZ axion by a factor $`c_{a\gamma \gamma }^2(\mathrm{KSVZ})/c_{a\gamma \gamma }^2(\mathrm{DFSZ})7`$.

In Fig. 5, we plot the signal and noise rates, $`R_s`$ and $`R_n`$, together depending on the atomic beam intensity $`I_{\mathrm{Ryd}}`$ (and the corresponding atom-photon coupling $`\mathrm{\Omega }_N`$), where the relevant parameters are taken as $`m_a=10^5\mathrm{eV}`$ ($`\omega _c=\omega _a`$), $`Q=2\times 10^4`$, $`L=0.2\mathrm{m}`$, $`v=350\mathrm{m}\mathrm{s}^1`$, $`\mathrm{\Omega }=5\times 10^3\mathrm{s}^1`$ and $`T_c=10\mathrm{m}\mathrm{K},12\mathrm{mK},15\mathrm{mK}`$. The noise rate is clearly proportional to the thermal photon number $`\overline{n}_c[m_a/T_c]`$, which is $`\overline{n}_c=8.9\times 10^6,6.2\times 10^5,4.3\times 10^4`$ at $`T_c=10\mathrm{m}\mathrm{K}`$, $`12\mathrm{m}\mathrm{K}`$, $`15\mathrm{m}\mathrm{K}`$, respectively, for the axion mass $`m_a=10^5\mathrm{eV}`$. In the weak beam intensity region, the signal and noise rates increase almost proportional to $`I_{\mathrm{Ryd}}=N/t_{\mathrm{tr}}`$. On the other hand, for sufficiently strong beam intensities the noise rate is saturated to certain asymptotic value. In this case, the atoms interact reversibly with the photons in the cavity so that the equilibrium with $`n_b[cb]\overline{n}_c`$ is realized. (In practice, due to the finite atomic lifetime the total number of excited atoms $`n_b[cb]=_0^L\overline{\rho }_b^{[\gamma ]}(x)𝑑x`$ obtained from the thermal photons is slightly different from $`\overline{n}_c`$.) As for the signal rate, we clearly observe in Fig. 5 that it is optimized for certain atomic beam intensity corresponding to the condition (128) for the atom-photon coupling. If the atomic beam intensity is too strong, the signal rate even decreases due to the Rabi splitting. It is hence important to adjust the beam intensity so as to attain the optimal condition for the signal.

It is also checked that the lines representing the noise rate $`R_s`$ versus beam intensity $`I_{\mathrm{Ryd}}`$ shifts horizontally by changes of the intrinsic atom-photon coupling $`\mathrm{\Omega }`$. Although the theoretical estimate of the electric dipole transition matrix element $`d`$ is available  for calculating $`\mathrm{\Omega }`$ with Eq. (48), some ambiguities would be present in such a naive theoretical estimate of $`\mathrm{\Omega }`$. Hence, it may rather be necessary to determine $`\mathrm{\Omega }`$ experimentally by fitting the thermal noise data in the weak beam intensity region to the expected lines, as shown in Fig. 5 for $`\mathrm{\Omega }=5\times 10^3\mathrm{s}^1`$, which are calculated by varying $`\mathrm{\Omega }`$ around some plausible value suggested by the theoretical calculation of $`d`$ .

Another aspect should be pointed out concerning the counting rates of the excited atoms. As given in Eqs. (121) and (123), the noise rate is calculated by the formula $`R_n=v\overline{\rho }_b^{[\gamma ]}(L)`$ with the excited atom density around the exit of cavity ($`x=L`$). Then, it may be noticed in Fig. 5 that the saturated value of $`R_n`$ in the strong beam intensity region, e.g., $`0.23\mathrm{s}^1`$ for $`T_c=12\mathrm{m}\mathrm{K}`$, is in fact significantly larger than the naive estimate $`\overline{n}_c/t_{\mathrm{tr}}`$, e.g., $`0.11\mathrm{s}^1`$ for $`\overline{n}_c[10^5\mathrm{eV}/12\mathrm{m}\mathrm{K}]=6.2\times 10^5`$, $`v=350\mathrm{m}\mathrm{s}^1`$ and $`L=0.2\mathrm{m}`$. This enhancement of the counting rates, which is expected for the signal as well as the noise, is due to the fact that the densities of excited atoms $`\overline{\rho }_b^{[a]}(L)`$ and $`\overline{\rho }_b^{[\gamma ]}(L)`$ around the exit of cavity become higher than the average values, as seen in Figs. 3 and 4. The nonuniform density of the excited atoms in the cavity is indeed brought by the fact that the atoms in the continuous beam are passing through the spatially varying electric field with finite transit time.

The behavior of the signal rate $`R_s`$ with respect to the axion detunig $`\mathrm{\Delta }\omega _a`$ is shown in Fig. 6. Here some typical values are taken for the atom-photon coupling, $`\mathrm{\Omega }_N/\gamma =0.1,0.3,0.7`$, and the other relevant parameters are chosen as $`m_a=10^5\mathrm{eV}`$, $`\beta _a=10^3`$, $`Q=2\times 10^4`$, $`L=0.2\mathrm{m}`$ and $`v=350\mathrm{m}\mathrm{s}^1`$. This behavior almost agrees with that of the form factor in Fig. 2 which is obtained for the simple case. (Note that the log-scale is taken for $`R_s`$ in Fig. 6.) It is therefore confirmed by these numerical calculations that the resonant axion-photon-atom interaction takes place if the axion detuning is small enough as $`|\mathrm{\Delta }\omega _a|\gamma _a`$ with the optimal atom-photon coupling $`\mathrm{\Omega }_N/\gamma (\gamma _a/\gamma )^{1/2}`$ (as long as $`\mathrm{}/t_{\mathrm{tr}}<\mathrm{}\gamma _a`$). It is also checked that if there is a small detuning $`\mathrm{\Delta }\omega _b`$ of the atomic transition frequency, this resonant condition is slightly modified as $`|\mathrm{\Delta }\omega _a\mathrm{\Delta }\omega _b|=|\omega _a\omega _b|\gamma _a`$.

As clearly observed in Fig. 6, a salient feature which is found by these calculations for the case with continous atomic beam is the appearance of a ripple structure in $`R_s`$ versus $`\mathrm{\Delta }\omega _a`$ with the relatively strong atom-photon coupling, $`\mathrm{\Omega }_N/\gamma >0.1`$ for the present choice of parameters. It is realized that this structure is brought from the narrow axion width $`\gamma _a\beta _a^2m_a/\mathrm{}`$. We have checked that this fine structure indeed disappears if the axion spectrum spreads much wider with higher mean velocity such as $`\beta _a=3\times 10^3`$.

The signal rate $`R_s`$ versus axion detuning $`\mathrm{\Delta }\omega _a`$ is calculated for some typical axion velocities, $`\beta _a=3\times 10^4,1\times 10^3,3\times 10^3`$. The results are shown in Fig. 7, where $`m_a=10^5\mathrm{eV}`$, $`Q=2\times 10^4`$, $`\mathrm{\Omega }_N/\gamma =0.1`$, $`L=0.2\mathrm{m}`$ and $`v=350\mathrm{m}\mathrm{s}^1`$ are taken. It is here clearly observed that this signal form factor exhibits the structure determined by the energy spread of axions $`\mathrm{}\gamma _a\beta _a^2m_a`$. The galactic axion spectrum can actually have some narrow peaks, as pointed out in . Then, by taking a small enough scanning frequency step $`\mathrm{\Delta }\omega _c`$, such peaks may be observed in the present detection scheme as well as the heterodyne method . The energy resolution is now limited by the qauntum uncertainty $`\mathrm{}/t_{\mathrm{tr}}`$, which could be improved by utilizing the lower velocity atomic beam.

The dependence of the axion signal $`R_s`$ on the atomic velocity $`v`$ is also shown in Fig. 8 for $`m_a=10^5\mathrm{eV}`$, $`\beta _a=10^3`$, $`\mathrm{\Omega }_N/\gamma =0.1`$, $`L=0.2\mathrm{m}`$ and $`v=350\mathrm{m}\mathrm{s}^1,2000\mathrm{m}\mathrm{s}^1,10000\mathrm{m}\mathrm{s}^1`$. Here, we note that for high atomic velocities such as $`v=10000\mathrm{m}\mathrm{s}^1`$ this axion signal form factor has a width much larger than that of axions. It is in fact determined by the energy uncertainty $`\mathrm{}/t_{\mathrm{tr}}`$ (e.g., $`0.1\mathrm{}\gamma `$ for $`v=10000\mathrm{m}\mathrm{s}^1`$) due to the atomic motion.

The measurement time $`\mathrm{\Delta }t`$ required for $`3\sigma `$ level at each scanning step is shown in Fig. 9 depending on the atomic beam intensity $`I_{\mathrm{Ryd}}`$. Here the axion detuning is taken for definiteness as $`|\mathrm{\Delta }\omega _a|\mathrm{\Delta }\omega _c/2`$ with the scanning frequency step

$$\mathrm{\Delta }\omega _c=0.05\gamma =6\mathrm{k}\mathrm{H}\mathrm{z}\left(\frac{m_a}{10^5\mathrm{eV}}\right)\left(\frac{2\times 10^4}{Q}\right),$$
(136)   
which is of the order of axion width $`\gamma _a`$ for $`\beta _a10^3`$ meeting the resonant condition (105) for the axion signal. Then, the total scanning time $`t_{\mathrm{tot}}`$ is estimated, as shown in Fig. 10. These calculations are made for the DFSZ axion by taking some typical values of $`T_c`$, $`Q`$ and $`v`$, which are indicated in Figs. 9 and 10. The other relevant parameters are chosen as $`m_a=10^5\mathrm{eV}`$, $`\beta _a=10^3`$, $`\mathrm{\Omega }_N/\gamma =0.1`$, $`L=0.2\mathrm{m}`$ and $`v=350\mathrm{m}\mathrm{s}^1`$. Here, we can see that the beam intensity to optimize the detection sensitivity changes with respect to $`Q`$ and $`v`$ according to the relation (129), as discussed before. The optimal estimates of the measurement time and the total scanning time also show roughly the dependence on $`Q`$ and $`v`$ as indicated in (134) and (135). For the KSVZ axion, the sensitivity is much better by a factor $`7^250`$ ($`R_s/R_n<1`$) or $`7`$ ($`R_s/R_n1`$).

The detection sensitivity for $`3\sigma `$ is finally estimated over the axion mass range $`m_a=3\times 10^6\mathrm{eV}3\times 10^5\mathrm{eV}`$, which can be searched with the present type of Rydberg atom cavity detector. In these calculations, the optimal atomic beam intensity is taken as $`I_{\mathrm{Ryd}}=4\times 10^5\mathrm{s}^1(m_a/10^5\mathrm{eV})^2(Q/3\times 10^4)^1`$ in accordance with the relation (129), which is also read from Figs. 9 and 10. The volume of the conversion cavity is taken as $`V_1=5000\mathrm{c}\mathrm{m}^3(m_a/10^5\mathrm{eV})^2`$ by considering the relation (131). The other relevant parameters are chosen as $`\beta _a=10^3`$, $`L=0.2\mathrm{m}`$, $`v=350\mathrm{m}\mathrm{s}^1`$, $`\mathrm{\Delta }\omega _c=0.05\gamma `$ and $`T_c=10\mathrm{m}\mathrm{K},12\mathrm{mK},15\mathrm{mK}`$. The results are shown in Figs. 11 and 12. It should here be noted that the $`Q`$ factor can actually increase for lower frequencies. In order to take into account this property of the $`Q`$ factor, we have assumed for example a relation

$$Q(m_a)=3\times 10^4(10^5\mathrm{eV}/m_a)^{2/3},$$
(137)   
providing $`Q(m_a)1.4\times 10^46.7\times 10^4`$ for $`m_a=3\times 10^5\mathrm{eV}3\times 10^6\mathrm{eV}`$. The results obtained by taking a fixed value and this relation for the $`Q`$ factor are plotted toghether by solid and dotted lines, respectively, in Figs. 11 and 12.

Here we clearly find that if the $`Q`$ factor increases as given in Eq. (137), the detection sensitivity can be improved significantly for the lower axion masses $`m_a10^6\mathrm{eV}`$. It is also observed in most cases that the the measurement time $`\mathrm{\Delta }t`$ and the total scanning time $`t_{\mathrm{tot}}`$ once become local minimum for certain axion mass around $`m_a=10^5\mathrm{eV}`$. As the axion mass gets smaller from this minimum point, $`\mathrm{\Delta }t`$ and $`t_{\mathrm{tot}}`$ increase until reaching a local maximum, and then turn to decrease again. This feature is understood by noting the behaviors of the signal and noise rates with respect to the axion mass. In fact, as shown in Eq. (132), the singal rate increases approximately as $`R_sm_a^3`$ (with fixed $`Q`$) when the axion mass gets smaller. This increase of the signal rate is eventually overwhelmed by the more rapid increase of the noise rate $`R_n\mathrm{e}^{m_a/k_\mathrm{B}T_c}`$ for $`m_a/k_\mathrm{B}T_c1`$ proportional to the thermal photon number $`\overline{n}_c`$. Then, as the axion mass gets smaller to be comparable to the cavity temperature $`T_c`$, the increase of the noise rate becomes gradually moderate, and the singnal rate begins to dominate again.

We can now conclude with these detailed numerical calculations. If the galactic axion search is made by utilizing the Rydberg atom cavity detecter, the DFSZ axion limit can be reached in frequency ranges of 10% around the axion mass $`m_a10^6\mathrm{eV}10^5\mathrm{eV}`$ for reasonable scanning times. The optimal condition for the detection sensitivity is attained by cooling the cavity down to a temperature $`T_c10\mathrm{m}\mathrm{K}`$ and adjusting the experimental parameters such as the atomic beam intensity $`I_{\mathrm{Ryd}}`$ and velocity $`v`$ and also the scanning frequency step $`\mathrm{\Delta }\omega _c`$.

## X Summary

We have developed quantum theoretical calculations on the dynamical system consisting of the cosmic axions, photons and Rydberg atoms which are interacting in the resonant cavity. The time evolution is determined for the number of Rydberg atoms in the upper state which are excited by absorbing the axion-converted and thermal background photons. The calculations are made, in particular, by taking into account the actual experimental situation such as the motion and uniform distribution of the Rydberg atoms in the incident beam and also the spatial variation of the electric field in the cavity. Some essential aspects on the axion-photon-atom interaction in the resonant cavity are clarified by these detailed calculations. Then, by using these results the detection sensitivity of the Rydberg atom cavity detector is estimated properly. The present quantum analysis clearly shows that the Rydberg atom cavity detector is quite efficient for the dark matter axion search.

###### Acknowledgements.

This research was partly supported by a Grant-in-Aid for Specially Promoted Research by the Ministry of Education, Science, Sports and Culture, Japan under the program No. 09102010.    

## A Analytic solution

Basic properties of the axion-photon-atom interaction in the cavity can be examined by considering the simple case where the atom-photon coupling does not depend on the time (i.e., the spatial variation of the electric field in the cavity is not considered). In this case, the analytic solution is obtained for the particle number matrix $`𝒩_{ij}(t)`$ as follows .

We first consider a more general case with $`K`$ atomic bunches in the incident beam where the atom-photon couplings $`\mathrm{\Omega }_1`$$`\mathrm{\Omega }_K`$, which may be defferent each other, are independent of the time. The effective Hamiltonian $``$ given in Eq. (107) with constant $`\mathrm{\Omega }_1`$$`\mathrm{\Omega }_K`$ is diagonalized by a nonsingular linear transformation $`𝒯`$ (not unitary for nonhermitian $``$):

$$\overline{}=𝒯^1𝒯=\mathrm{diag}.(\lambda _1,\lambda _2,\mathrm{},\lambda _{K+2}).$$
(A1)   
Accordingly, the relevant matrices are introduced by 

$`\overline{𝒩}^{}(t)`$ $`=`$ $`𝒯^1𝒩^{}(t)𝒯^{1\mathrm{T}},`$ (A2)   
$`\overline{𝒟}`$ $`=`$ $`𝒯^1𝒟𝒯^{1\mathrm{T}},`$ (A3)   
where $`𝒩^{}(t)=𝒰^1(t)𝒩(t)𝒰^{1\mathrm{T}}(t)`$, as given in Eq. (75), with $`𝒰(t)=\mathrm{exp}[it]`$ for the time-independent $``$. Then, the master equation (74) is reduced to a set of unmixed equations,

$$\frac{d\overline{𝒩}_{ij}^{}}{dt}=\mathrm{e}^{i(\lambda _i^{}\lambda _j)t}\overline{𝒟}_{ij}.$$
(A4)   
These equations are readily solved with an appropriate initial value $`𝒩(0)=𝒩^{}(0)`$ as

$$\overline{𝒩}_{ij}^{}(t)=\overline{𝒩}_{ij}^{}(0)+\overline{𝒞}_{ij}(t),$$
(A5)   
where

$$\overline{𝒞}_{ij}(t)=i\frac{\overline{𝒟}_{ij}}{\lambda _i^{}\lambda _j}\left[\mathrm{e}^{i(\lambda _i^{}\lambda _j)t}1\right].$$
(A6)   
By using this solution for $`\overline{𝒩}^{}(t)`$, the original $`𝒩(t)`$ is caluclated as

$$𝒩(t)=𝒰^{}(t)\left[𝒩(0)+𝒞(t)\right]𝒰^\mathrm{T}(t),$$
(A7)   
where 

$`𝒞(t)`$ $`=`$ $`𝒯^{}\overline{𝒞}(t)𝒯^\mathrm{T},`$ (A8)   
$`𝒰(t)`$ $`=`$ $`𝒯\mathrm{exp}[i\overline{}t]𝒯^1.`$ (A9)   
The above procedure with Eqs. (A5) – (A9) to find the solution may be viewed as an inhomogeneous linear transformation from $`𝒩(0)`$ to $`𝒩(t)`$ depending on the form of $``$:

$$𝒩(t)=𝒰_{}[𝒩(0)].$$
(A10)   
The analytic solution obtained in this way can be used in a good approximation for a sufficiently small time interval where the variation of the atom-photon coupling is neglected. (This is equivalent to approximate the electric field profile $`f(x)`$ with a relevant step-wise function, if the time interval is taken to be $`\delta t=\delta x/v`$.) Then, by applying the transformation (A10) for a time sequence $`tt_{z+1}>t_z>t_{z1}>\mathrm{}>t_0`$ with small intervals $`t_kt_{k1}`$ ($`1kz+1`$) the approximate solution for $`𝒩(t)`$ may be obtained as

$$𝒩(t_0)\stackrel{𝒰_{\left(t_0\right)}}{}𝒩(t_1)\stackrel{𝒰_{\left(t_1\right)}}{}\mathrm{}\stackrel{𝒰_{\left(t_z\right)}}{}𝒩(t).$$
(A11)   
We now consider the specific case of $`K=1`$ with $`\mathrm{\Omega }(t)=\mathrm{\Omega }_N`$, where the Hamiltonian is given by

$$=\left(\begin{array}{ccc}\omega _b\frac{i}{2}\gamma _b& \mathrm{\Omega }_N& 0\\ \mathrm{\Omega }_N& \omega _c\frac{i}{2}\gamma _c& \kappa \\ 0& \kappa & \omega _a\frac{i}{2}\gamma _a\end{array}\right).$$
(A12)   
Then, the transformation matrix $`𝒰(t)`$ is calculated from Eq. (A9) as

$$𝒰_{ij}(t)=\underset{k}{}g_{ij}^k\mathrm{e}^{i\lambda _kt},$$
(A13)   
with

$$g_{ij}^k=𝒯_{ik}𝒯_{kj}^1,$$
(A14)   
where the sum is not taken over the index $`k`$. This transformation matrix may also be expressed in terms of the Laplace transform $``$  as

$$𝒰(t)=^1[(s\mathrm{𝟏}+i)^1].$$
(A15)   
The coefficients $`g_{ij}^k`$ are readily found by using the Laplace transform as

$$g_{ij}^k=\underset{si\lambda _k}{lim}(s+i\lambda _k)(s\mathrm{𝟏}+i)_{ij}^1,$$
(A16)   
where in the leading orders of $`\kappa `$

$`(s\mathrm{𝟏}+i)^1`$ (A17)   
$`\mathrm{\Lambda }^1(s)\left(\begin{array}{ccc}s_as_c& i\mathrm{\Omega }_Ns_a& \kappa \mathrm{\Omega }_N\\ i\mathrm{\Omega }_Ns_a& s_as_b& i\kappa s_b\\ \kappa \mathrm{\Omega }_N& i\kappa s_b& s_cs_b+\mathrm{\Omega }_N^2\end{array}\right)`$ (A21)   
with

$$s_is+i\omega _i+\frac{1}{2}\gamma _i,$$
(A22)   
and 

$`\mathrm{\Lambda }(s)`$ $`=`$ $`\mathrm{det}(s\mathrm{𝟏}+i)`$ (A23)   
$`=`$ $`(s+i\lambda _1)(s+i\lambda _2)(s+i\lambda _3).`$ (A24)   
The respective particle numbers are then obtained with the initial value $`𝒩(0)=\mathrm{diag}.(0,\overline{n}_c,\overline{n}_a)`$ as

$$n_i(t)=𝒩_{ii}(t)=r_{ic}(t)\overline{n}_c+r_{ia}(t)\overline{n}_a,$$
(A25)   
where

$$r_{ij}(t)=\underset{m,n}{}g_{ij}^mg_{ij}^n\left[\left(1\frac{\gamma _j}{\mathrm{\Lambda }_{mn}}\right)\mathrm{e}^{\mathrm{\Lambda }_{mn}t}+\frac{\gamma _j}{\mathrm{\Lambda }_{mn}}\right]$$
(A26)   
with

$$\mathrm{\Lambda }_{mn}=i(\lambda _m^{}\lambda _n).$$
(A27)   
The eigenmodes of $``$ are readily found as follows in the limit of $`\kappa 0`$, which is indeed sufficient for $`\kappa /\gamma <10^{15}`$ or so:

$`\lambda _1`$ $`=`$ $`\omega _b{\displaystyle \frac{i}{2}}\gamma _b\mathrm{\Delta }_{bc}+\left[(\mathrm{\Delta }_{bc})^2+\mathrm{\Omega }_N^2\right]^{1/2},`$ (A28)   
$`\lambda _2`$ $`=`$ $`\omega _c{\displaystyle \frac{i}{2}}\gamma _c+\mathrm{\Delta }_{bc}\left[(\mathrm{\Delta }_{bc})^2+\mathrm{\Omega }_N^2\right]^{1/2},`$ (A29)   
$`\lambda _3`$ $`=`$ $`\omega _a{\displaystyle \frac{i}{2}}\gamma _a,`$ (A30)   
where $`\mathrm{\Delta }_{bc}\frac{1}{2}(\omega _b\omega _c)\frac{i}{4}(\gamma _b\gamma _c)`$, and $`z^{1/2}|z|^{1/2}\mathrm{exp}[i\mathrm{arg}(z)/2]`$. The eigenmodes $`\lambda _1`$ and $`\lambda _2`$, which mainly consist of the atom and photon, are determined by the equation

$$\left(\lambda \omega _b+\frac{i}{2}\gamma _b\right)\left(\lambda \omega _c+\frac{i}{2}\gamma _c\right)\mathrm{\Omega }_N^2=0.$$
(A31)   
The eignemode $`\lambda _3`$ is, on the other hand, almost identical with the axion mode. These eigenmodes of $``$ satisfy the conditions

$$\mathrm{Re}[\mathrm{\Lambda }_{mn}]=(\mathrm{Im}[\lambda _m]+\mathrm{Im}[\lambda _n])>0,$$
(A32)   
so that the respective particle numbers $`n_i(t)`$ converge to certain asymptotic values for $`t\mathrm{}`$.

## B Langevin picture

We here describe the quantum evolution of the axion-photon-atom system in the Langevin picture , which is supplementary to the treatment in the Liouville picture reproducing the same results for the quantum averages of the relevant quantities. The Langevin equation for the damped oscillators is given by

$$\frac{dq_i}{dt}=i_{ij}(t)q_j+F_i.$$
(B1)   
The external forces $`F_i`$ are introduced for the Liouvillian relaxations. These external operators obey the following correlation property :

$`F_i^{}(\tau )F_j(\tau ^{})`$ $`=`$ $`\delta _{ij}\gamma _i\overline{n}_i\delta (\tau \tau ^{}),`$ (B2)   
$`F_i^{}(\tau )q_j(0)`$ $`=`$ $`0(\tau ,\tau ^{}>0).`$ (B3)   
The solution of the Langevin equation (B1) is readily found with the transformation matrix $`𝒰(t)`$ representing the time evolution due to $`(t)`$, as given in Eq. (76):

$$q_i(t)=𝒰_{ij}(t)q_j(0)+_0^t\left[𝒰(t)𝒰^1(\tau )\right]_{ij}F_j(\tau )𝑑\tau .$$
(B4)   
Then, the time evolution of the particle numbers is calculated with this solution (B4) of the Langevin equation and the correlations (B3) as

$`𝒩_{ij}(t)=q_i^{}(t)q_j(t)`$ (B5)   
$`=𝒰_{ik}^{}(t)\left[𝒩(0)+{\displaystyle _0^t}𝒰^1(\tau )𝒟𝒰^{1\mathrm{T}}(\tau )𝑑\tau \right]_{kl}𝒰_{lj}^\mathrm{T}(t),`$ (B6)   
where $`𝒟_{ij}=\gamma _i\overline{n}_i\delta _{ij}`$. This solution for $`𝒩(t)`$ just coincides with that given in Eq. (80) which is obtained in the Liouville picture.