text:
# Effects of the trapping potential on a superfluid atomic Fermi Gas

## Abstract

We examine a dilute two-component atomic Fermi gas trapped in a harmonic potential in the superfluid phase. For experimentally realistic parameters, the trapping potential is shown to have crucial influence on various properties of the gas. Using an effective hamiltonian, analytical results for the critical temperature, the temperature dependence of the superfluid gap, and the energy of the lowest collective modes are derived. These results are shown to agree well with numerical calculations. We furthermore discuss in more detail a previous proposed method to experimentally observe the superfluid transition by looking at the collective mode spectrum. Our results are aimed at the present experimental effort to observe a superfluid phase transition in a trapped atomic Fermi gas.

Pacs Numbers: 03.75.Fi, 05.30.Fk,

With the increasingly impressive experimental ability to trap and cool atoms with fermionic statistics, a new and very active area within the field of ultracold atomic gases has emerged. Several groups are now involved in trapping and cooling atoms with Fermi statistics and temperatures as low as $`0.2T_F`$ have been obtained for 40K and 6Li . One major motivation for this experimental effort is to observe a phase transition to a superfluid state predicted to occur below a certain critical temperature $`T_c`$ for a two-component gas .

In this paper, we examine various properties of such a two-component Fermi gas with an equal number of atoms $`N_\sigma `$ of mass $`m`$ in two hyperfine states $`|\sigma =,`$ trapped in a spherically symmetric harmonic potential of the form $`U_0(r)=m\omega _T^2r^2/2`$. Assuming that the gas is in the dilute regime, the interaction between atoms in the two different hyperfine states can be well described by a contact potential $`g\delta (𝐫)`$ with $`g=4\pi \mathrm{}^2a/m`$ and $`a`$ the $`s`$-wave scattering length. The dilute regime is defined as $`k_F|a|1`$ with $`k_F=\sqrt{2m\mu _F}/\mathrm{}`$ and $`\mu _F`$ the chemical potential. For $`a<0`$ the effective atom-atom interaction is attractive and the gas is unstable toward the formation of Cooper pairs below a critical temperature $`T_c`$ . However, the specific formula for the critical temperature was derived for a homogeneous system and is not necessarily valid for a trapped gas as we shall see. Indeed, it was shown in ref.  that the nature of the superfluid state depends on the coupling strength and the number of particles trapped. For many particles trapped and/or stronger coupling strengths (keeping $`k_F|a|1`$), the size of the Cooper pairs given by the coherence length $`\xi _{BCS}=\mathrm{}^2k_F/(\pi m\mathrm{\Delta })`$, with $`\mathrm{\Delta }`$ being the $`BCS`$ pairing gap (for $`T=0`$), is much smaller than the size $`R=\sqrt{2\mu _F/m}/\omega _T`$ of the system. In this limit the system can be described as a quasi-homogeneous bulk superconductor with a locally varying chemical potential $`\mu _F(r)=\mu _FU_0(r)`$ . For weak atom-atom interactions or for few particles trapped on the other hand, the size $`\xi _{BCS}`$ of the Cooper pairs is comparable to the size of the system. The pairs are then formed between particles with time-reversed angular momentum quantum numbers $`(l,m)`$ and $`(l,m)`$ residing in the same harmonic oscillator shell and the gas is in the so-called intrashell regime . It can in this limit in many ways be regarded as a large nucleus with some properties qualitatively different from the bulk regime. The transition between the two different regimes is roughly determined by $`\xi _{BCS}R`$ or equivalently $`\mathrm{\Delta }\mathrm{}\omega _T`$. Using the bulk theory prediction $`\mathrm{\Delta }=(2/e)^{7/3}\mu _F\mathrm{exp}(\pi /2k_F|a|)`$  and the Thomas-Fermi result $`\mu _F=(6N_\sigma )^{1/3}\mathrm{}\omega _T`$, the equation $`\mathrm{\Delta }=\mathrm{}\omega _T`$ yields $`N_\sigma ^{}=\frac{1}{6}\left(\frac{e}{2}\right)^7\mathrm{exp}(3\pi /2k_F|a|)`$. For $`N_\sigma N_\sigma ^{}`$ atoms trapped, the gas is in the intrashell regime and for $`N_\sigma N_\sigma ^{}`$ atoms trapped it is in the bulk regime. With $`k_F|a|=0.1`$ this condition yields $`N_\sigma ^{}=4\times 10^{20}`$ and for $`k_F|a|=0.3`$ we get $`N_\sigma ^{}=9\times 10^6`$. Current experiments have $`N𝒪(10^5)`$ particles trapped  suggesting that in the dilute limit one is likely to be in the intrashell regime. The experiments seem to be aimed at using optical traps as the confining potential since these offer more flexibility with respect to which atoms can be trapped . These optical traps are rather tight with trapping frequencies $`\omega _T𝒪(10^3Hz)`$ yielding $`\mathrm{}\omega _T𝒪(10^2nK)`$. Thus, the condition $`k_BT_c\mathrm{}\omega _T`$ \[$`k_BT_ck_BT_F=(6N_\sigma )^{1/3}\mathrm{}\omega _T`$\] for the intrashell regime is less restrictive than for magnetic traps. However, cooling schemes as the one described in , which enables the lowering of $`T/T_F`$ by orders of magnitude, will probably have to be employed.

For $`k_F|a|0.3`$ the gas is no longer in the dilute regime. Recently, a number of papers have been published examining the gas in this strongly correlated limit . It is presently unclear however, whether the gas is stable toward spinodal decomposition or other instabilities and it is outside the scope of the present paper to describe the gas in this regime.

One property of great experimental importance is the critical temperature $`T_c`$ which we will now estimate for the gas in the intrashell regime. In this limit there is only pairing between particles with quantum numbers $`(n,l,m)`$ and $`(n,l,m)`$ with the unperturbed harmonic oscillator shell energy $`ϵ_n=(n+3/2)\mathrm{}\omega _T`$. This fact, combined with a number of semiclassical approximations allows one to derive an effective Hamiltonian of the form

$$\widehat{H}_{eff}=\underset{nlm\sigma }{}\xi _n\widehat{a}_{nlm\sigma }^{}\widehat{a}_{nlm\sigma }$$
(1)
$$\frac{2\mathrm{}\omega _TG}{\mathrm{\Omega }_{n_F}}\underset{nlm,n^{}l^{}m^{}}{}(1)^{m+m^{}}\widehat{a}_{nlm}^{}\widehat{a}_{nlm}^{}\widehat{a}_{n^{}l^{}m^{}}\widehat{a}_{n^{}l^{}m^{}}$$
(2)
where $`\mathrm{\Omega }_{n_F}=(n_F+1)(n_F+2)/2`$ is the unperturbed degeneracy of the highest harmonic oscillator shell occupied with energy $`ϵ_{n_F}`$ (ignoring the Hartree field for the moment), $`\widehat{a}_{nlm}`$ is the operator removing a particle in the shell with energy $`ϵ_n`$ and angular momentum quantum numbers $`l,m`$ and $`\xi _n=ϵ_n\mu _F`$. The effective coupling strength giving rise to the pairing correlations is given by $`G=32\pi ^2k_F|a|/15`$. From Eq. (1) a number of analytical results can be derived. The gap equation becomes :

$$\frac{\mathrm{\Delta }}{\mathrm{}\omega _T}=G\underset{n}{\overset{n2n_F}{}}\frac{\mathrm{\Delta }}{E_n}[12f(E_n)],$$
(3)
where we have as a first approximation to the pseudopotential scheme introduced the cutoff $`n2n_F`$. The quasiparticle energy is $`E_n=(\xi _n^2+\mathrm{\Delta }^2)^{1/2}`$ and $`f(x)=[\mathrm{exp}(\beta x)+1]^1`$ with $`\beta =1/k_BT`$. As $`\mathrm{\Delta }\mathrm{}\omega _T`$ in the intrashell regime, shell effects are pronounced; the behaviour of the system depends crucially on whether the system is in a filled shell configuration or whether the highest occupied shell is only partially filled . Since the filled shell configuration occurs only for certain “magic” numbers of particles trapped, we expect the typical experimental situation to be that of a partially filled highest shell \[$`\mu _F(n_F+3/2)\mathrm{}\omega _T`$\]. In this case, and for $`k_BT_c\mathrm{}\omega _T`$ appropriate for the intrashell regime, Eq. (3) yields the critical temperature

$$\frac{k_BT_c}{\mathrm{}\omega _T}\frac{G}{24\mathrm{ln}(e^\gamma n_F)G}=\frac{\mathrm{\Delta }_0}{2\mathrm{}\omega _T}.$$
(4)
Here $`\gamma 0.577`$ is Euler’s constant and $`\mathrm{\Delta }_0`$ is the gap for $`T=0`$ in the intrashell regime . For $`G0`$, we have $`k_BT_c/\mathrm{}\omega _T=G/2`$. In this very weak interaction limit, there is only pairing in the partially filled shell right at the chemical potential ($`n=n_F`$) and the shells away from the chemical potential are unaffected by the pairing. For increasing $`G`$, the term $`\mathrm{ln}(e^\gamma n_F)G`$ in Eq. (4) becomes important and $`k_BT_c/\mathrm{}\omega _T>G/2`$. Physically, this corresponds to the interaction now being so strong that there is significant pairing in the shells away from the chemical potential; each Cooper pair is however still only formed between particles within the same shell. We see that $`T_c`$ in the intrashell regime given by Eq. (4) is very different from $`k_BT_c\mu _F\mathrm{exp}(\pi /2k_F|a|)`$ relevant for the bulk regime ; for $`G1`$, Eq. (4) yields a critical temperature which is orders of magnitude higher than the bulk result. This stems from the fact that the shell structure of the trap levels increases the density of states for certain energies. For $`\mathrm{ln}[\mathrm{exp}(\gamma )n_F]G=1/2`$, Eq. (4) diverges but this is clearly unphysical as the intrashell limit is characterized by $`k_BT_c\mathrm{}\omega _T`$. For such large values of $`G`$ the intrashell pairing ansatz is invalid and Cooper pairs are formed between particles residing in different shells. It should be noted that Eq. (4) does not include the effect of induced interactions coming from density oscillations of the gas. In the homogeneous case, this effect is known to reduce $`T_c`$ by a factor $`2.2`$  and it should be investigated which influence this effect has in the intrashell regime as well. However, it must be expected that the critical temperature including the induced interactions will be of the same order of magnitude as Eq. (4) and thus very different from the bulk theory prediction.

To check the validity of Eq. (4), we plot in Fig. 1 $`k_BT_c/\mathrm{}\omega _T`$ as a function of $`G`$. The solid curve is Eq. (4) and the $`\times `$’s are obtained from a numerical calculation solving the Bogoliubov-de Gennes (BdG) equations derived from the full Hamiltonian using a pseudopotential scheme . The dashed line depicts the bulk prediction $`k_BT_c\mu _F\mathrm{exp}(\pi /2k_F|a|)`$. We have chosen $`n_F30`$ giving $`N=2N_\sigma 10^4`$ particles trapped. There is excellent agreement between the numerical results and Eq. (4) for $`G0.05`$. We conclude in agreement with ref.  that the intrashell description works well for $`k_BT_c\mathrm{}\omega _T`$. In this limit on the other hand, the bulk theory is qualitatively wrong as expected. For stronger coupling, the nature of the superfluid state undergoes a transition and the intrashell description breaks down. There is now pairing between particles residing in different harmonic oscillator shells. This breakdown of the intrashell theory occurring for $`k_BT_c𝒪(\mathrm{}\omega _T)`$ for larger coupling strengths $`G`$ is not depicted in Fig. 1. Eventually the coherence length $`\xi _{BCS}`$ becomes much smaller than the size of the system. As explained above, for $`N10^4`$ trapped the gas is no longer in the dilute limit when $`\xi _{BCS}R`$.

The Hartree field introduces a dispersion on the Harmonic oscillator levels with respect to the angular momentum $`l`$ lifting the shell degeneracy $`\mathrm{\Omega }_n=(n+1)(n+2)/2`$. When $`n_F30`$ (depending on the size of the coupling strength), this splitting is small compared to the pairing energy $`\mathrm{\Delta }`$ and it can be ignored . For such systems, Eq. (4) for $`k_BT_c`$ should be accurate. For very large systems with $`n_F200`$ on the other hand, the Hartree splitting of the shell is much larger than $`\mathrm{\Delta }`$ and there is for weak coupling only pairing within states with a fixed $`n`$ and $`l`$ . In order to examine the effect of the Hartree field, we plot in Fig. 1 as $`+`$’s $`k_BT_c/\mathrm{}\omega _T`$ as a function of $`G`$ obtained from a numerical calculation where the Hartree field is included. We see that effect of the Hartree field is to suppress $`T_c`$ for small $`G`$ due to the fact that the lifting of the shell degeneracy reduces the density of states at the chemical potential. For stronger coupling, the Hartree field increases $`T_c`$ as it in this limit enhances the spatial overlap between the orbitals of the Cooper pairs. However, from Fig. 1 we see that Eq. (4) agrees qualitatively with the numerical results of the critical temperature. This is because $`n_F30`$ is not quite in the limit where the Hartree field can be neglected, but also far away from the regime $`n_F200`$ where there is only pairing within levels with a single $`l`$-value. Again, these results are in agreement with earlier calculations on the $`T=0`$ properties of the gas .

We now examine the temperature dependence of the pairing energy $`\mathrm{\Delta }`$ in the intrashell regime. Solving Eq. (3) for $`\mathrm{\Delta }(T)/\mathrm{}\omega _T1`$, we obtain

$$\frac{T}{T_c}=\frac{\mathrm{\Delta }(T)}{k_BT_c\mathrm{ln}\left[\frac{1+\mathrm{\Delta }(T)/2k_BT_c}{1\mathrm{\Delta }(T)/2k_BT_c}\right]}$$
(5)
from which $`\mathrm{\Delta }(T)`$ can be found by inversion. For $`(T_cT)/T_c1`$ such that $`\mathrm{\Delta }(T)/k_BT_c1`$ we have

$$\frac{\mathrm{\Delta }(T)}{k_BT_c}=\sqrt{12}\sqrt{1T/T_c}.$$
(6)
In Fig. 2, we plot the lowest quasiparticle energy $`E(T)/K_BT_c=\mathrm{\Delta }(T)/k_BT_C`$ (since $`\xi _{n_Fl}=0`$) as a function of $`T/T_c`$. The solid line is Eq. (5) and the $`\times `$’s and $`+`$’s are numerical results from a solution of the BdG Eq. with $`10^4`$ particles trapped and $`G=0.067`$ where the effect of the Hartree field has been excluded for the $`\times `$’s.

We see that there is good agreement between the numerical results and Eq. (5) in the intrashell regime. In particular, we have $`\mathrm{\Delta }(T=0)/k_BT_c2`$ as was predicted in Eq. (4). Again, the effect of the Hartree field does not change the qualitative behaviour of the system for this set of parameters since $`n_F200`$. So we conclude that Eq. (5) describes $`\mathrm{\Delta }(T)`$ well in the intrashell regime.

An important question is how to experimentally observe the superfluid transition. One way is to look at the collective mode spectrum of the gas on which superfluidity should have several observable effects . When the gas is in the intrashell regime, two such effects are the emergence of low energy density oscillation modes of monopole and quadrupole symmetry. The energy of these new modes should scale as $`\mathrm{}\omega 2\mathrm{\Delta }\mathrm{}\omega _T`$ . To examine this, we plot in Fig. 3 as $`+`$’s the frequency $`\omega _Q`$ of the low energy quadrupole mode as a function of $`T`$. The collective response is calculated using a RPA-scheme including the superfluid correlations as described in ref. . The solid line gives $`2\mathrm{\Delta }(T)`$ obtained from Eq. (5) and we also for comparison plot as $`\times `$’s $`2\mathrm{\Delta }(T)`$ obtained from a numerical solution to the BdG eqn. The parameters used are the same as for Fig. 2 and the effect of the Hartree field is included.

We see that there is fairly good agreement between the numerical results from the RPA calculation and Eq. (5). The energy of the collective mode given within the RPA is approximately $`2\mathrm{\Delta }`$ indicating that the mode consists of breaking a Cooper pair. The fact that the collective response frequency in general is lower than the pair breaking prediction $`2\mathrm{\Delta }`$ can be understood as follows: The excitation consists of the breaking of a Cooper pair of monopole symmetry costing $`2\mathrm{\Delta }`$ in energy. However, the broken pair of particles then form a new Cooper pair with quadrupole symmetry; the difference $`2\mathrm{\Delta }\mathrm{}\omega _Q`$ can be interpreted as the binding energy of the quadrupole Cooper pair. The excitation can therefore be understood as the breaking of a Cooper pair with monopole symmetry and the subsequent formation of a Cooper pair of quadrupole symmetry, i.e. a pair vibration mode.

The frequency $`\omega _M`$ of the lowest monopole mode is approximately the same as for the lowest quadrupole mode shown in Fig. 3. Again, the excitation frequency is lower than the simple pair breaking prediction $`2\mathrm{\Delta }`$ because the broken pair form a new Cooper pair in an excited monopole state. We will in a future publication examine in more detail the binding energy of these excited Cooper pairs of monopole and quadrupole symmetry giving rise to the collective mode frequencies being somewhat lower than $`2\mathrm{\Delta }`$. For the purpose of the present paper however, it is adequate to conclude that the onset of superfluidity below $`T_c`$ indeed does give rise to new density oscillation modes of monopole and quadrupole symmetry with energy $`\mathrm{}\omega 2\mathrm{\Delta }`$ in the intrashell regime.

The collective modes can be excited by modulating the trapping frequency producing a perturbation with the right symmetry . There is however one important point concerning the monopole pair vibration mode: In order to excite this mode, one should apply a perturbation of the trapping potential $`\delta \widehat{H}=\mathrm{sin}(\omega t)_\sigma d^3rF(𝐫)\widehat{\psi }_\sigma ^{}(𝐫)\widehat{\psi }_\sigma (𝐫)`$ with $`F(𝐫)\alpha r^2`$. This is because in the intrashell limit, the Cooper pairs are only formed between states within the same harmonic oscillator shell with unperturbed energy $`(n+3/2)\mathrm{}\omega _T`$. Neglecting the effect of the Hartree field, the states within one shell all have the same expectation value for $`r^2`$. Thus, due to the orthogonality between the eigenstates of $`\widehat{H}`$, a perturbation $`r^2`$ (such as the modulation of the trapping frequency) will not probe states which correspond to excitations of pairs within the same harmonic oscillator shell and the monopole pair vibration mode $`\mathrm{}\omega _M2\mathrm{\Delta }`$ will not be excited. To excite this mode, the scheme using off-resonant lasers demonstrated might be useful . For the numerical calculations, the monopole mode was excited by a perturbation with $`F(r)r^4`$. A related issue is the spectral weight with which the pair vibration modes are excited as compared to the normal phase modes at $`\omega 2\omega _T`$. The spectral weight $`S(\omega )_\nu |0|\delta \widehat{H}|0|^2\delta (\mathrm{}\omega E_0E_\eta )`$ with $`|0`$ and $`|\nu `$ being the ground and excited states with energies $`E_0`$ and $`E_\nu `$ respectively, is for $`T=0`$ in the quasiparticle approximation

$$S(\omega )\underset{\eta \eta ^{}}{}|d^3r[u_\eta v_\eta ^{}+v_\eta u_\eta ^{}]F|^2\delta (\mathrm{}\omega E_\eta E_\eta ^{}),$$
(7)
where the elementary quasiparticles with excitation energies $`E_\eta `$ are described by the Bogoliubov wave functions $`u_\eta (𝐫)`$ and $`v_\eta (𝐫)`$ . $`\eta =(n,l,m)`$ is shorthand for the quantum numbers characterizing the quasiparticle. The collective pair vibration is predominantly built out of quasiparticle excitations with $`n=n^{}=n_F`$ whereas the normal phase excitation $`\omega 2\omega _T`$ is built out of excitations with $`n=n_F2,n_F1,n_F`$ and $`n^{}=n_F+2`$ or vice versa in Eq. (7). Now, for the $`n=n_F`$ shell, even a minute pairing will mix the particles and holes yielding $`u_{n_Flm}(𝐫)=v_{n_Flm}(𝐫)=2^{1/2}\varphi _{n_Flm}^0(𝐫)`$ where $`\varphi _{n_Flm}^0(𝐫)`$ are the non-interacting harmonic trap eigenfunctions with energy $`(n+3/2)\mathrm{}\omega _T`$ . This makes the $`u_\eta v_\eta ^{}`$ in Eq. (7) large and depending on the spatial overlap between $`F(𝐫)`$ and the excitations involved, the strength of the pairing vibration will be of the same order as the normal phase excitation $`\omega 2\omega _T`$, although in general somewhat smaller since it is built from excitations from only the $`n=n_F`$ shell whereas the $`\omega 2\omega _T`$ mode is built from the shells $`n_F2,n_F1,n_F`$. This is illustrated in the inset in Fig. 3, where we plot $`S(\omega )`$ (calculated using the RPA) for $`T=0`$ with the parameters chosen above for the quadrupole mode excited by $`F(𝐫)=r^2Y_{20}(\theta ,\varphi )`$. We see that the strength of the pairing vibration is comparable to the normal phase mode thereby making it experimentally observable.

When the trap is not spherical, the single particle energies and wave functions are different making the quantitative predictions of this paper invalid. However, if all trapping frequencies $`\omega _i`$ ($`i=x,y,z`$) are such that $`\mathrm{}\omega _ik_BT_c`$, the qualitative conclusions of this paper still holds: $`T_c`$ can be orders of magnitude higher than the bulk result due to increased density of states for certain energies and there will be pair vibrations with $`\mathrm{}\omega 2\mathrm{\Delta }`$. If, on the other hand, one frequency is loose (say $`\mathrm{}\omega _zk_BT_c`$), then the gas will be in some interesting mixture of the bulk and the intrashell regime which warrants further study.

In conclusion, we have shown that for experimentally realistic parameters, the trapping potential has a significant influence on the properties of a trapped atomic Fermi gas in the superfluid phase in the dilute limit. In the intrashell regime, we used an effective Hamiltonian to derive analytical results concerning the critical temperature, the temperature dependence of the pairing gap and the lowest collective oscillations of monopole and quadrupole symmetry. These results were shown to agree well with numerical calculations. We also discussed in further detail a previously proposed scheme to detect the onset of superfluidity with decreasing temperature by looking at the collective mode spectrum. Much of the analysis of this paper should be relevant for the present experimental effort aimed at observing the superfluid transition in a two-component trapped atomic Fermi gas.

We acknowledge valuable discussions with C. Salomon, B. Mottelson, and J. R. Krumrine.


format: text

text:
# Investigation of the role of neutron transfer in the fusion of 32,34S with 197Au,208Pb using quasi-elastic scattering

## 1 Introduction

Measured cross sections for heavy ion fusion at energies below the Coulomb barrier show strong isotopic dependences and exceed theoretical predictions based on a single barrier penetration model by several orders of magnitude . This has been observed for a wide range of systems and is understood to arise from the coupling of the relative motion of the interacting nuclei to their rotational and vibrational states or to particle transfer channels. The coupling gives rise to a distribution of fusion barriers $`D(E)`$ . Experimentally, a representation $`D^{fus}(E)`$ of this distribution can be extracted from precision measurements of fusion excitation functions $`\sigma ^{fus}(E)`$ using :

$$D^{fus}(E)=\frac{d^2(E\sigma ^{fus})}{dE^2}$$
(1)
At energies below the Coulomb barrier quasi-elastic scattering excitation functions $`d\sigma ^{qel}`$/$`d\sigma ^R(E)`$, measured at backward angles, have been found  to be another suitable means to extract representations $`D^{qel}(E)`$ of the distribution $`D(E)`$ with:

$$D^{qel}(E)=\frac{d}{dE}\left(\frac{d\sigma ^{qel}}{d\sigma ^R}(E)\right)$$
(2)
In this technique, which is generally less complex than detailed fusion measurements, quasi-elastic scattering is understood to comprise elastic and inelastic scattering and also particle transfer channels. Experiments which have employed this approach have recently been carried out by several other groups .

In previous work  it has been clearly demonstrated that the quasi-elastic scattering representation $`D^{qel}(E)`$ is not identical to the representation $`D^{fus}(E)`$ extracted from fusion data, although it appears that this has not been appreciated in all studies. In particular, $`D^{qel}(E)`$ has been shown to be insensitive to the high energy part of the barrier distribution $`D(E)`$. The quasi-elastic scattering representations are therefore most useful for investigating couplings, which produce signatures in the low energy part of $`D(E)`$. This is the case when the relative motion of the two nuclei couples to positive Q-value channels .

Indeed, the comparison of the representations $`D^{fus}(E)`$ and $`D^{qel}(E)`$ for the systems 40Ca + 90,96Zr has shown that the effect of positive Q-value neutron transfer channels on the fusion dynamics are clearly seen in the representation $`D^{qel}(E)`$ . The low-lying collective states in the two Zr isotopes have very similar excitation energies and deformation parameters $`\beta _2`$, and the main differences between these two systems are in their Q-values for neutron transfer. In the heavier system the calcium nucleus can pick-up as many as eight neutrons in transfer reactions with positive Q-value, whereas the equivalent channels in the lighter system all have negative Q-values. This pronounced difference has been found to be reflected in both types of representations, $`D^{fus}(E)`$ and $`D^{qel}(E)`$, measured for these systems . The straight-forward measurement of quasi-elastic scattering excitation functions thus appears to be a promising tool to investigate the role of positive Q-value transfer channels in fusion.

The fusion reactions of the sulphur projectiles 32,34S with 197Au and 208Pb are a suitable test case for this new experimental approach to the dynamics of fusion. This is apparent from Table 1, which shows the Q-values for the pick-up of one and two neutrons for these systems. Also shown are the equivalent Q-values for the reactions 36S + 197Au, 208Pb. It is apparent that with decreasing projectile mass the Q-values progessively favour the neutron pick-up channels. The 32S and 34S projectile nuclei have similar structure, with the lowest $`2^+`$ states not being very different in terms of excitation energy and deformation parameter (see Table 2).

The systems 32,34S + 197Au, 208Pb are intermediate in mass between a large number of lighter fusion reactions, which have been well-studied using the coupled-channels framework , and the more massive systems employed for the synthesis of super-heavy elements . Results may thus indicate, if the representation $`D^{qel}(E)`$ is also applicable to this important latter group of fusion reactions, for which multiple neutron transfer may lead to macroscopic effects such as neutron-flow or neck-formation.

This paper presents detailed measurements at backward angles of quasi-elastic scattering excitation functions for the two pairs of reactions 32,34S + 197Au and 32,34S + 208Pb, from which representations $`D^{qel}(E)`$ of the fusion barrier distribution have been extracted.

## 2 Experimental Method

The experiments were performed with 32,34S-beams from the 14UD Pelletron accelerator at the Australian National University in the energy range $`E_{lab}=`$ 90.0–180.0 MeV. Three different self-supporting Au targets were used, with thicknesses in the range 140–170 $`\mu `$g/cm2. The 208PbS target was 140 $`\mu `$g/cm2 thick, evaporated onto a $``$ 20 $`\mu `$g/cm2 carbon backing. The target thickness was determined by measuring the energy loss of elastically scattered projectiles in the target at a backward angle.

A schematic diagram of the experimental setup is shown in Figure 1. In the initial study of the 32S + 208Pb reaction  quasi-elastic scattering was detected at a scattering angle of $`\theta _{lab}=170^{}`$. An energy loss signal $`\mathrm{\Delta }E`$ was measured with a gas ionisation detector. The gas detector was backed by a silicon surface barrier detector which detected the residual energy $`E_{res}`$ of the scattered nuclei. The $`\mathrm{\Delta }E`$ signal allowed the separation of the charged particle transfer contributions to the quasi-elastic scattering yield. Since this separation was not required for the interpretation of the data, in the other experiments the $`(\mathrm{\Delta }EE_{res})`$ detector telescope was replaced with a single silicon surface barrier detector at $`\theta _{lab}=159^{}`$.

For normalisation purposes, two silicon surface barrier detectors were placed at scattering angles $`\theta _{lab}=\pm 30^{}`$ to measure Rutherford scattering of projectiles. For some preliminary measurements, instead of these two monitor detectors, a readily available gas-ionisation detector  at $`\theta _{lab}=30.35^{}`$ was employed. Both setups gave consistent results, so that the data have been combined.

Figure 2 shows typical energy spectra from the backward silicon detector for the system 32S + 197Au. The spectra for the other three systems are similar. At low beam energies the spectrum only shows a well-defined peak of elastically scattered sulphur nuclei. With increasing energy the peak gradually develops a low energy tail as the yield of non-elastic scattering events rises and elastic scattering is diminished. At the higher energies fission fragments from quasi-fission and fission following compound nucleus formation are also detected. The energy window chosen for the integration of the quasi-elastic scattering yield is indicated in the figure.

For all experiments the quasi-elastic scattering yield, measured at the backward angle, comprising the sum of elastic, inelastic and transfer events, was divided by the Rutherford scattering yield detected at forward angles. These ratios have been normalised to unity at energies well below the Coulomb barrier, where only Rutherford scattering is observed. The normalised ratios thus represent the quasi-elastic scattering excitation function $`\frac{d\sigma ^{qel}}{d\sigma ^R}(E)`$, where the indices $`qel`$ and $`R`$ indicate quasi-elastic scattering and Rutherford scattering, respectively.

The measured quasi-elastic scattering excitation functions are shown in Figure 3. With the exception of the highest energies, where the quasi-elastic scattering yield is low, the statistical uncertainty is better than $`1\%`$. The energy scale has been adjusted for energy loss in the target, and the centrifugal energy corresponding to the respective detection angle has been subtracted, as described in . The excitation functions decrease smoothly with energy. The apparent energy shifts between the four sets of data reflect the expected differences in Coulomb barrier height.

## 3 Discussion of the Experimental Data

The barrier distribution representations $`D^{qel}(E)`$ have been extracted from the quasi-elastic scattering excitation functions $`\frac{d\sigma ^{qel}}{d\sigma ^R}(E)`$ by differentiation with respect to energy according to Equation (2). A point-difference formula with discrete energy steps in the range $`\mathrm{\Delta }E_{lab}=`$3–6 MeV was used to evaluate the differential. The data sets obtained for the different energy steps are consistent, so that they have been combined. The resulting barrier distribution representations $`D^{qel}(E)`$ are shown for all four systems in Figure 4. In order to facilitate a direct comparison, the energy scales of the barrier distribution representations have been normalised by dividing by an average barrier $`B_0`$, which was chosen as the energy where $`\frac{d\sigma ^{qel}}{d\sigma ^R}(E)=0.5`$. Although this determination of the average barriers may seem somewhat arbitrary, the values obtained are realistic. For example, the average barrier for the system 32S + 208Pb has been determined as 144.4 MeV by fitting the high energy part of the fusion excitation function using a single barrier penetration model . This compares well with the value of $`B_0`$=144.2 MeV used here.

As emphasized in the introduction, above the average barrier energy $`B_0`$ the representation $`D^{qel}(E)`$ is not sensitive to the fusion dynamics. Indeed, at these high energies $`D^{qel}(E)`$ is the same for all four systems. The low energy parts of the measured representations $`D^{qel}(E)`$ are discussed below.

For the system 34S + 208Pb (open circles in Figure 4 (top)) the slope of $`D^{qel}(E)`$ over the energy range $`0.92<E/B_0<0.99`$ is steeper than that for 32S + 208Pb (filled circles in Figure 4 (top)). Also, the maximum of $`D^{qel}(E)`$ for the reaction 34S + 208Pb is 0.08 MeV-1, whereas the maximum for the lighter system is only about 0.06 MeV-1. Since the integral of $`D(E)`$ is unity, this implies that the barrier distribution for 34S + 208Pb is narrower than that for 32S + 208Pb. This is consistent with significant coupling to positive Q-value neutron transfer channels in the 32S + 208Pb fusion reaction. Indeed, both the one neutron ($`Q=+1.3`$ MeV) and two neutron transfer ($`Q=+6.0`$ MeV) Q-values for this system are positive.

The equivalent data for the 197Au target (Figure 4 (bottom)) do not show such a pronounced difference. It is apparent from Figure 5 (top) that the representations $`D^{qel}(E)`$ for 32S + 197Au and 32S + 208Pb agree, which would be consistent with these systems having similar barrier distribution and thus equivalent coupling interactions. However, the comparison of the representations $`D^{qel}(E)`$ in Figure 5 (bottom) for 34S + 197Au and 34S + 208Pb demonstrates that for the gold systems $`D^{qel}(E)`$ is already broad for the heavier sulphur projectile, for which coupling to neutron transfer is less favoured. This suggests that coupling to states in the 197Au nucleus generates barrier strength at low energies, which is absent for the reactions with 208Pb.

The experimental results for the two lead systems are consistent with those reported from fusion measurements for the two reactions 32,36S + 110Pd , where additional barrier strength at low energies was also found for the lighter projectile 32S. While the new data support an important role of positive Q-value neutron transfer channels in the fusion of 32S + 208Pb, such an interpretation is only unique, if the properties of the collective states in 32S and 34S are identical, or at least can be assumed to be very similar. Recent results for the fusion of the sulphur nuclei 32,34S with 89 show that in that case the different collectivity of their quadrupole excitations (Table 2) results in a broader fusion barrier distribution for 32S than for 34S. Thus the differences observed in this work between the barrier distributions for 32S + 208Pb and 34S + 208Pb may not be solely due to coupling to the positive Q-value neutron pick-up channels. Measurements for the heavier system 36S + 208Pb may shed additional light on the fusion mechanism. In this latter system the Q-values for both one neutron and two neutron transfer are negative (see Table 1), so that any effects due to neutron pick-up can be ruled out.

## 4 Conclusions

The experiments reported here have demonstrated that precision measurements of quasi-elastic scattering at backward angles are able to probe the fusion barrier distribution of heavy systems below the average barrier. Such measurements are thus in principle sensitive to the effects of positive Q-value transfer channels. Indeed it was found that neutron transfer may affect the fusion of 32S with 208Pb. The results, however, are also consistent with the observed additional barrier strengths at low energies being due to the increased collectivity of 32S, when compared with 34S. For the fusion reactions of 32,34S with 197Au couplings related to the nuclear structure of 197Au appear to be dominant in shaping the low energy part of the barrier distibution.

Since quasi-elastic scattering experiments are generally not as complex as fusion measurements, they are well suited to survey a number of reactions to determine good candidates for detailed studies of the fusion dynamics. The extracted representations of the barrier distribution can be indicative of important coupling interactions, however, the conclusive identification of these couplings may require the measurement and interpretation of the fusion excitation function.

## Acknowledgements

The authors are grateful to the late Prof. Trevor Ophel for his contributions to these experiments and would like to thank Dr David Hinde for indepth discussions of the results. The support of Dr Jack Leigh and Dr Clyde Morton is also acknowledged.


text:
# Untitled Document

THE ASYMPTOTIC NUMBER OF ATTRACTORS

IN THE RANDOM MAP MODEL

David Romero and Federico Zertuche

Instituto de Matemáticas, UNAM

Unidad Cuernavaca, A.P. 273-3

62251 Cuernavaca, Morelos, Mexico.

david@matcuer.unam.mx zertuche@matcuer.unam.mx

Abstract

The random map model is a deterministic dynamical system in a finite phase space with $`n`$ points. The map that establishes the dynamics of the system is constructed by randomly choosing, for every point, another one as being its image. We derive here explicit formulas for the statistical distribution of the number of attractors in the system. As in related results, the number of operations involved by our formulas increases exponentially with $`n`$; therefore, they are not directly applicable to study the behavior of systems where $`n`$ is large. However, our formulas lend themselves to derive useful asymptotic expressions, as we show.

Short Title: The Asymptotic Number of Attractors

Keywords: Cellular Automata, Random Graphs, Random Maps, Binary Systems.

PACS Numbers: 89.75.-k, 45.05.+x

1. Introduction

Since the 70’s the random map model has attracted the attention of physicists . Indeed, in 1987 Derrida and Flyvbjerg have shown that the random map model is equivalent to the Kauffman model of cellular automata when the number of connections among the automata goes to infinity . This enlarged its application in the realm of theoretical biology, disordered systems, and cellular automata, for possible approaches of DNA replication, cell differentiation, and evolution theory .

On the other hand, one half of a century ago, some mathematicians had approached the random map model in the context of random graphs. First, in 1953, Metropolis & Ulam posed the problem of determining the number $`\mathrm{\Theta }(n)`$ of expected connected components (i.e., attractors) in random graphs with $`n`$ nodes ; at the time, $`\mathrm{\Theta }(n)`$ was estimated to be of order $`\mathrm{log}n`$. Only one year later, Kruskal elegantly solved the problem obtaining an exact formula together with its corresponding asymptotic behavior :

$$\mathrm{\Theta }(n)=\underset{k=1}{\overset{n}{}}\frac{n!}{\left(nk\right)!kn^k},$$
$`(1.1)`$
$$\mathrm{\Theta }(n)\frac{1}{2}\left(\mathrm{ln}2n+\gamma \right)+ϵ,\mathrm{for}n>>1,$$
$`(1.2)`$
where $`ϵ`$ vanishes for $`n\mathrm{}`$, and $`\gamma `$ is the Euler-Mascheroni constant.

The statistical distribution of the number of connected components was addressed both by Rubin & Sitgreaves , and Folkert in his Ph.D. thesis under Leo Katz supervision . Unfortunately, these two practically out of reach works were never published even though their relevance and an offer by Katz to do so . Later on, Harris partially reviewed and enlarged these results, proposing a new combinatorial expression for the distribution of connected components, as well as the complete solution of the simpler case in which the random map is one-to-one . The mathematical expressions of a general random map found by Folkert and Harris, however correct, are based on constrained sums with a number of terms of order en, and involve Stirling numbers of the first kind \[11-13\]. Therefore, unfortunately, their straightforward use in physics or biology appears quite limited (for example, it is typical to deal with $`n2^{100}`$ in models such as those of cellular automata).

In spite of its importance, and as far as we know, a study of the variance for the distribution of connected components has not yet been undertaken.

In this paper we propose a still new combinatorial formula —equivalent to the previous ones, of course— for the statistical distribution of the number of connected components. As in earlier results, it also relies on a constrained sum, and the involved computational effort increases exponentially with $`n`$. Nevertheless, it has the advantage of allowing us the derivation of the long-needed asymptotic formula for the statistical distribution. Furthermore, we easily deduce from it asymptotic formulas for the corresponding average and variance.

The paper is organized as follows. Section 2 is devoted to the main definitions of the random map model, and to settle our conventions. Then, in Section 3, we determine an exact combinatorial expression for the statistical distribution of the number of connected components in the model. The corresponding asymptotic formula for this distribution is derived in Section 4, as well as asymptotic formulas for the average and the variance for the number of connected components. To end, we present our conclusions in Section 5.

2. The Random Map Model and Functional Graphs

Let $`\mathrm{\Omega }=\{1,2,\mathrm{},n\}`$ be a set of $`n`$ points. To each point in $`\mathrm{\Omega }`$ assign at random one point in $`\mathrm{\Omega }`$ with uniform probability distribution, thus defining a function $`f:\mathrm{\Omega }\mathrm{\Omega }`$. In this way a dynamical system has been established on the so-called phase space $`\mathrm{\Omega }`$ through the iterations of $`f`$; this is the random map model .

Since $`\mathrm{\Omega }`$ is finite every orbit of $`f`$ will eventually end in a periodic attractor, and several questions are in order, like: what is the expected number of attractors in the system? which is the statistical distribution of the number of attractors? how large is the dispersion of the distribution? Here, we answer these questions starting from combinatorial arguments.

For each function $`f`$ on $`\mathrm{\Omega }`$ define a functional graph whose nodes are precisely the elements of $`\mathrm{\Omega }`$; moreover, if $`f(i)=j`$ then draw a directed link from node $`i`$ to node $`j`$, and whenever $`f(i)=i`$ a loop on node $`i`$ is drawn. As an example, Figure 1 shows a functional graph with three connected components (i.e., three attractors) in a set with $`n=11`$.

Note that each function $`f`$ on $`\mathrm{\Omega }`$ (or functional graph) can be represented by an $`n\times n`$ binary matrix $`M=\{M_{ij}\}`$, where $`M_{ij}=1`$ whenever $`f(i)=j`$ and $`M_{ij}=0`$ otherwise. Every row of matrix $`M`$ has $`n1`$ zeros and one ‘1’. Clearly, there exist $`n^n`$ such matrices, and thus $`n^n`$ is the number of functions $`f`$ that can be defined on $`\mathrm{\Omega }`$, as well as the number of distinct functional graphs on $`n`$ nodes.

 A functional Graph with three connected components for $`n=11`$.

3. The Distribution of Connected Components

What is the number $`a_n`$ of connected functional graphs (i.e., having precisely one connected component) that can be found among the $`n^n`$ functional graphs on $`n`$ nodes? Through clever combinatorial arguments Katz obtained

$$a_n=(n1)!\underset{k=0}{\overset{n1}{}}\frac{n^k}{k!},$$
$`(3.1)`$
and showed that, for $`n>>1`$, this figure grows asymptotically like

$$a_nn^n\sqrt{\frac{\pi }{2n}}\left(1+𝒪\left(\frac{1}{n}\right)\right).$$
$`(3.2)`$
Using (3.1) we derive now a new expression for the statistical distribution of the number of connected components in functional graphs. Consider first a partition $`𝒫`$ of $`\mathrm{\Omega }`$, in $`k`$ disjoint subsets $`\mathrm{\Omega }_1,\mathrm{},\mathrm{\Omega }_k`$, with $`n_1,\mathrm{},n_k`$ points, respectively; then

$$\underset{j=1}{\overset{k}{}}n_j=n,\text{with}\text{ }1n_jn,\text{for}\text{ }j=1,\mathrm{},k.$$
$`(3.3a)`$
Obviously, $`_{j=1}^ka_{n_j}`$ is the number of functional graphs $`G`$ with $`k`$ connected components, such that the subset $`\mathrm{\Omega }_j`$ defines a connected component of $`G`$, for $`j=1,\mathrm{},k`$. Moreover, for given integers $`n_1,\mathrm{},n_k`$, satisfying (3.3a), the multinomial coefficient $`\left(\genfrac{}{}{0pt}{}{n}{n_1,\mathrm{},n_k}\right)=\frac{n!}{n_1!\mathrm{}n_k!}`$ yields the number of distinct ways we can distribute $`n`$ objects in distinguishable boxes $`B_1,\mathrm{},B_k`$, of sizes $`n_1,\mathrm{},n_k`$, respectively . This leads to

$$\frac{1}{k!}\underset{\{n_1,\mathrm{},n_k\}^{}}{}\left(\genfrac{}{}{0pt}{}{n}{n_1,\mathrm{},n_k}\right)\underset{j=1}{\overset{k}{}}a_{n_j},$$

as the number of functional graphs with $`k`$ connected components, where $`\{n_1,\mathrm{},n_k\}^{}`$ means that the sum is over all vectors $`(n_1,\mathrm{},n_k)`$ satisfying (3.3a), and the factor $`\frac{1}{k!}`$ stands to account for the unavoidable repetitions incurred by our above assumption of ‘distinguishable boxes’.

Therefore, for $`k=1,\mathrm{},n`$, the distribution for the number of connected components can be expressed as

$$\rho _n(k)=\frac{1}{n^nk!}\underset{\{n_1,\mathrm{},n_k\}^{}}{}\left(\genfrac{}{}{0pt}{}{n}{n_1,\mathrm{},n_k}\right)\underset{j=1}{\overset{k}{}}a_{n_j}.$$
$`(3.3b)`$
It is not difficult to see that the sum in (3.3b) has as much as $`\left(\genfrac{}{}{0pt}{}{n1}{k1}\right)`$ terms , a figure growing exponentially with $`n`$ when $`k\frac{n}{2}`$. This fact dooms to failure any numerical application of (3.3b) as it stands, and hinders further analytical work.

The approaches of Folkert and Harris lead to increasing complications. The former yields

$$\rho _n(k)=\frac{1}{n^n}\underset{\mu =k}{\overset{n}{}}\frac{𝒮_\mu ^{(k)}}{\mu !}\underset{\{n_1,\mathrm{},n_\mu \}^{\prime \prime }}{}\left(\genfrac{}{}{0pt}{}{n}{n_1,\mathrm{},n_\mu }\right)n_1^{n_1}\mathrm{}n_\mu ^{n_\mu },$$
$`(3.4a)`$
where $`𝒮_\mu ^{(k)}`$ are the Stirling numbers of the first kind \[11-13\] (see Appendix), and $`\{n_1,\mathrm{},n_\mu \}^{\prime \prime }`$ means that the sum is over all vectors $`n_1,\mathrm{},n_\mu `$, constrained by

$$\underset{j=1}{\overset{\mu }{}}n_j=n\text{with}\text{ }1n_jn.$$
$`(3.4b)`$
Harris managed to propose

$$\rho _n\left(k\right)=\frac{n!}{n^nk!}\underset{\mu =k}{\overset{n}{}}𝒮_\mu ^{\left(k\right)}\underset{\{n_1,\mathrm{},n_n\}^{\prime \prime \prime }}{}\frac{1}{n_1!,\mathrm{},n_n!}\left(\frac{1^1}{1!}\right)^{n_1}\left(\frac{2^2}{2!}\right)^{n_2}\mathrm{}\left(\frac{n^n}{n!}\right)^{n_n},$$
$`(3.5a)`$
in which $`\{n_1,\mathrm{},n_n\}^{\prime \prime \prime }`$ means that the sum is over all vectors $`n_1,\mathrm{},n_n`$, constrained by

$$\underset{j=1}{\overset{n}{}}n_j=\mu \text{and}\text{ }\underset{j=1}{\overset{n}{}}jn_j=n\text{with}\text{ }0n_jn.$$
$`(3.5b)`$
Expressions (3.4) and (3.5) appear more difficult to handle than (3.3) because of the added extra terms via the summation involving Stirling numbers.

To our knowledge, the statistical moments of the distribution for the number of connected components have not yet been obtained in closed form. Further, the exact calculation done by Kruskal for the expected number of attractors (1.1) is not derived from a distribution (see Ref. ). For application purposes ($`n>>1`$) it is then worthwhile to derive manageable asymptotic formulas for $`\rho _n(k)`$; this is the subject of the next section.

4. Asymptotic Expressions

Let us start by defining

$$\beta _m=2e^m\frac{a_m}{\left(m1\right)!}$$
$`(4.1)`$
and

$$\alpha _m=1\beta _m.$$
$`(4.2)`$
Due to the asymptotic relation (3.2) it happens that

$$\beta _m=1+𝒪\left(\frac{1}{m}\right),\text{for}\text{ }m>>1,$$
$`(4.3a)`$
and

$$\alpha _m=𝒪\left(\frac{1}{m}\right),\text{for}\text{ }m>>1.$$
$`(4.3b)`$
Now, express (3.3b) in terms of $`\beta _m`$ and use (3.3a) to get

$$\rho _n\left(k\right)=\frac{n!e^n}{n^nk!2^k}\underset{\{n_1,\mathrm{},n_k\}^{}}{}\underset{j=1}{\overset{k}{}}\frac{\beta _{n_j}}{n_j}.$$

Constraint (3.3a) can be broken by introducing a Kronecker delta inside the summation. Using the integral representation

$$\delta _{n,m}=\frac{1}{2\pi i}_\mathrm{\Gamma }\frac{z^m}{z^{n+1}}𝑑z,$$

where $`\mathrm{\Gamma }`$ is any closed contour in the complex plane of $`z`$ such that the origin is inside, we go back to the expression for $`\rho _n(k)`$ to obtain

$$\rho _n\left(k\right)=\frac{n!e^n}{n^nk!2^k2\pi i}_\mathrm{\Gamma }\frac{dz}{z^{n+1}}g_k\left(z\right)$$

with a modicum of algebra, where

$$g_k(z)=\left[\phi _p\left(z\right)\right]^k\text{with}\text{ }\phi _p\left(z\right)=\underset{m=1}{\overset{p}{}}\frac{\beta _m}{m}z^m\text{and}\text{ }pn.$$
$`(4.4a)`$
Since $`g_k(z)`$ is an analytic function we can apply the Cauchy integral theorem to find

$$\rho _n(k)=\frac{e^n}{n^nk!2^k}g_k^{(n)}(0)$$
$`(4.4b)`$
as an alternate way to compute $`\rho _n(k)`$.

The computation of (4.4) may be as difficult as that of (3.3); however, the former is manageable for $`n>>1`$, as we shall see. First note that (4.4b) does not depend on the particular value we give to $`p`$ in (4.4a), as long as $`pn`$, a feature that we shall use to our convenience. Taking $`p=n`$ in (4.4a), $`g_k(z)`$ becomes a polynomial of degree $`nk`$, and its $`n`$-th derivative may be calculated in terms of finite differences with the help of Stirling numbers (see formulas ($`A`$.6) in Appendix), obtaining

$$g_k^{\left(n\right)}\left(0\right)=n!\underset{r=n}{\overset{nk}{}}\frac{𝒮_r^{\left(n\right)}}{r!}\mathrm{\Delta }^rg_k\left(0\right),$$
$`(4.5a)`$
where $`𝒮_r^{(n)}`$ are the Stirling numbers of the first kind, and

$$\mathrm{\Delta }^rg_k(0)=\underset{m=1}{\overset{r}{}}(1)^{rm}\left(\genfrac{}{}{0pt}{}{r}{m}\right)g_k(m).$$
$`(4.5b)`$
Now we provide asymptotic approximations for $`\phi _p\left(z\right)`$ in (4.4a) with $`p=n`$ and $`n>>1`$. Consider first $`z>1`$, and note that, due to (4.3),

$$\phi _n\left(z\right)=\underset{k=1}{\overset{n}{}}\frac{\beta _k}{k}z^k\underset{k=1}{\overset{n}{}}\frac{1}{k}z^k\left(1+𝒪\left(\frac{1}{n}\right)\right).$$

Now, the Euler-Maclaurin formula states

$$\begin{array}{ccc}\hfill \underset{k=1}{\overset{n}{}}f\left(k\right)& _1^nf\left(x\right)𝑑x+\frac{1}{2}\left(f\left(n\right)+f\left(1\right)\right)\hfill & \\ & +\underset{k=1}{\overset{\mathrm{}}{}}\frac{B_{2k}}{\left(2k\right)!}\left(f^{\left(2k1\right)}\left(n\right)f^{\left(2k1\right)}\left(1\right)\right),\hfill & (4.6a)\hfill \end{array}$$

where $`f\left(x\right)`$ is a $`C^{\mathrm{}}`$ function over the interval $`[1,n]`$, and $`B_k`$ are the Bernoulli numbers defined by the generating function

$$\frac{t}{e^t1}=\underset{k=0}{\overset{\mathrm{}}{}}\frac{B_k}{k!}t^k\text{for}\text{ }t<2\pi ,$$
$`(4.6b)`$
with $`B_{2k+1}=0`$, for $`k=1,2,\mathrm{}`$

Let us take $`f(x)=\frac{z^x}{x}`$ with $`z>1`$, and substitute in (4.6a) obtaining

$$\underset{k=1}{\overset{n}{}}\frac{\beta _k}{k}z^k\frac{z^n}{n}\left[\frac{1}{\mathrm{ln}z}+\frac{1}{2}+\underset{k=1}{\overset{\mathrm{}}{}}\frac{B_{2k}}{\left(2k\right)!}\left(\mathrm{ln}z\right)^{2k1}\right]\left(1+𝒪\left(\frac{1}{n}\right)\right).$$

Using (4.6b) yields

$$\underset{k=1}{\overset{n}{}}\frac{\beta _k}{k}z^k\frac{z}{z1}\frac{z^n}{n}\left(1+𝒪\left(\frac{1}{n}\right)\right)\text{for}\text{ }1<z<e^{2\pi },$$

which can be continued analytically for any $`z>1`$, thus obtaining

$$\phi _n\left(z\right)=\underset{m=1}{\overset{n}{}}\frac{\beta _m}{m}z^m\frac{z}{z1}\frac{z^n}{n}\zeta _n\text{for}\text{ }z>1,$$
$`(4.7a)`$
where

$$\zeta _n=1+𝒪\left(\frac{1}{n}\right).$$
$`(4.7b)`$
Instead, for $`z=1`$, clearly

$$\phi _n\left(1\right)=\underset{m=1}{\overset{n}{}}\frac{\beta _m}{m}\left(\mathrm{ln}n+\gamma \tau _n\right)\eta _n,$$
$`(4.7c)`$
where

$$\eta _n=1+𝒪\left(\frac{1}{n\mathrm{ln}n}\right),$$
$`(4.7d)`$
$$\tau _n=\underset{m=1}{\overset{n}{}}\frac{\alpha _m}{m},$$
$`(4.7e)`$
and $`\gamma =\mathrm{0.577\; 215\; 66}\mathrm{}`$ is the Euler-Mascheroni constant. Note that, due to (4.3b), when $`n\mathrm{}`$ in (4.7e), the limit $`\tau `$ exists. Substituting back (4.7a,b) into (4.4a), expanding then in Taylor series the term $`\frac{z}{z1}`$ for $`z>1`$, and cutting the series beyond order $`nk`$ we come to

$$g_k\left(z\right)\frac{z^{nk}}{n^k}\zeta _n^k\underset{m=0}{\overset{nk}{}}\frac{\left(k\right)_m}{m!}\frac{1}{z^m}\left(1+𝒪\left(\frac{1}{z^{nk}}\right)\right)\text{for}\text{ }z>1,$$
$`(4.8a)`$
with

$$(k)_m=k(k+1)\mathrm{}(k+m1)=\frac{\mathrm{\Gamma }(k+m)}{\mathrm{\Gamma }(k)}.$$

While substituting (4.7c-e) into (4.4a) gives us

$$g_k\left(1\right)\left(\mathrm{ln}n+\gamma \tau _n\right)^k\eta _n^k.$$
$`(4.8b)`$
Replace (4.8) into (4.5b) to get

$$\mathrm{\Delta }^rg_k\left(0\right)\left(1\right)^{r1}r\left[g_k\left(1\right)\lambda _n^{\left(k\right)}\right]+\frac{r!}{n^k}\zeta _n^k\underset{l=0}{\overset{nk}{}}\frac{\left(k\right)_l}{l!}𝒯_{nkl}^{\left(r\right)},$$
$`(4.9)`$
where we have used ($`A`$.4) to introduce the Stirling numbers of the second kind $`𝒯_p^{(l)}`$, and

$$\lambda _n^{(k)}\frac{\zeta _n^k}{n^k}\underset{l=0}{\overset{nk}{}}\frac{(k)_l}{l!}.$$

The summation for $`\lambda _n^{(k)}`$ can be found in summation tables yielding

$$\underset{l=0}{\overset{nk}{}}\frac{(k)_l}{l!}=\frac{\left[k\left(n+1\right)\right]!}{k!(nk)!}.$$

Next, applying Stirling’s approximation ($`s!\sqrt{2\pi s}s^se^s`$ for $`s>>1`$) to the factorials for $`n>>1`$, one finds $`\lambda e^k`$, while $`g_k(1)e^{k\left\{\mathrm{ln}[\mathrm{ln}n+\gamma \tau _n]+\mathrm{ln}\eta _n\right\}}`$; thus, $`\lambda `$ is exponentially small as compared to $`g_k\left(1\right)`$, and can be neglected in (4.9). Replacing back (4.9) into (4.5a), using properties ($`A`$.5), and substituting in (4.4b), we find for the asymptotic distribution of the number of attractors

$$\rho _n\left(k\right)\frac{\sqrt{2\pi n}}{k!2^k}g_k\left(1\right)\underset{l=n}{\overset{nk}{}}\frac{\left(1\right)^{l1}}{\left(l1\right)!}𝒮_k^{\left(n\right)}+𝒩_n\left(k\right),$$
$`(4.10)`$
where

$$𝒩_n\left(k\right)=\frac{\sqrt{2\pi n}}{k!2^k}\frac{1}{n^k}\zeta _n^k\frac{\left(k\right)_{\left[n\left(k1\right)\right]}}{\left[n\left(k1\right)\right]!}.$$
$`(4.11)`$
We show now that $`𝒩_n(k)`$ is a null term, that is to say

$$\underset{n\mathrm{}}{lim}\underset{k=1}{\overset{n}{}}k^m𝒩_n(k)=0$$
$`(4.12)`$
for any $`m0`$, and therefore it does not contribute to the averages taken with $`\rho _n(k)`$, for $`n>>1`$. First note that $`𝒩_n\left(k\right)\frac{1}{\sqrt{n}}`$ for any $`k`$, so as in the summation (4.12) the terms with $`k𝒪(1)`$ do not contribute. Hence, we can make the approximation for $`n>>1`$ and $`k>>1`$ in (4.11), and calculate (4.12) with it, obtaining

$$𝒩_n\left(k\right)\frac{e^{\left(1+\frac{1}{2n}\right)}}{\sqrt{n}}\frac{1}{k!}\left[\frac{1}{2}\zeta _ne^{1+\frac{1}{2n}}\right]^k.$$

Since

$$\underset{k=1}{\overset{\mathrm{}}{}}k^m\frac{x^k}{k!}=𝒫_m(x)e^x,$$

where $`𝒫_m(x)`$ is a polynomial in $`x`$ of degree $`m`$, we have

$$\underset{k=1}{\overset{n}{}}k^m𝒩_n\left(k\right)\frac{e^{\left(1+\frac{1}{2n}\right)}}{\sqrt{n}}𝒫_m\left(\frac{1}{2}\zeta _ne^{1+\frac{1}{2n}}\right)\mathrm{exp}\left(\frac{1}{2}\zeta _ne^{1+\frac{1}{2n}}\right),$$

and thus the expression goes to zero as $`\frac{1}{\sqrt{n}}`$, which proves (4.12).

Then, (4.10) can be expressed very simply in terms of Stirling numbers of the first kind by noting the following: let us take for a moment $`\alpha _m=0`$ in (4.4a), and also set $`p=\mathrm{}`$ with $`z<1`$ (so the infinite sum converges), obtaining

$$\stackrel{~}{g}_k\left(z\right)=\left(\underset{m=1}{\overset{\mathrm{}}{}}\frac{1}{m}z^m\right)^k=\left\{\mathrm{ln}\left(1z\right)\right\}^k,$$

where the tilde indicates that we have set $`\alpha _m=0`$. Using ($`A`$.2) and substituting in (4.4b) we arrive to

$$\stackrel{~}{\rho }_n\left(k\right)=\frac{e^n}{n^n2^k}\left(1\right)^{nk}𝒮_n^{\left(k\right)}.$$

Comparing with (4.10), we obtain, up to the null term, that

$$\rho _n\left(k\right)\frac{e^n}{n^n}\left(1\right)^{nk}𝒮_n^{\left(k\right)}\left(\frac{\mu _n}{2}\right)^k,$$
$`(4.13a)`$
where

$$\mu _n=1\frac{\tau _n}{\gamma +\mathrm{ln}n}.$$
$`(4.13b)`$
Since $`\rho _n\left(k\right)`$ is a statistical distribution,

$$\underset{n\mathrm{}}{lim}\underset{k=1}{\overset{n}{}}\rho _n(k)=1.$$

Which allows us to compute $`\tau =lim_n\mathrm{}\tau _n`$ in (4.13b), by using property ($`A`$.1) of Stirling numbers, giving

$$\tau =\mathrm{ln}2.$$
$`(4.13c)`$
We can replace back (4.13c) into (4.13a) and (4.13b), yielding and error of order $`𝒪\left(1/\mathrm{ln}n\right)`$, for $`1kn`$. So we have, for the asymptotic approximation of the distribution of the number of attractors, the following expression:

$$\rho _n\left(k\right)\frac{e^n}{n^n}\left(1\right)^{nk}𝒮_n^{\left(k\right)}\left(\frac{\mu }{2}\right)^k\left(1+𝒪\left(\frac{1}{\mathrm{ln}n}\right)\right),$$
$`(4.14a)`$
where

$$\mu =1\frac{\mathrm{ln}2}{\mathrm{ln}n}.$$
$`(4.14b)`$
Now we can use (4.14) to calculate averages. Using expression ($`A`$.3), we have

$$\mathrm{\Theta }\left(n\right)\frac{1}{2}\left(\mathrm{ln}2n+\gamma +𝒪\left(\frac{1}{\mathrm{ln}n}\right)\right)$$
$`(4.15)`$
for the average number of attractors $`\mathrm{\Theta }\left(n\right)`$, which coincides, as it should be, with Kruskal asymptotic approximation (1.2). The variance can also be computed by using the derivative of ($`A`$.3):

$$\sigma ^2\mathrm{\Theta }\left(n\right)\left(1+𝒪\left(\frac{1}{\mathrm{ln}n}\right)\right),$$
$`(4.16)`$
which is a new result.

Conclusion

We have proposed a new and simpler expression for the distribution of the number of attractors in the random map model (Eq. (3.3)). The number of operations involved for the numerical evaluation of the distribution grows exponentially with $`n`$, thus making expression (3.3) useless for direct calculations for $`n>>1`$. To overcome this difficulty, however, we have derived an asymptotic formula (Eq. (4.14)), from which we directly deduced asymptotic values both for the average number of attractors and for its variance (Eqs. (4.15) and (4.16), resp.).

In the random map model, additional statistical figures are of interest. Among them, the average attractor size, the average period (or average length) of the cycles; also, given a point $`x`$ in the phase space, the expected length of its orbit, the expected length of the cycle in the attractor containing $`x`$, the expected number of points from which $`x`$ can be attained, etc. Some of these figures have already been computed (for example, see Ref. ), and others are subject of recent studies .

Acknowledgments

This work is supported in part by CONACyT project number U40004-F. The authors would like to thank V. Domínguez, E. Sacristán, F. Toledo for computational advice. The second author (FZ) thanks Alberto Verjovsky for his invaluable clarification of some aspects of asymptotic series.

Appendix: Properties of Stirling Numbers

Here we list the equations and definitions on Stirling numbers necessary to follow the calculations in this article. For a more extensive treatment see Refs. \[11-13\]. Stirling numbers of the first kind are generated by the functions

$$z(z+1)\mathrm{}(z+n1)=\underset{m=1}{\overset{n}{}}(1)^{nm}𝒮_n^{(m)}z^m$$
$`(A.1)`$
and

$$\left\{\mathrm{ln}\left(1+z\right)\right\}^k=k!\underset{r=k}{\overset{\mathrm{}}{}}\frac{𝒮_r^{(k)}}{r!}z^r\text{for}\text{ }z<1.$$
$`(A.2)`$
Applying the operator $`z\frac{d}{dz}`$ to equation ($`A`$.1) we get the important relation

$$z\left(z+1\right)\mathrm{}\left(z+n1\right)\underset{m=1}{\overset{n}{}}\frac{z}{z+m1}=\underset{m=1}{\overset{n}{}}\left(1\right)^{nm}𝒮_n^{\left(m\right)}mz^m.$$
$`(A.3)`$
Stirling numbers of the first kind may be expressed in closed form by

$$𝒮_k^{\left(n\right)}=\underset{l=0}{\overset{kn}{}}\left(1\right)^l\left(\genfrac{}{}{0pt}{}{k1+l}{kn+l}\right)\left(\genfrac{}{}{0pt}{}{2kn}{knl}\right)𝒯_{kn+l}^{\left(l\right)},$$

where

$$𝒯_p^{\left(l\right)}=\frac{1}{l!}\underset{k=0}{\overset{l}{}}\left(1\right)^{lk}\left(\genfrac{}{}{0pt}{}{l}{k}\right)k^p$$
$`(A.4)`$
are the Stirling numbers of the second kind. It is known that Stirling numbers verify the relations

$$𝒮_k^{\left(n\right)}=0\text{if}\text{ }k<n,$$
$`(A.5a)`$
$$𝒯_p^{\left(l\right)}=0\text{if}\text{ }p<l,$$
$`(A.5b)`$
and

$$\underset{k=m}{\overset{n}{}}𝒮_k^{\left(m\right)}𝒯_n^{\left(k\right)}=\underset{k=m}{\overset{n}{}}𝒮_n^{\left(k\right)}𝒯_k^{\left(m\right)}=\delta _{m,n}.$$
$`(A.5c)`$
By means of Stirling numbers of the first kind it is possible to express derivatives of a function in terms of finite differences by the formula

$$\frac{d^m}{dz^m}f\left(z\right)=m!\underset{k=m}{\overset{\mathrm{}}{}}\frac{𝒮_k^{\left(m\right)}}{k!}\mathrm{\Delta }^kf\left(z\right),$$
$`(A.6a)`$
if the summation is convergent, and where $`\mathrm{\Delta }f\left(x\right)=f\left(x+1\right)f\left(x\right)`$ and

$$\mathrm{\Delta }^kf\left(z\right)=\underset{l=0}{\overset{k}{}}\left(1\right)^{kl}\left(\genfrac{}{}{0pt}{}{k}{l}\right)f\left(z+l\right).$$
$`(A.6b)`$
References

 Kauffman, S.A., Metabolic stability and epigenesis in randomly connected nets, J. Theoret. Biol. 22 (1969) 437; Cellular homeostasis, epigenesis and replication in randomly aggregated macromolecular systems, J. Cybernetics 1 (1971) 71; Wolfram, S., Statistical mechanics of cellular automata, Rev. Mod. Phys. 55 (1983) 601; Guanajuato Lectures on Complex Systems and Binary Networks. Springer Verlag Lecture Notes series. Eds. R. López Peña, R. Capovilla, R. García-Pelayo, H. Waelbroeck and F. Zertuche. (1995); Aldana M., Coppersmith S. and Kadanoff L. Boolean Dynamics with Random Couplings. nlin.AO/020406 (2002).

 Coste, J. and Henon, M., In Disordered systems and biological organization, Eds. Bienenstock, M.Y., Fogelman Soulié F., and Weisbuch, G., p. 361, Springer Verlag, Heidelberg (1986).

 Derrida, B. and Flyvbjerg, H., The random map model: a disordered model with deterministic dynamics, J. Physique 48 (1987) 971.

 Kauffman, S.A., The origins of order: Self-organization and selection in evolution. Oxford University Press (1993).

 Metropolis, N. and Ulam, S., A property of randomness of an arithmetical function, Amer. Math. Monthly 60 (1953) 252.

 Kruskal, M.D., The expected number of components under a random mapping function, Am. Math. Monthly 61 (1954) 392.

 Rubin, H. and Sitgreaves, R., Probability distributions related to random transformations on a finite set, Tech. Rept. No. 19A Applied Mathematics and Statistics Laboratory, Stanford University (1954). Unpublished.

 Folkert, J.E., The distribution of the number of components of a random mapping function, (1955) Unpublished Ph. D. dissertation, Michigan State University.

 Katz, L., Probability of indecomposability of a random mapping function, Ann. Math. Stat. 26 (1955) 512.

 Harris, B., Probability distributions related to random mappings, Ann. Math. Stat. 31 (1960) 1045.

 Abramowitz, M. and Stegun, I.A., Handbook of mathematical functions. Dover Publications, New York (1972).

 Jordan, C., Calculus of finite differences. Chelsea Publishing Company, New York (1947).

 Graham, R.L., Knuth, D.E. and Patashnik, O., Concrete Mathematics. Addison-Wesley Publishing Company, New York (1994).

 Arfken, G., Mathematical Methods for Physicists. Academic Press, New York (1970) Chap. 10.

 Mangulis, V., Handbook of Series for Scientists and Engineers. Academic Press. New York and London (1965). Page 60.

 Romero, D. and Zertuche, F., Grasping the Connectivity of Functional Graphs, submitted for publication (2002).


text:
# Gravitational wave detection using electromagnetic modes in a resonance cavity

##

During the last decades the quest for detecting gravitational waves has intensified. The efforts have been inspired by the indirect evidence for gravitational radiation , advances in technology and the prospects of obtaining new useful astrophysical information through the development of gravitational wave astronomy . A number of ambitious detector projects are already in operation or being built all over the world, for example Ligo and ALLEGRO in USA, Virgo, AURIGA and GEO 600 in Europe, Tama 300 in Japan, AIGO and NIOBE in Australia . Furthermore, there are well developed plans to use space based gravitational wave detectors, i.e., the LISA project . The detection mechanisms are basically of a mechanical nature in the cases above, but there have also been several proposals for electromagnetic detection mechanisms .

In the present paper we will investigate a detection mechanism based on the interaction of electromagnetic modes and gravitational radiation in a cavity with highly conducting walls. The main feature of our proposed gravitational wave detector is that it supports two electromagnetic eigenmodes with nearby eigenfrequencies, a possibility that has previously been discussed in Refs. . If one eigenmode is excited initially (called the pump-mode), and a quasi-monochromatic gravitational wave with a frequency equal to the eigenmode frequency difference reaches our system, a new electromagnetic eigenmode can be excited due to the gravitational-electromagnetic wave interaction. The coupling mechanism is similar in principle to the wave interaction processes described in, e.g., Ref. .

The low-frequency nature of the gravitational modes, as compared to the electromagnetic resonance frequencies, at first seem to greatly limit the efficiency of a cavity with nearby frequencies. To get a large gravitationally induced mode-coupling for a simple cavity geometry, the estimated cavity dimensions become prohibitively large, i.e., comparable to the wavelength of the gravitational wave. Solutions to this problem has been found by Refs. , who has considered a gravitational wave detector consisting of two coupled cavities. Cavities based on these principles have been built, and experimental results are presented in Refs. . In this letter we consider a single cavity with a variable crossection. The main purpose of varying the crossection is the following. For a cavity with dimensions much smaller than the wavelength of the gravitational wave, the suggested geometry greatly magnifies the gravitational wave induced mode-coupling. In our present work we have simulated the effect of a varying crossection, by considering a cavity filled with three different dielectrics, in order to be able to perform most of the calculations analytically. It is straightforward to make a semi-quantitative translation of our results to the case of a vacuum cavity with a varying crossection. Using current data on the latter type of cavity , we estimate the minimum detection level of the metric perturbation to the value $`h_{\mathrm{min}}7\times 10^{23}`$, where we have considered an inspiraling neutron star or black hole pair as a gravitational wave source. If such a level of sensitivity can be reached, neutron star or black hole binaries close to collapse could be detected at a distance $`r_{\mathrm{max}}10^8\mathrm{ly}`$. Adopting data from Ref.  for the occurrence of compact binary mergers, we obtain the estimate 1–20 detection events per year.

In vacuum, a linearized gravitational wave can be represented by $`g_{\mu \nu }=\eta _{\mu \nu }+h_{\mu \nu }(\xi )`$ where $`\xi xct`$, and $`h_{xx}=h_{yy}h_+,h_{xy}=h_{yx}h_\times `$ in standard notation.

Neglecting terms proportional to derivatives of $`h_+`$ and $`h_\times `$, the wave equation for the magnetic field is

$$\left[\frac{n^2}{c^2}\frac{^2}{t^2}^2\right]𝑩=\left[h_+\left(\frac{^2}{y^2}\frac{^2}{z^2}\right)+h_\times \frac{^2}{yz}\right]𝑩,$$
(1)
and similarly for the electric field. Here $`n`$ is the index of refraction. For the moment, we will neglect mechanical effects, i.e., effects which are associated with the varying coordinates of the walls due to the restoring forces of the cavity.111This is true neglecting thermal noise and assuming a long cavity, as compared to the speed of sound over the gravitational wave frequency.

The coupling of two electromagnetic modes and a gravitational wave in a cavity will depend strongly on the geometry of the electromagnetic eigenfunctions. We can greatly magnify the coupling, as compared to a rectangular prism geometry, by varying the cross-section of the cavity, or by filling the cavity partially with a dielectric medium. The former case is of more interest from a practical point of view, since a vacuum cavity implies better detector performance, but we will consider the latter case since it can be handled analytically. *However*, we will show how to make a semi-quantitative translation of our results to the case of a varying cavity cross-section.

Specifically, we choose a rectangular cross-section (side lengths $`x_0`$ and $`y_0`$), and we divide the length of the cavity into three regions. Region 1 has length $`l_1`$ (occupying the region $`0<z<l_1`$) and a refractive index $`n_1`$. Region 2 has length $`l_2`$ (occupying the region $`l_1<z<l_1+l_2`$), with a refractive index $`n_2`$, while region 3 consists of vacuum and has length $`l_3`$ (occupying the region $`l_1+l_2<z<l_1+l_2+l_3l`$). The cavity is supposed to have positive coordinates, with one of the corners coinciding with the origin. Furthermore, we require that $`l_2l_1`$, and that the wave number in region 2 is less than in region 1. The reason for this arrangement is twofold. Firstly, we want to obtain a large coupling between the wave modes, and secondly we want an efficient filtering of the eigenmode with the lower frequency in region three.

The first step is to analyze the linear eigenmodes in this system. The simplest modes are of the type

$`E_y={\displaystyle \frac{i\omega x_0}{m\pi }}\stackrel{~}{B}_{zj}\mathrm{sin}\left({\displaystyle \frac{m\pi x}{x_0}}\right)\mathrm{sin}[k_jz+\phi _j]e^{i\omega t},`$ (2a)
$`B_z=\stackrel{~}{B}_{zj}\mathrm{cos}\left({\displaystyle \frac{m\pi x}{x_0}}\right)\mathrm{sin}[k_jz+\phi _j]e^{i\omega t},`$ (2b)
$`B_x={\displaystyle \frac{k_jx_0}{m\pi }}\stackrel{~}{B}_{zj}\mathrm{sin}\left({\displaystyle \frac{m\pi x}{x_0}}\right)\mathrm{cos}[k_jz+\phi _j]e^{i\omega t},`$ (2c)
in regions $`j=1`$, $`2`$ and $`3`$, respectively, where the wave in region 3 is a standing wave, and $`m`$ is the mode number. In region 3 we may also have a decaying wave

$`E_y={\displaystyle \frac{i\omega x_0}{m\pi }}\stackrel{~}{B}_{z3}\mathrm{sin}\left({\displaystyle \frac{m\pi x}{x_0}}\right)\mathrm{sinh}[k_3z+\phi _3]e^{i\omega t},`$ (2ca)
$`B_z=\stackrel{~}{B}_{z3}\mathrm{cos}\left({\displaystyle \frac{m\pi x}{x_0}}\right)\mathrm{sinh}[k_3z+\phi _3]e^{i\omega t},`$ (2cb)
$`B_x={\displaystyle \frac{k_3x_0}{m\pi }}\stackrel{~}{B}_{z3}\mathrm{sin}\left({\displaystyle \frac{m\pi x}{x_0}}\right)\mathrm{cosh}[k_3z+\phi _3]e^{i\omega t}.`$ (2cc)
Using standard boundary conditions, the wave numbers are calculated for an eigenmode, and the relation between the amplitudes in the three regions is found, and thereby the mode profile. We are interested in the shift from decaying to oscillatory behavior in region 3. Denote the highest frequency mode which is decaying in region 3 with index $`a`$, and the wave number and decay coefficient with $`k_{1a}`$, $`k_{2a}`$ and $`k_{3a}`$ respectively. Similarly, the next frequency, which is oscillatory in both regions, is denoted by index $`b`$. If we have $`lx_0`$ (and $`m`$ the same) these two frequencies will be very close, and a gravitational wave which has a frequency equal to the difference between the electromagnetic modes causes a small coupling between these modes. An example of two such eigenmodes is shown in fig. 1.

We define the eigenmodes to have the form $`\stackrel{~}{B}_{za,b}=\stackrel{~}{B}_{a,b}(t)B_{\mathrm{norm}a,b}(z,y,z)`$, where $`\stackrel{~}{B}_{a,b}`$ is a time-dependent amplitude and the normalized eigenmodes $`B_{\mathrm{norm}a,b}`$ satisfy $`_V\left|B_{\mathrm{norm}a,b}\right|^2𝑑V=V`$. We let all electromagnetic field components be of the form $`A=A_a(𝒓)\mathrm{exp}(i\omega _at)+A_b(𝒓)\mathrm{exp}(i\omega _bt)+\mathrm{c}.\mathrm{c}.`$, where c.c. stands for complex conjugate, and the indices stand for the eigenmodes discussed above. The gravitational perturbation can be approximated by $`h_{+,\times }=\widehat{h}_{+,\times }\mathrm{exp}(i\omega _gt)+\mathrm{c}.\mathrm{c}`$, where we neglect the spatial dependency, since the gravitational wavelength is assumed to be much longer than all of the cavity dimensions. During a certain interval in time, the frequency matching condition $`\omega _b=\omega _a+\omega _g`$ will be approximately fulfilled. Given the wave equation (1), and the above ansatz we find after integrating over the length of the cavity

$$\frac{2i\omega _b}{c^2}\left(\frac{}{t}\gamma \right)\stackrel{~}{B}_b=h_+k_{1a}^2\stackrel{~}{B}_aI_c,$$
(2cd)
where

$$I_c=\frac{1}{Vk_{1a}^2}_V\frac{^2B_{\mathrm{norm}a}}{z^2}B_{\mathrm{norm}b}𝑑V,$$
(2ce)
and we have added a phenomenological linear damping term represented by $`\gamma `$. Thus we note that for the given geometry, only the $`h_+`$-polarization gives a mode-coupling. Calculations of the eigenmode parameters show that $`I_c`$ may be different from zero when $`n_{1,2}1`$, and generally $`I_c`$ of the order of unity can be obtained, see fig. 1 for an example. From Eq. (2cd), we find that the saturated value of the gravitationally excited mode is

$$\left|\stackrel{~}{B}_{b\mathrm{sat}}\right|=\frac{h_+k_{1a}^2c^2\stackrel{~}{B}_a}{2\gamma \omega _b}I_c.$$
(2cf)
In fig. 1 it is shown that we can get an appreciable mode-coupling constant $`I_c`$ for a cavity filled with materials with different dielectric constants, and it is of interest whether or not this can be achieved in a vacuum cavity. As seen by Eq. (2ce), the coupling is essentially determined by the wave numbers of the modes, given by $`k=(m^2\pi ^2/x_0^2n^2\omega ^2/c^2)^{1/2}`$. Thus by adjusting the width $`x_0`$ in a vacuum cavity, we may get the same variations in the wave numbers as when varying the index of refraction $`n`$. The translation of our results to a vacuum cavity with a varying width is not completely accurate, however. When varying $`x_0`$, the mode-dependence on $`x`$ and $`z`$ does not exactly factorize, in particular close to the change in width. Moreover, the contribution to the coupling $`I_c`$ in each section becomes proportional to the corresponding volume, and thereby also to the cross-section. However, since most of the contribution to the integral in Eq. (2ce) comes from region 1, our results can still be approximately translated to the case of a vacuum cavity, by varying $`x_0`$ instead of $`n`$ such as to get the same wavenumber as in our above example.

We denote the minimum detection level of the excited mode with $`\stackrel{~}{B}_{b\mathrm{min}}`$ and the maximum allowed field in the cavity with $`\stackrel{~}{B}_{a\mathrm{max}}`$. The magnetic field $`\stackrel{~}{B}_{b\mathrm{min}}`$ is related to the minimum number of photons $`N_{bm}`$ needed for detection by $`\stackrel{~}{B}_{b\mathrm{min}}\left(N_{bm}\mu _0\omega _b\mathrm{}/x_0y_0l\right)^{1/2}`$. Furthermore, we have $`2\gamma \omega _b/k_{1a}^2c^2`$ $`(1/3Q)(\omega _b^2/k_{1a}^2c^2)2/3Q`$, where $`Q=\omega _b/2\pi \gamma `$ is the quality factor of the cavity, and thus we obtain, using Eq. (2cf),

$$h_{\mathrm{min}}\frac{2}{3QI_c\stackrel{~}{B}_{a\mathrm{max}}}\left(\frac{\mu _0\omega _b\mathrm{}N_{bm}}{x_0y_0l}\right)^{1/2}.$$
(2cg)
Before giving estimates of the parameter values, we will investigate certain other effects that may limit the detection efficiency.

Massive binaries produce monochromatic radiation to a good approximation, but close to merging, the frequency will be increasing rather rapidly which means that the phase matching will be lost. The cavity is designed to detect the frequency $`\omega _g`$, and we have assumed that the signal varies as $`\mathrm{exp}(\mathrm{i}\omega _gt)`$, but in reality we have $`\mathrm{exp}(\mathrm{i}\psi (t))`$, where for simplicity we assume $`\psi =0`$, $`d\psi /dt=\omega _g`$ at $`t=0`$. The coherence time $`t_{\mathrm{coh}}`$ is roughly defined by $`\omega _gt_{\mathrm{coh}}2\psi (t_{\mathrm{coh}}/2)=\pi `$. Thus, $`t_{\mathrm{coh}}[2\pi /(d\omega _g/dt)]^{1/2}`$ provided $`\omega _gt_{\mathrm{coh}}1`$. Using Newtonian calculations of two masses $`M`$ in circular orbits around the center of mass, complemented by the quadrupole formula for gravitational radiation, we find

$$\frac{d\omega _g}{dt}=\frac{388\omega _g^2}{5c^5}\left(\frac{GM\omega _g}{8}\right)^{5/3}.$$
(2ch)
Close to binary merging, the time of coherent interaction will be shorter than the photon life time, and in that regime the growth of the excited mode is limited by decoherence rather than the damping due to a finite quality factor of the cavity. However, formula (2cg) can still be applied if we simply replace $`Q`$ by $`\omega _bt_{\mathrm{coh}}`$.

To be able to estimate the number of photons needed for detection, we must study various sources of noise. The simplest effect is direct thermal excitation of photons in mode $`b`$. Since each mode has a thermal energy level of order $`KT`$, the number of such photons is $`KT/\mathrm{}\omega _b`$. However, we will also have a contribution associated with thermal variations in the length of the cavity. A standard model for the variations in length $`\delta l`$ is

$$\frac{d^2\delta l}{dt^2}+\frac{\omega _l}{Q_{mec}}\frac{d\delta l}{dt}+\omega _l^2\delta l=l\omega _g^2h_+\mathrm{exp}(i\omega _gt)+a_{th},$$
(2ci)
where $`\omega _l`$ is the eigenfrequency of longitudinal oscillations (of the order of the acoustic velocity in the cavity divided by the length), $`Q_{mec}`$ is the mechanical quality factor associated with these oscillations, and $`a_{th}`$ is the stochastic acceleration due to the thermal motion, giving rise to a random walk in the oscillation amplitude. First we note that the amplitude of the length variations with the gravitational frequency is $`\delta l_g=\omega _g^2h_+/(\omega _g^2\omega _l^2)`$. We will consider the case when $`\omega _l^22\omega _g^2`$, which implies that the amplitude of the gravitational length perturbation is essentially unaffected by the restoring force of the cavity, as is the number of gravitationally generated photons.

The thermal fluctuation contribution to the right hand side of Eq. (2cd), via the coupling to the pump wave, becomes proportional to $`\delta l_{th}\mathrm{exp}[i(\omega _b\omega _a)t]`$. Since $`\delta l_{th}\mathrm{exp}(i\omega _lt)`$, the longitudinal oscillations give rise to slightly off-resonant (i.e., driven) fields with a frequency difference $`\delta \omega =\omega _l\omega _g`$ compared to mode $`b`$. However, due to the stochastic changes in amplitude and/or phase of the longitudinal oscillation (where the time-scale for significant changes is given by $`\tau Q_{mec}/\omega _l`$), a contribution to mode $`b`$ of the order $`\delta l_{th}\sqrt{\omega _l/\delta \omega Q_{mec}}dt\sqrt{KT/m\omega _l^2}\sqrt{1/Q_{mec}}(2\pi /\delta \omega )`$ is also made during a single oscillation period $`dt=2\pi /\delta \omega `$. Here $`m`$ is the mass taking part in the longitudinal oscillation. During a time of the order $`t_{coh}`$, this contribution adds up to approximately $`(\sqrt{2\pi }/\delta \omega )\sqrt{KT/m\omega _l^2}\sqrt{\delta \omega t_{coh}/Q_{mec}}`$ by a random walk process. Using $`\delta \omega \omega _g`$, the condition for the gravitational contribution $`hlt_{coh}`$ to be larger than that of the thermal fluctuations can be written

$$h_{min}>\frac{2\pi }{\omega _ll}\sqrt{\frac{KT}{m}}\frac{1}{\sqrt{Q_{mec}\omega _gt_{coh}}}.$$
(2cj)
Assume that we want to reach a sensitivity $`h_{min}7\times 10^{23}`$. We let $`x_0=0.1\mathrm{m}`$, $`y_0=0.2\mathrm{m}`$, $`l=4\mathrm{m}`$ , and $`m=4`$ which gives $`\omega _b4\pi c/x_0=4\times 10^{10}\mathrm{rad}/\mathrm{s}`$. Furthermore, we take $`\omega _g=1000\mathrm{rad}/\mathrm{s}`$ and let $`\omega _l=1500\mathrm{rad}/\mathrm{s}`$. We assume $`T=1\mathrm{K}`$ together with $`m=2500\mathrm{kg}`$. From (2cj), we then find that we need a mechanical quality factor $`Q_{mec}6\times 10^7`$ to reach the desired sensitivity.222Note that mechanical quality factors can be several orders of magnitude better, see e.g., Ref. . Moreover, assuming the necessary number of photons for detection to be $`N_{bm}=20`$ (well above the direct electromagnetic noise level of a few photons), $`I_c0.3`$ and $`\stackrel{~}{B}_{a\mathrm{max}}=0.2\mathrm{T}`$ , we need an electromagnetic quality factor $`Q=5\times 10^9`$ to reach the desired sensitivity $`h_{\mathrm{min}}7\times 10^{23}`$. Note that the coherence time is slightly longer than the photon life-time, as needed.

Following the example given above for calculating the coherence time, i.e., a binary consisting of two compact objects, each of one solar mass $`M_{}`$, separated by a distance $`d`$, the amplitude of the metric perturbation at a distance $`r`$ is given by $`h2G^2M^2/c^2rd`$. For definiteness we choose $`d=80\mathrm{km}`$ corresponding to $`\omega _g=1000\mathrm{rad}/\mathrm{s}`$, and thus for $`h_{\mathrm{min}}`$ of the order of $`7\times 10^{23}`$ we obtain the maximum observational distance $`r_{\mathrm{max}}8\times 10^7\mathrm{ly}`$. Using data from Ref. , we deduce that the number of galaxies within the observational distance $`r_{\mathrm{max}}10^8\mathrm{ly}`$ is of the order $`N=10^5`$. Combining these figures with the expected number of compact binary mergers per galaxy and million years , we obtain a detection rate in the interval $`120`$ events per year (the uncertainty is due to different models used for the birth of compact binaries).

Our proposal for a gravitational wave detector is to a large extent based on currently available technology , and our requirements are moderate given the performance of some existing microwave cavities. For example, values of $`Q4\times 10^{10}`$ has been reported in Ref.  and quantum non-demolition measurements of single microwave photons have been made in Ref. . Furthermore, the key performance parameters of the superconducting niobium cavities, i.e., the quality factor and the maximum allowed field strength before field emission, have been improving over the years, suggesting that the detection sensitivity can be increased even further.

A sensitivity $`h_{\mathrm{min}}10^{22}`$ seems extremely good, but, on the other hand, idealizations have been made when making the estimate. In addition to any effects induced by the gravitational wave, the walls generally vibrate slightly due to the electromagnetic forces exerted by the large pump field. While these later oscillations clearly will be larger than the variations in length that are directly due to the gravitational wave, the associated nonlinearities will be harmonics of the pump frequency, and thus such effects do not couple to the other eigenmodes of the cavity. Furthermore, we have assumed that the detection of the excited mode is not much affected by the presence of the pump signals. Even though mode $`a`$ is partially filtered out in region 3, the small frequency shift between the two electromagnetic modes may pose a certain difficulty in this respect: In particular, a very narrow bandwidth (pump) antenna signal must be used, in order to exclude the slightest initial perturbation at the gravitationally excited frequency $`\omega _b`$.

However, although there are technical difficulties in constructing an electromagnetic detector, the real advantage is the possibility to reduce the size of the devise. In the example presented above, the length of the cavity has been taken to be $`4\mathrm{m}`$, and the cross section roughly $`1/50\mathrm{m}^2`$. This alone could prove to be useful when trying to set up new gravitational wave observatories.

We thank Ulf Jordan for helpful discussions.

## References