format: text

text:
# 1 Introduction

## 1 Introduction

The proposal to extract the values of light-quark current masses from the QCD sum rules for the correlators of quark scalar and pseudoscalar currents is rather old (see e.g. the works of Refs.  and the reviews of Refs. ). However, in view of the appearance of new experimental data and new theoretical results and ideas, the investigations in this direction are continuing (for a summary of the modern situation, see forinstance the reviews of Refs. ).

When studies began, it was not clear what normalization point should be associated with the definite values of the light-quark current masses. Starting from the mid-90,s it became common wisdom, accepted by the Particle Data Group, to consider as phenomenolgical parameters the running quark masses of the $`\overline{\mathrm{MS}}`$ scheme, normalized at either 1 GeV or 2 GeV. Another quantity, probably more suitable for the calculations on the lattice, is the renormalization-scheme-invariant mass $`\widehat{m}`$, connected with the running one by the following solution of the renormalization group equation

$`\overline{m}(Q^2)`$ $`=`$ $`\widehat{m}\mathrm{exp}\left[{\displaystyle _0^{\alpha _\mathrm{s}(\mathrm{Q}^2)}}{\displaystyle \frac{\gamma _\mathrm{m}(\mathrm{x})}{\beta (\mathrm{x})}}\mathrm{dx}{\displaystyle _0^{2\beta _0}}{\displaystyle \frac{\gamma _0}{\beta _0\mathrm{x}}}\mathrm{dx}+2{\displaystyle \frac{\gamma _0}{\beta _0}}\mathrm{ln}(2\beta _0)\right]`$
$`=`$ $`\widehat{m}(2\beta _0\alpha _s(Q^2))^{\gamma _0/\beta _0}[1+\mathrm{}.],`$
where the perturbative QCD $`\beta `$-function and the mass-anomalous dimension function $`\gamma _m`$ are now known at the 4-loop approximations (see Ref.  and Refs.  respectively) with the one-loop scheme-independent terms $`\beta _0=[11(2/3)N_f]/4`$ and $`\gamma _0=1`$.

Another important definition of the current masses is the on-shell or pole quark masses. Despite the fact that, contrary to the case of charged leptons, quarks do not exist as free particles, this definition is commonly accepted and widely used when heavy quarks are considered. Within perturbation theory and the $`\overline{\mathrm{MS}}`$ scheme, the running quark masses and pole quark masses can be straightforwardly connected. At the next-to-leading order (NLO) this relation was calculated in Ref.  and confirmed later on in Ref. . Beyond the NLO level, higher order terms to this relation were estimated in the large-$`\beta _0`$ approximation , which can be obtained from the gauge-independent large-$`N_f`$ expansion after renormalon-inspired substitution $`N_f6\beta _0`$. At the NNLO level the correction to the relation between pole and running quark masses was estimated in Ref.  within the framework of the scheme-invariant approach developed in the works of Ref. .

The results of the subsequent approximate NNLO calculations, made by evaluating expansions of the 3-loop quark self-energy for small and large external momentum and further Padé improvements of the corresponding mass-dependent series , turned out to be in reasonable agreement with estimates, obtained by both large-$`\beta _0`$ resummation  and scheme-invariant methods . In its turn, the numerical result of the approximate calculation of Ref.  was confirmed by the recent explicit analytical calculations, performed in Ref. .

In this report the survey of the uncertainties of recent QCD sum-rule determinations of light-quark masses from the correlators of scalar and pseudoscalar currents at 3-loop  and 4-loop  perturbative QCD approximations is given. Indeed, a single-chain renormalon structure of the Borel images of the functions, related to the correlator of scalar quark currents was recently studied in collaboration with Broadhurst and Maxwell in Ref. . In this note special attention is paid to the results obtained in Ref. . The higher-order perturbative QCD uncertainties for the studied functions are estimated with the help of the large-$`\beta _0`$ approximation. The bridge between the analysis of similar questions in the case of the relation between pole and running quark masses, which contain the renormalon singularity , will also be constructed. The reason for the limitations of the scheme-invariant estimates of Ref.  in the scalar channel will be clarified.

## 2 Renormalon calculus

The studies of renormalon effects (for reviews see, ) usually start from the consideration of the contributions of the internal single chain of fermion loops to some physical quantity $`D`$, say to the anomalous magnetic moment of the electron  or to the $`e^+e^{}`$-annihilation Adler function . From a perturbative point of view, these insertions form the basis of the large-$`N_f`$ expansions, which are widely used at present. In principle, the coefficients of these expansions are scheme-dependent (for the discussion of the scheme dependence of renormalon contributions, see Ref. ), but in what follows this problem will not be analysed since the $`\overline{\mathrm{MS}}`$ scheme will be essentially used.

The next step in the renormalon calculus is the construction, from a perturbative expansion of the physical quantity $`D(a)`$

$$D(a)=\underset{n0}{}d_na^{n+1}$$
(2)
of the Borel transform

$$B[D](\delta )=\underset{n0}{}d_n\frac{\delta ^n}{n!},$$
(3)
and of the Borel integral

$$\stackrel{~}{D}(a)=_0^{\mathrm{}}e^{\delta /a}B[D](\delta )𝑑\delta .$$
(4)
The basic step in defining this integral is related to the analysis of the singularities of the Borel transform $`B[D(\delta )]`$. In the case when it has no singularities on the real positive axis $`\delta `$, the integral of Eq. (3) is well defined. The singularities, located on the negative axis $`\delta `$, are called ultraviolet (UV) renormalons; after Borel resummation, they generate the sign-alternating perturbative QCD series with factorially increasing coefficients. The singularities, located on the positive axis $`\delta `$ are called infrared  (IR) renormalons. These are associated with the sign constant QCD perturbative expansions with coefficients, increasing like $`n!`$. In this case the Borel integral is ill-defined, but its Cauchy Principle Value (PV) can be used (for one of the first applications to the resummation of the renormalon-chain, perturbative contributions to the $`(g2)/2`$ of electron, in QED  with coefficients growing like $`(+1)^nn!`$, see Ref. ). In fact, it is possible to show that the IR renormalon singularities can be associated to power-suppressed corrections of order $`(\mathrm{\Lambda }/Q)^{(2p)}`$, where $`p`$ is the number, related to the pole singularity, at the points $`\delta `$ of the positive axis of the Borel plane (see e.g. Ref. ).

Up to now all discussions were made within the framework of large-$`N_f`$ expansion, which is well-defined in the cases of both QED and QCD. It should be stressed that, since the first coefficient of the QED $`\beta `$-function is proportional to $`N_f`$, the large-$`N_f`$ limit of QED is identical to the large-$`\beta _0`$ approximation. In the case of QCD the situation is more ambiguous. Indeed, in order the rewrite the large-$`N_f`$ perturbative series in the form of the large-$`\beta _0`$ approximation the naive nonabelianization ansatz , namely the substitution $`N_fN_f11N_c/2`$=$`6\beta _0`$ should be used. This procedure turned out to be rather useful in the analysis of both renormalon effects and estimates of the uncertainties due to perturbative corrections beyond the calculated ones (see e.g. Ref. ).

## 3 Definitions of the basic quantities

Consider now the correlator of scalar quark currents $`J_S=m(\mu )\overline{\psi }\psi `$, namely the Green function

$$\mathrm{\Pi }_S(Q^2)=ie^{iqx}J_S(x)J_S(0)_0d^4x.$$
(5)
Its imaginary part $`\mathrm{Im}\mathrm{\Pi }_S=2\pi sR_S\left(1+O(1/s)\right)`$ is related to the spectral function $`R_S(s)`$ of two Euclidean constructs, considered in the work of Ref. :

$$\stackrel{~}{D}_S(Q^2)=Q^2_0^{\mathrm{}}\frac{R_S(s)ds}{(s+Q^2)^2}$$
(6)
$$\overline{D}_S(Q^2)=2Q^2_0^{\mathrm{}}\frac{sR_S(s)ds}{(s+Q^2)^3}.$$
(7)
The first one is defined through the dispersion relation for the first derivative of $`\mathrm{\Pi }_S(Q^2)/Q^2`$, while the second one is defined by the second derivative for $`\mathrm{\Pi }_S(Q^2)`$. The $`\stackrel{~}{D}(Q^2)`$-function was used in the process of calculating the 2-loop and 3-loop QCD corrections to $`R_S(s)`$ in the $`\overline{\mathrm{MS}}`$ scheme in Refs. -. The 4-loop QCD logarithmic corrections to $`\mathrm{\Pi }_S`$, and therefore the $`\alpha _s^3`$ contributions to its imaginary part $`R_S(s)`$, were evaluated in Ref. . In the process of these calculations the following inhomogeneous renormalization group equation was used:

$$\left(\frac{}{\mathrm{log}\mu ^2}+\beta (\alpha _s)\frac{}{\alpha _s}+2\gamma _m(\alpha _s)\right)\mathrm{\Pi }_S=[m(\mu ^2)]^2Q^2\left[\gamma _{SS}(\alpha _s)+O\left(\frac{m^2}{Q^2}\right)\right].$$
(8)
In leading order in $`1/Q^2`$, and next-to-leading order in $`1/N_f`$, we have :

$`\mathrm{\Pi }_S`$ $`=`$ $`[m(\mu ^2)]^2Q^2d_F\left(2L4+{\displaystyle \frac{C_Fb}{T_fN_f}}H(L,b)+O\left({\displaystyle \frac{1}{N_f^2}}\right)\right)`$ (9)
$`\gamma _{SS}`$ $`=`$ $`d_F\left(2+{\displaystyle \frac{C_Fb}{T_fN_f}}h(b)+O\left({\displaystyle \frac{1}{N_f^2}}\right)\right),`$
where $`d_F=3`$, $`C_F=4/3`$, $`L=\mathrm{ln}(\mu ^2/\mathrm{Q}^2)`$ and $`b=T_fN_f\alpha _s/(3\pi )`$. The expression of the coefficients $`H_n(L)`$ of the large-$`N_f`$ expansion of $`H(L,b)`$, namely $`H(L,b)`$=$`_{n>1}H_n(L)b^{n2}`$, can be obtained from the all-order solution of the following relation :

$`n(n1)H_n(L)`$ $`=`$ $`n\left(h_{n+1}4(L+2)g_n\right)+4g_{n+1}9(1)^nD_n(L)`$ (10)
$`{\displaystyle \underset{n}{}}{\displaystyle \frac{D_n(L)\delta ^n}{n!}}`$ $`=`$ $`\left[1+\delta G_D(\delta )\right]\mathrm{exp}\left(\left(\mathrm{L}+{\displaystyle \frac{5}{3}}\right)\delta \right),`$
where $`h_n`$ are defined by the large-$`N_f`$ expansion of the $`h(b)`$-function in the expression for $`\gamma _{SS}`$:

$$h(b)=\underset{n1}{}\left(T_fN_f\frac{\alpha _s}{3\pi }\right)^{n2}\left[h_n+O\left(\frac{1}{N_f}\right)\right]$$
(11)
and $`g_n`$ are the numbers that are generated through the $`ϵ`$-expansion of the following identity

$$\underset{n}{}g_nϵ^{n1}=\left[4\underset{n>1}{}\left(\frac{3}{2^n}+\frac{n}{2}\right)ϵ^{n2}\right]\mathrm{exp}\left(\underset{l>2}{}\frac{2^l3(1)^l}{l}\zeta (l)ϵ^l\right).$$
(12)
## 4 Unravelling the $`\delta =1`$ IR renormalon knot

The contributions of both UV and IR renormalons to $`\mathrm{\Pi }_S`$ is given by the corresponding expression for $`G(\delta )`$ :

$$G_D(\delta )=\frac{2}{1\delta }\frac{1}{2\delta }+\frac{8(1\delta )}{3}\underset{k=2}{\overset{\mathrm{}}{}}\frac{(1)^kk}{(k^2(1\delta )^2)^2}$$
(13)
It can be separated into two pieces as $`G_D(\delta )=G_{}(\delta )+G_+(\delta )`$. The Borel resummable UV renormalons are contained in the first term

$$G_{}(\delta )=\frac{2}{3}\underset{k>0}{}\frac{(1)^k}{(k+\delta )},$$
(14)
with singularities located at $`\delta <0`$. At the large-$`N_f`$, limit it generates the sign-constant asymptotic series. After performing a naive nonabelianization, namely after the replacement $`b\beta _0\alpha _s(\mu ^2)/\pi `$, one obtains sign-alternating asymptotic perturbative QCD series for $`\mathrm{\Pi }_S`$.

The IR renormalons are defined by $`G_+(\delta )`$, which has the following expression

$$G_+(\delta )=\frac{2}{1\delta }\frac{1}{2\delta }+\frac{2}{3}\underset{k>2}{}\frac{(1)^k}{(k\delta )^2}.$$
(15)
Contary to the vector case (see e.g. Ref. ) one can observe the appearance in Eq. (15) of the first IR renormalon at $`\delta =1`$. In accordance with the general rules of renormalon calculus, this IR renormalon generates the $`\mathrm{\Lambda }^2/Q^2`$ contribution to the $`\stackrel{~}{D}_S(Q^2)`$-function of Eq. (6). It should be stressed that this term cannot appear within the operator-product expansion machinery. Indeed, the first higher-twist operator, contributing to the $`\stackrel{~}{D}_S(Q^2)`$-function in the massless scheme comes from the gluon condensate $`(\alpha _s/\pi )G^2_0`$ , related to the IR renormalon pole at $`\delta =2`$ (note that the perturbative QCD correction to the gluon condensate was calculated in Ref. ).

The reason for the appearance of the $`\delta =1`$ IR renormalon contribution to the $`\stackrel{~}{D}_S`$-function is related to the fact that the dispersion relation of Eq. (6) is ill-defined. Indeed, contrary to the vector case, where the Ward identity allows us to fix the expression for $`\mathrm{\Pi }_V(Q^2)`$ at zero transferred momentum, the expression for $`\mathrm{\Pi }_S(0)`$ is infinite. Therefore, the dispersion relation for the first derivative of $`\mathrm{\Pi }_S(Q^2)/Q^2`$ of Eq. (6) contains an ambiguity of order $`\mathrm{\Lambda }^2/Q^2`$. For large momentum transfer this problem can be neglected. However, one should avoid the application of Eq. (6) in the process of long-distance studies, e.g. the determination of light-quark masses values. In this case it is necessary to use the dispersion relation for the second derivative of $`\mathrm{\Pi }_S`$ (see Eq. (7)), which is free from unphysical $`\mathrm{\Lambda }^2/Q^2`$ term. This renormalon rediscovery of the advantages of the twice-differentiated Euclidean construct of Eq. (7), which was originally used in Ref.  and later on in Refs. , is setting its application for the extraction of light-quark masses on a more solid background. It should also be mentioned, that the spectral density $`R_S(s)`$ of both Eq. (6) and Eq. (7) is also free from contributions of $`\delta =1`$ IR renormalon. Moreover, since the factor $`\pi \delta /\mathrm{sin}(\pi \delta )`$ of analytic continuation of the Euclidean constructs to the Minkowskian region is removing all single poles in Eq. (15) and transforming double poles in $`\delta `$ into single poles, the Borel image of $`R_S(s)`$ contains $`\delta >2`$ IR renormalons only, which are associated with the contributions to $`R_S(s)`$ of the operators $`O_k`$ with dimension $`d_k6`$.

## 5 Correlator of scalar quark currents and <br>large-$`\beta _0`$ expansion

The large-$`\beta _0`$ expansion in the $`\overline{\mathrm{MS}}`$ scheme is widely used to generate perturbative series from the direct calculations of the renormalon-chain diagrams, which contain the information on both UV and IR renormalons structure of the related Borel transforms. Let us summarize the results of its application for modelling the behaviour of the perturbative expansions for basic quantities, related to the correlators of the scalar and pseudoscalar quark currents (see Ref. ). The quantities defined above in Eqs. (6),(7) have the following expansions:

$`\stackrel{~}{D}_S(Q^2)`$ $`=`$ $`3[m(Q^2)]^2\left(1+{\displaystyle \underset{n>0}{}}d_na^n\right)`$ (16)
$`\overline{D}_S(Q^2)`$ $`=`$ $`3[m(Q^2)]^2\left(1+{\displaystyle \underset{n>0}{}}\overline{d}_na^n\right)`$ (17)
$`R_S(s)`$ $`=`$ $`3[m(s)]^2\left(1+{\displaystyle \underset{n>0}{}}s_na^n\right).`$ (18)
In leading order of the large-$`\beta _0`$ expansion, the coefficients of these series, obtained using the naive nonabelianization procedure, can be presented as

$`d_n^{NNA}`$ $`=`$ $`\beta _0^{n1}\stackrel{~}{\mathrm{\Delta }}_{n+1}`$ (19)
$`s_n^{NNA}`$ $`=`$ $`\beta _0^{n1}\mathrm{\Delta }_{n+1}`$ (20)
$`\overline{d}_n^{NNA}`$ $`=`$ $`\beta _0^{n1}\overline{\mathrm{\Delta }}_{n+1}`$ (21)
where the coefficients $`\stackrel{~}{\mathrm{\Delta }}_{n+1}`$, $`\mathrm{\Delta }_{n+1}`$ and $`\overline{\mathrm{\Delta }}_{n+1}`$ were calculated in Ref.  in the $`\overline{\mathrm{MS}}`$ scheme. The concrete results of calculations are presented in Table 1.

Table 1: The values of coefficients in Eqs. (19),(20),(21).

$$\begin{array}{cccc}& & & \\ n\hfill & \hfill \stackrel{~}{\mathrm{\Delta }}_n& \hfill \mathrm{\Delta }_n& \hfill \overline{\mathrm{\Delta }}_n\\ & & & \\ 2\hfill & \hfill 7.33333& \hfill 7.33333& \hfill 5.33333\\ 3\hfill & \hfill 10.4696& \hfill 7.17968& \hfill 3.13622\\ 4\hfill & \hfill 32.6145& \hfill 8.48885& \hfill 11.6754\\ 5\hfill & \hfill 97.9534& \hfill 4.36402& \hfill 0.10978\\ 6\hfill & \hfill 503.887& \hfill 2.96849& \hfill 112.074\\ 7\hfill & \hfill 2194.28& \hfill 54.2101& \hfill 325.157\\ 8\hfill & \hfill 16465.8& \hfill 123.639& \hfill 3300.11\end{array}$$

One can see that the coefficients $`\stackrel{~}{\mathrm{\Delta }}_{n+1}`$ are increasing very fast. This feature is related to the manifestation of $`\delta =1`$ IR renormalon contribution to Eq. (13). The suppression of this singularity in the expressions for $`R_S(s)`$ and $`\overline{D}_S(Q^2)`$ leads to a rather well-behaved perturbative series up to $`n=6`$-loop order for $`R_S(s)`$ and $`n=5`$-loop order in the latter case, where the corresponding coefficient $`\mathrm{\Delta }_5`$ is remarkably small. The factorial growth of the related series starts to manifest itself at the level of higher-order corrections.

Consider now the uncertainties of the large-$`\beta _0`$ series, comparing the numerical values for the coefficients of Eqs. (19),(20),(21) with the explicit numbers, obtained as the result of analytical calculations, performed in Refs.  (see Table 2).

Table 2: Ratios of NNA estimates to exact results.

$$\begin{array}{cccc}& & & \\ Coeff.\hfill & \hfill N_f=3& \hfill N_f=4& \hfill N_f=5\\ & & & \\ d_2^{\mathrm{NNA}}/d_2\hfill & \hfill 0.514& \hfill 0.496& \hfill 0.477\\ s_2^{\mathrm{NNA}}/s_2\hfill & \hfill 0.507& \hfill 0.490& \hfill 0.472\\ \overline{d}_2^{\mathrm{NNA}}/\overline{d}_2\hfill & \hfill 0.498& \hfill 0.484& \hfill 0.469\\ & & & \\ d_3^{\mathrm{NNA}}/d_3\hfill & \hfill 0.354& \hfill 0.346& \hfill 0.339\\ s_3^{\mathrm{NNA}}/s_3\hfill & \hfill 0.482& \hfill 0.565& \hfill 0.747\\ \overline{d}_3^{\mathrm{NNA}}/\overline{d}_3\hfill & \hfill 0.764& \hfill 0.871& \hfill 1.081\end{array}$$

One can see that contrary to the applications of the scheme-invariant approach for estimates of higher-order perturbative QCD corrections in the scalar channel , the large-$`\beta _0`$ approximation gives the correct signs for the coefficients of all considered quantities. The reason for the failure of the scheme-invariant estimates at the 4-loop level will be clarified below. It should also be stressed that although the large-$`\beta _0`$ predictions are of the same order of magnitude as the exactly calculable coefficients, the application of the naive nonabelianization procedure in the scalar channel has the tendency to underestimate the concrete perturbative terms by a factor of order 0.5. Keeping this in mind, and taking into account the results from Table 1, we arrive at the conclusion that the large-$`\beta _0`$ approximation indicates, that the 5-loop corrections to the perturbative series for the $`\overline{D}_S`$ function and $`R_S(s)`$, which enter in the determination of the strange-quark mass within QCD sum rules method (see Refs. ,\- ), are extremely small. Indeed, fixing $`N_f=3`$ and $`\alpha _s0.3`$ we obtain

$`\overline{D}_S`$ $`=`$ $`1+3.67\left({\displaystyle \frac{\alpha _s}{\pi }}\right)+14.17\left({\displaystyle \frac{\alpha _s}{\pi }}\right)^2+77.36\left({\displaystyle \frac{\alpha _s}{\pi }}\right)^3+2\times 1.26\left({\displaystyle \frac{\alpha _s}{\pi }}\right)^4`$
$`=`$ $`1+0.350+0.129+0.067+0.0002`$
and

$`R_S(s)`$ $`=`$ $`1+5.667\left({\displaystyle \frac{\alpha _s}{\pi }}\right)+31.864\left({\displaystyle \frac{\alpha _s}{\pi }}\right)^2+89.156\left({\displaystyle \frac{\alpha _s}{\pi }}\right)^3+98\left({\displaystyle \frac{\alpha _s}{\pi }}\right)^4`$
$`=`$ $`1+0.541+0.291+0.078+0.008.`$
Therefore, large-$`\beta _0`$ estimates indicate that the doubts about a serious underestimate of the errors of the QCD strange-quark-mass extraction from the pseudoscalar correlator at the $`(\alpha _s^2)`$–order , formulated in Ref.  using the positivity of a spectral function for a pseudoscalar correlator, cannot be related to the non-applicability of the perturbative QCD expansions at rather small energy scale. Moreover, the $`O(\alpha _s^3)`$ result of Ref.  $`\overline{m}_s(1\mathrm{GeV})=205\pm 19\mathrm{M}\mathrm{e}\mathrm{V}`$ is closer to the upper bound from Ref. , namely $`140\overline{m}_s(1\mathrm{GeV})254`$ MeV. At the energy scale of 2 GeV the results of Ref.  are equivalent to the value $`\overline{m}_s(2\mathrm{GeV})`$ 148 MeV, which are higher than the current lattice results (for a review, see Ref. ). Note, however, that the current determination of light-quark masses from QCD sum rules for the scalar and pseudoscalar two-point functions does not touch the problem of the possible manifestation of the instanton contributions in the theoretical part of these sum rules (see Ref. ). This is why the determination of the s-quark mass in Cabibbo-suppressed $`\tau `$-lepton decay characteristics, which seem to be free from instanton contributions, attracts more and more attention (see e.g. Ref.  and in particular Ref. ). Note that the central value of the s-quark mass result of Ref. , namely $`\overline{m}_s(1\mathrm{GeV})=176\pm 37_{\mathrm{exp}}\pm 13_{\mathrm{th}}`$ MeV, which is equivalent to $`\overline{m}_s(2\mathrm{GeV})`$ 128 MeV, is lower than of the one of Ref.  and is closer to the lattice results, summarized in the review of Ref. . One of the reasons can be related to the approach applied in Ref. , of resummation of the infinite subset of $`\pi ^2`$-dependent effects that arize from analytical continuation of the mass-dependent quantities from the Euclidean to the Minkowskian region. Within the large-$`\beta _0`$ approximation the special features of this procedure were studied in Ref. .

## 6 Pole vs running quark masses and <br>the large-$`\beta _0`$ approximation

Another example where the naive nonabelianization procedure  and the large-$`\beta _0`$ estimates are working reasonably well, is the relation between pole and running quark masses. Indeed, this relation was recently evaluated analytically at the $`\alpha _s^3`$-order . The results of these calculations, which confirm those obtained in Ref.  in a semi-analytical way, have the following numerical form :

$`M`$ $`=`$ $`\overline{m}(\overline{m})[1+{\displaystyle \frac{4}{3}}\left({\displaystyle \frac{\overline{\alpha }_s}{\pi }}\right)+\left({\displaystyle \frac{\overline{\alpha }_s}{\pi }}\right)^2(1.0414N_f+13.4434)`$
$`+`$ $`\left({\displaystyle \frac{\overline{\alpha }_s}{\pi }}\right)^3(0.6527N_f^226.655N_f+190.595)]`$
with $`\overline{\alpha }_s=\overline{\alpha }_s(\overline{m})`$. This series can be rewritten in terms of the scheme-independent coefficients of the QCD $`\beta `$-function as :

$`M`$ $`=`$ $`\overline{m}(\overline{m})[1+{\displaystyle \frac{4}{3}}\left({\displaystyle \frac{\overline{\alpha }_s}{\pi }}\right)+\left({\displaystyle \frac{\overline{\alpha }_s}{\pi }}\right)^2(6.248\beta _03.739)`$
$`+`$ $`\left({\displaystyle \frac{\overline{\alpha }_s}{\pi }}\right)^3(23.497\beta _0^2+6.248\beta _1+1.019\beta _029.94)],`$
where $`\beta _0`$ is as defined above and $`\beta _1=[102(38/3)N_f]/16`$. It is possible now to deduce that the large-$`\beta _0`$ approximation works reasonably well in this case as well, thanks to the cancellations between the term proportional to the $`\beta _1`$ coefficient and the conformal $`\beta `$-independent one, in the order correction of order $`\alpha _s^3`$ to Eq. (6) the large–$`\beta _0`$ approximation is working reasonably well <sup>1</sup><sup>1</sup>1I wish to thank E. Gardi for discussions on questions related to this topic.. This feature is also illustrated in Table 3, taken from Ref. , where the comparison of the results of semi-analytic calculations of the $`O(\alpha _s^3)`$ correction to Eq. (24) , presented in the second column, are compared with the estimates of these terms, made in Ref.  using the scheme-invariant methods (namely the effective-charges approach (ECH) of Ref.  and the PMS approach of Ref. ) and the large-$`\beta _0`$ estimates of Ref. .

Table 3: Comparison of the results of Ref.  paper with estimates based on ECH, PMS and the large-$`\beta _0`$ approximation for $`M/\overline{m}(\overline{m})`$.

| $`N_f`$ |  |  (ECH) |  (PMS) |  (large-$`\beta _0`$) |
| --- | --- | --- | --- | --- |
| $`2`$ | $`143(3)`$ | $`152.71`$ | $`153.76`$ | $`137.23`$ |
| $`3`$ | $`119(3)`$ | $`124.10`$ | $`124.89`$ | $`118.95`$ |
| $`4`$ | $`96(2)`$ | $`97.729`$ | $`98.259`$ | $`101.98`$ |
| $`5`$ | $`75(2)`$ | $`73.616`$ | $`73.903`$ | $`86.318`$ |

One can see that the scheme-invariant methods, developed in other cases in Ref. , also give reasonable estimates. However, there are cases, where the comparison of the results of estimates, produced by scheme-invariant methods or the large $`\beta _0`$-approach, with the explicitly calculated corrections of the perturbative series under consideration, demonstrate larger ambiguities.

## 7 Limitations of estimating procedures: when and why?

It should be mentioned that unlike the case of the qualitative success of the large-$`\beta _0`$ estimates in Ref.  of the coefficients of the perturbative series of quantities related to the Green function of the scalar and pseudoscalar quark currents, the application of the scheme-invariant estimates in this case result in the appearance of large and negative- order $`\alpha _s^4`$ corrections in the expression for $`R_S(s)`$ . It is thus worth-while trying to understand the reason for the appearance of these suspicious numbers and, therefore, the limitations of the applications of the scheme-invariant procedure of Ref. , used in Ref. . These studies were done in Ref.  and their main conclusions will be summarized here.

It is worth-while recalling that the approach of Ref.  is constructed for the estimates of higher-order perturbative corrections to renormalization- group-invariant quantities of Eq. (2), which obey the renormalization-group equation without anomalous dimension term, with the $`\beta `$-function defined as $`\beta (a)=_{i0}\beta _ia^{i+2}=\beta _0a^2(1+_{i1}c_ia^i)`$, where $`c_i=\beta _i/\beta _0`$. It is known that the method of the effective charges of Ref.  prescribes to define the scheme where all higher-order corrections to $`D(a)`$ will be nullified and the corresponding expression will have the following form: $`D(\stackrel{~}{a})=d_0\stackrel{~}{a}`$. The new expansion parameter obeys the renormalization-group equation with the effective-charges $`\beta `$-function, defined as

$$\stackrel{~}{\beta }(\stackrel{~}{a})=\underset{n0}{}\stackrel{~}{\beta }_n\stackrel{~}{a}^{i+2},$$
(26)
where the coefficients $`\stackrel{~}{\beta }_n`$ are scheme-invariant. The basic formula for generating the estimates from the scheme-invariant procedure of Ref.  is:

$$\frac{\stackrel{~}{\beta }_n\beta _n}{(n1)}=\beta _0(d_n\mathrm{\Omega }_n),$$
(27)
where $`\mathrm{\Omega }_n`$ are defined as

$`\mathrm{\Omega }_2`$ $`=`$ $`d_1(c_1+d_1)`$ (28)
$`\mathrm{\Omega }_3`$ $`=`$ $`d_1\left(c_2{\displaystyle \frac{1}{2}}c_1d_12d_1^2+3d_2\right)`$
$`\mathrm{\Omega }_4`$ $`=`$ $`d_1\left(c_3{\displaystyle \frac{4}{3}}c_2d_1+{\displaystyle \frac{2}{3}}c_1d_2+{\displaystyle \frac{14}{3}}d_1^3{\displaystyle \frac{28}{3}}d_1d_2+4d_3\right)`$ (29)
$`+{\displaystyle \frac{1}{3}}\left(c_2d_2c_1d_3+5d_2^2\right).`$
at 5 loops. The idea that lies beyond the procedure of Ref.  is that reasonable estimates of the higher-order coefficients $`d_n`$ of the perturbative series for physical quantities can be obtained if the coefficients $`\stackrel{~}{\beta }_n`$ of the process- dependent, effective-charges $`\beta `$-functions are of the same sign and magnitude as the coefficients $`\beta _n`$ of the $`\beta `$-finction in the $`\overline{\mathrm{MS}}`$ scheme. This assumption turns out to be true at the 4-loop level in QCD for the perturbative series of deep-inelastic scattering sum rules and the $`e^+e^{}`$-annihilation $`D`$-function, in the cases $`N_f=4`$ and $`N_f=3`$ in particular. Indeed the numerical analogues of the analytical expressions for the 4-loop QCD $`\beta `$-function  have the following form

$`\beta (N_f=3)`$ $`=`$ $`2.250a^24.000a^310.060a^447.228a^5`$ (30)
$`\beta (N_f=4)`$ $`=`$ $`2.083a^23.208a_s^36.349a^431.387a^5`$ (31)
$`\beta (N_f=5)`$ $`=`$ $`1.917a^22.417a^32.827a^418.852a^5`$ (32)
with $`a\alpha _s/\pi `$. It should be stressed that at $`N_f=3`$ and $`N_f=4`$ the effective-charges $`\beta `$-functions of the polarized and unpolarized Bjorken sum rules and of the $`e^+e^{}`$-annihilation $`D`$-function have a similar 4-loop behaviour, provided the explicitly unknown coefficients $`d_3`$ in $`\mathrm{\Omega }_4`$ are modelled by the numerical expressions for $`\mathrm{\Omega }_3`$. However, this is not true in the case of renormalization-group-invariant quantities, related to the correlators of scalar and pseudoscalar currents. Indeed, for the renormalization-group-invariant quantity related to the $`\overline{D}_s`$-function of Eq. (7) as

$$\overline{R}_D(Q^2)=\frac{1}{2}\frac{d\mathrm{log}\overline{\mathrm{D}}_\mathrm{S}(\mathrm{Q}^2)}{d\mathrm{logQ}^2}$$
(33)
the behaviour of the effective-charges $`\beta `$-function in the $`O(\stackrel{~}{a}^3)`$ and $`O(\stackrel{~}{a}^3)`$-approximations differs significantly from the behaviour of Eqs. (30),(31),(32) . The coresponding expressions

$`\stackrel{~}{\beta }(N_f=3)`$ $`=`$ $`2.250\stackrel{~}{a}^24.000\stackrel{~}{a}^3+58.920\stackrel{~}{a}^42148.503\stackrel{~}{a}^5`$ (34)
$`\stackrel{~}{\beta }(N_f=4)`$ $`=`$ $`2.083\stackrel{~}{a}^23.208\stackrel{~}{a}^3+53.852\stackrel{~}{a}^41687.191\stackrel{~}{a}^5`$ (35)
$`\stackrel{~}{\beta }(N_f=5)`$ $`=`$ $`1.917\stackrel{~}{a}^22.417\stackrel{~}{a}^3+49.356\stackrel{~}{a}^41303.490\stackrel{~}{a}^5`$ (36)
indicate the appearance of a spurious perturbative infrared fixed point. This spurious zero is compensated by the appearance of large and negative $`\stackrel{~}{a}^5`$ contributions. In the case of the $`R_D`$-function, similar features were already observed at the 3-loop level in Ref.  and at the 4-loop in Ref.  both in the Minkowskian region; they were studied in detail in the Euclidean region in Ref. . Here we conclude that these features demonstrate the limitations of the application of the scheme-invariant methods to the estimates of the O($`\alpha _s^4`$)-corrections to the $`\stackrel{~}{D}_S`$ function of Eq. (6), the $`\overline{D}_S`$ function of Eq. (7), and their spectral density $`R_S(s)`$.

However, the approximation of the large $`\beta _0`$-expansion, which we think is working reasonably well in the case of the correlator of the scalar and pseudoscalar quark currents , also has definite limitations. The inapplicability of this method for the estimates of higher-order corrections to the DGLAP non-singlet kernel, studied in detail in Ref. , can be explained by the absence of renormalon-type contributions to the anomalous dimensions of the non-singlet operators. A more interesting fact is related to the contradiction between the large-$`\beta _0`$ estimates of the high-order perturbative corrections to the coefficient functions of $`N=2`$ and $`N=4`$ non-singlet Mellin moments, for the $`F_2`$ structure function , and the result of explicit calculations performed in Ref. . The limitations of the predictive possibilities of the large-$`\beta _0`$ method in this case are still unclear.

Acknowledgements

I am grateful to D.J. Broadhurst and C.J. Maxwell for the productive and pleasant collaboration, which resulted in the publication of Ref. . Different aspects of this work were discussed at the Quarks-2000 International Seminar in the talks of D.J.B. and A.L.K., and are described in part above. It is a pleasure to thank G. Altarelli, V.M. Braun, J. Ellis, J. Fischer, N.V. Krasnikov, S.A. Larin, N. Paver, A.A. Pivovarov, A. Schafer and F.J. Yndurain for discussions and constructive questions and comments.

Part of the work on the problems approached in this paper was done in Russia within the framework of scientific projects of RFBR Grants N 99-01-00091, N 00-02-1742.


text:
# Interferometry with Entangled Atoms

## Abstract

A quantum gravity-gradiometer consists of two spatially separated ensembles of atoms interrogated by pulses of a common laser beam. Laser pulses cause the probability amplitudes of atomic ground-state hyperfine levels to interfere, producing two motion-sensitive phase shifts which allow the measurement of the average acceleration of each ensemble, and, via simple differencing, of the acceleration gradient. Here I propose entangling the quantum states of the two ensembles prior to the pulse sequence, and show that entanglement encodes their relative acceleration in a single interference phase which can be measured directly, with no need for differencing.

PACS numbers: 39.20.+q, 03.75.Dg, 03.67.-a, 03.67.Lx

1. Introduction

Inertial sensors based on matter-wave interferometry are a rapidly developing technology \[1—4\]. One class of such sensors are gradiometers designed to measure linear acceleration gradients, typically gradients caused by inhomogeneities in the gravitational potential. The state-of-the-art design for an atom-wave gravity gradiometer consists of a pair of atom-wave accelerometers separated spatially by a fixed distance and direction. Each of the two accelerometers contains an ensemble of laser-cooled atoms, through which a common pair of (counter-propagating) Raman laser beams are pulsed in a carefully controlled sequence to drive Rabi oscillations between atomic ground-state hyperfine levels. The sequence and timing of the laser pulses are adjusted so that the hyperfine ground and excited states of any single atom interfere during the atom’s motion through the laser field. At the end of the sequence, the probability of finding an atom in its excited hyperfine state is a simple function of a relative phase, which accumulates between successive pulses in proportion to the atom’s acceleration (relative to the inertial frame defined by the laser beam) averaged over its flight time. Subsequently, the resulting pair of observations of this phase (one for each ensemble) can be subtracted to obtain a measurement of the relative acceleration, and, upon dividing the result by the separation length, a measurement of the acceleration gradient. I will give a more quantitative review of this quantum interferometry technique in Sect. 2 below. A detailed description of the experimental setup can be found in  (see esp. Fig. 1 there).

Even though subtracting the two acceleration measurements allows a number of correlated noise sources to cancel out as common modes, it is still desirable to avoid such differencing of two nearly-equal measured quantities which are corrupted by (uncorrelated) noise. In this letter, I will show that by prior-entangling (prior to the laser pulsing sequence) the quantum states of atoms in the two ensembles (such that every entangled pair has one atom in each ensemble) it becomes possible to encode the relative-acceleration information directly in a single interference phase, thereby eliminating the need for the differencing of two separate phase measurements. This idea is inspired by the recent discovery of a similar quantum algorithm which uses entangled states to synchronize atomic clocks non-locally. The Quantum Clock Synchronization (QCS) algorithm relies on preparing the atoms “stabilizing” each atomic clock in a special entangled quantum state whose time evolution reduces to a pure multiplicative phase as long as each atom in the pair evolves under the same unitary transformation. By contrast, the Entangled Quantum Interferometry (EQI) algorithm to be described below relies on transforming the pair’s quantum state into a sensor which senses the difference between unitary evolutions of the entangled atoms. Put another way, the EQI algorithm is a spatial analogue of the QCS algorithm, in exactly the same sense as atomic clocks are temporal analogues of quantum (atom-wave) interferometers.

2. Overview of matter-wave gravity gradiometry

Consider an atom in its internal ground state $`|0`$, moving through the inertial frame defined by the Raman beams of an atom interferometer. \[Here “inertial frame” has the same meaning as in general relativity: in the presence of a gravitational field, my coordinate system will always be locally-Minkowski (free-falling, or geodesic-normal coordinates) so that the gravitational potential and its first derivatives vanish at the origin.\] The wave function can be written as a tensor product of the atom’s internal and external Hilbert-space states:

$`\mathrm{\Psi }_{0:\stackrel{}{p}}(\stackrel{}{x},t)`$ $`=`$ $`|0{\displaystyle a(\stackrel{}{k},\omega )e^{i(\stackrel{}{k}\stackrel{}{x}\omega t)}d^3k}`$ (1)
$``$ $`|0\psi (\stackrel{}{x},t),`$ (2)
where $`a(\stackrel{}{k},\omega )`$ defines the wave packet representing the external (space-time) part $`\psi (\stackrel{}{x},t)`$ of the quantum state, $`\omega `$ is related to $`\stackrel{}{k}`$ via the usual dispersion relation (plus an additive constant corresponding to the atom’s internal energy level), and $`\stackrel{}{p}`$ denotes the average (in general time-dependent) momentum:

$$\stackrel{}{p}\psi |i\mathrm{}|\psi =(2\pi )^3\mathrm{}\stackrel{}{k}|a(\stackrel{}{k},\omega )|^2d^3k.$$
(3)
The need for this cumbersome subscript notation $`\mathrm{\Psi }_{0:\stackrel{}{p}}`$ will become clear in a moment. Since the two Raman beams are detuned in frequency by an amount $`\mathrm{\Omega }`$ which corresponds to the energy difference between the ground and excited states $`|0`$ and $`|1`$, they might stimulate the ground-state atom to absorb a photon from one beam and emit a lower frequency photon into the other, resulting in a transition to the excited state $`|1`$. This stimulated transition gives the atom a net momentum kick in the amount $`\mathrm{}(\stackrel{}{k}_1\stackrel{}{k}_2)\mathrm{}\stackrel{}{K}`$, where $`\stackrel{}{k}_1`$, $`\stackrel{}{k}_2`$ ($`|\stackrel{}{k}_1|>|\stackrel{}{k}_2|`$) denote the wave vectors corresponding to the (counter-propagating) Raman beams. \[As the vectors $`\stackrel{}{k}_i`$, $`i=1,\mathrm{\hspace{0.33em}2}`$, point in opposite directions, the magnitude of $`\stackrel{}{K}`$ is typically twice that of the $`\stackrel{}{k}_i`$. Note that $`\mathrm{\Omega }`$ is the total energy absorbed by the atom, which includes not only the transition energy (within the appropriate bandwidth), but also the recoil contribution.\] The wave function of the same atom excited in this way from its ground state can then be decomposed, similarly to Eq. (1), in the form

$`\mathrm{\Psi }_{1:\stackrel{}{p}+\mathrm{}\stackrel{}{K}}(\stackrel{}{x},t)`$ $`=`$ $`|1{\displaystyle a(\stackrel{}{k},\omega )e^{i[(\stackrel{}{k}+\stackrel{}{K})\stackrel{}{x}(\omega +\mathrm{\Omega })t]}d^3k}`$ (4)
$`=`$ $`|1\psi (\stackrel{}{x},t)e^{i(\stackrel{}{K}\stackrel{}{x}\mathrm{\Omega }t)}.`$ (5)
Since the wave packet $`\psi (\stackrel{}{x},t)`$ is common to both the ground and excited state wave functions, I will suppress it (as an overall “normalization” factor) in what follows. Note that whenever an atom interrogated in the interferometer is excited from its ground state $`|0:\stackrel{}{p}`$ to the internal excited state $`|1`$ by Raman pulses, it always picks up an associated excess energy-momentum $`(\stackrel{}{K},\mathrm{\Omega })`$ as specified in Eq. (3); in other words, the excited atom’s total quantum state becomes $`|1:\stackrel{}{p}+\mathrm{}\stackrel{}{K}`$. Conversely, a stimulated transition in the reverse direction results in the transformation $`|1:\stackrel{}{p}|0:\stackrel{}{p}\mathrm{}\stackrel{}{K}`$. This correlation between an atom’s internal state and its momentum during the laser interrogation process is the key feature enabling atom interferometry.

In the interferometer, atoms are manipulated by laser pulses of two basic kinds: A Hadamard pulse $`H_{\frac{\pi }{2}}`$ (a $`\pi /2`$ pulse followed by the spin operator $`\sigma _z`$), whose action is

$`H_{\frac{\pi }{2}}|0:\stackrel{}{p}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}(|0:\stackrel{}{p}+|1:\stackrel{}{p}+\mathrm{}\stackrel{}{K}),`$ (6)
$`H_{\frac{\pi }{2}}|1:\stackrel{}{p}+\mathrm{}\stackrel{}{K}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}(|0:\stackrel{}{p}|1:\stackrel{}{p}+\mathrm{}\stackrel{}{K}),`$ (7)
and a “double” Hadamard pulse $`H_\pi `$ (a $`\pi `$ pulse followed by the spin operator $`\sigma _z`$), whose action is

$`H_\pi |0:\stackrel{}{p}`$ $`=`$ $`|1:\stackrel{}{p}+\mathrm{}\stackrel{}{K},`$ (8)
$`H_\pi |1:\stackrel{}{p}+\mathrm{}\stackrel{}{K}`$ $`=`$ $`|0:\stackrel{}{p}.`$ (9)
Consider a single atom’s flight across the laser field during the interrogation phase: Initially, at times $`t<t_1`$, say, the atom is in the ground state $`|0:\stackrel{}{p}`$. At time $`t=t_1`$ the first $`H_{\frac{\pi }{2}}`$ pulse hits, and the atom’s state becomes (up to an overall phase factor which I will always ignore):

$$\mathrm{\Psi }(t_1)=\frac{1}{\sqrt{2}}\left(|0:\stackrel{}{p}+|1:\stackrel{}{p}+\mathrm{}\stackrel{}{K}\right).$$
(10)
Since atoms are well-localized spatially, their wave functions $`\mathrm{\Psi }(\stackrel{}{x},t)`$ are sharply peaked around their average position $`\mathrm{\Psi }(t)|\stackrel{}{x}|\mathrm{\Psi }(t)`$ at all times $`t`$. When an atom is excited by Raman pulses to a superposition of its internal states, its distribution becomes bimodal, but still sharply peaked around the two positions (modes). Therefore, to a very good approximation, I can follow the evolution of the atom’s wave function during its flight in the interferometer by examining it in the vicinity of the atom’s (average) position(s) as a function of time $`t`$. (See  for a more rigorous treatment.) Accordingly, at any time $`t`$ after $`t_1`$ and before the next Raman pulse hits, the state Eq. (6) evolves, according to Eqs. (1)–(3), as

$`\mathrm{\Psi }(t)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left[|0:\stackrel{}{p}+e^{i\stackrel{}{K}(\stackrel{}{x}_t\stackrel{}{x}_1)i\mathrm{\Omega }(tt_1)}|1:\stackrel{}{p}+\mathrm{}\stackrel{}{K}\right],`$ (12)
$`t:t_1t<t_2,`$
where $`\stackrel{}{x}_t\stackrel{}{x}(t)`$, $`\stackrel{}{x}_1\stackrel{}{x}(t_1)`$ etc. is shorthand notation for the atom’s position at the different times, and an overall phase factor as well as the common wave packet $`\psi (\stackrel{}{x},t)`$ are suppressed as advertised. At a time $`t=t_2>t_1`$, an $`H_\pi `$ pulse hits, and the atom’s new state becomes

$`\mathrm{\Psi }(t_2)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\text{[}|1:\stackrel{}{p}+\mathrm{}\stackrel{}{K}`$ (13)
$`+`$ $`e^{i\stackrel{}{K}(\stackrel{}{x}_2\stackrel{}{x}_1)i\mathrm{\Omega }(t_2t_1)}|0:\stackrel{}{p}].`$ (14)
At any time $`t`$ after $`t_2`$ and before the next pulse, the state Eq. (8) evolves as

$`\mathrm{\Psi }(t)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}[e^{i\stackrel{}{K}(\stackrel{}{x}_t\stackrel{}{x}_2)i\mathrm{\Omega }(tt_2)}|1:\stackrel{}{p}+\mathrm{}\stackrel{}{K}`$ (15)
$`+`$ $`e^{i\stackrel{}{K}(\stackrel{}{x}_2\stackrel{}{x}_1)i\mathrm{\Omega }(t_2t_1)}|0:\stackrel{}{p}],`$ (17)
$`t:t_2t<t_3.`$
At time $`t=t_3`$, the second and last $`H_{\frac{\pi }{2}}`$ pulse hits, and transforms $`\mathrm{\Psi }(t)`$ into the state

$`\mathrm{\Psi }(t_3)`$ $`=`$ $`{\displaystyle \frac{1}{2}}[e^{i\stackrel{}{K}(\stackrel{}{x}_3\stackrel{}{x}_2)i\mathrm{\Omega }(t_3t_2)}(|0:\stackrel{}{p}|1:\stackrel{}{p}+\mathrm{}\stackrel{}{K})`$ (18)
$`+`$ $`e^{i\stackrel{}{K}(\stackrel{}{x}_2\stackrel{}{x}_1)i\mathrm{\Omega }(t_2t_1)}(|0:\stackrel{}{p}+|1:\stackrel{}{p}+\mathrm{}\stackrel{}{K})].`$ (19)
So as a result of the pulse sequence $`H_{\frac{\pi }{2}}H_\pi H_{\frac{\pi }{2}}`$, an atom which enters the interferometer in the ground state $`|0:\stackrel{}{p}`$ ends up, at the end of the last $`H_{\frac{\pi }{2}}`$ pulse (at $`t=t_3`$), in the state Eq. (10) which can be rewritten as

$`\mathrm{\Psi }(t_3)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\text{[}\left(e^{i\mathrm{\Theta }_{21}}+e^{i\mathrm{\Theta }_{32}}\right)|0:\stackrel{}{p}`$ (20)
$`+`$ $`(e^{i\mathrm{\Theta }_{21}}e^{i\mathrm{\Theta }_{32}})|1:\stackrel{}{p}+\mathrm{}\stackrel{}{K})],`$ (21)
where for $`i,j=1,2,3,\mathrm{}`$

$$\mathrm{\Theta }_{ij}\stackrel{}{K}(\stackrel{}{x}_i\stackrel{}{x}_j)\mathrm{\Omega }(t_it_j).$$
(22)
After absorbing a factor of $`e^{i\mathrm{\Theta }_{21}}`$ into the overall phase multiplying $`\mathrm{\Psi }(t_3)`$, Eq. (11) takes the more familiar form

$`\mathrm{\Psi }(t_3)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\text{[}\left(1+e^{i\mathrm{\Theta }}\right)|0:\stackrel{}{p}`$ (23)
$`+`$ $`(1e^{i\mathrm{\Theta }})|1:\stackrel{}{p}+\mathrm{}\stackrel{}{K})],`$ (24)
where $`\mathrm{\Theta }\mathrm{\Theta }_{32}\mathrm{\Theta }_{21}`$. The phase $`\mathrm{\Theta }`$ can now be observed by measuring the relative abundance of ground vs. excited state atoms in the state $`\mathrm{\Psi }(t_3)`$. Specifically, the fraction $`P_1`$ of excited-state atoms in the state Eq. (13) is given by

$$P_1=\frac{1}{4}\left|1e^{i\mathrm{\Theta }}\right|^2=\mathrm{sin}^2\left(\frac{\mathrm{\Theta }}{2}\right),$$
(25)
and, similarly, the fraction $`P_0`$ of ground state atoms is $`P_0=|1+e^{i\mathrm{\Theta }}|^2/4=\mathrm{cos}^2(\mathrm{\Theta }/2)`$. On the other hand, substituting Eqs. (12) in $`\mathrm{\Theta }=\mathrm{\Theta }_{32}\mathrm{\Theta }_{21}`$ gives

$$\mathrm{\Theta }=\stackrel{}{K}[\stackrel{}{x}(t_3)2\stackrel{}{x}(t_2)+\stackrel{}{x}(t_1)]\mathrm{\Omega }(t_32t_2+t_1).$$
(26)
If the Raman pulses making up the sequence $`H_{\frac{\pi }{2}}H_\pi H_{\frac{\pi }{2}}`$ are aligned in time such that $`t_2=t_1+T`$ and $`t_3=t_2+T`$ for some common interrogation time $`T`$, the second term in parenthesis in Eq. (15) vanishes, and, using the standard Taylor expansion of $`\stackrel{}{x}(t)`$, we can rewrite the first term in the form

$$\mathrm{\Theta }=\stackrel{}{K}\ddot{\stackrel{}{x}}(t_2)T^2+O(T^3).$$
(27)
Here the magnitude of the $`O(T^3)`$ remainder is down by a factor of order $`T/T_a`$ relative to the first (acceleration) term, where $`T_a\stackrel{}{a}/\dot{\stackrel{}{a}}`$ denotes the timescale over which the atom’s acceleration $`\stackrel{}{a}(t)\ddot{\stackrel{}{x}}(t)`$ varies. Therefore, as long as the time scale $`T_a`$ is much larger (as is typically the case) than the interrogation time $`T`$, measurement of the phase $`\mathrm{\Theta }`$ following the pulse sequence $`H_{\frac{\pi }{2}}{}_{\overline{T}}{}^{}H_{\pi }^{}{}_{\overline{T}}{}^{}H_{\frac{\pi }{2}}^{}`$ yields a direct measurement of the inertial acceleration component $`\stackrel{}{K}\stackrel{}{a}`$.

3. Interferometry with entangled atoms

For a gradient measurement, ordinarily it is necessary to apply two separate but simultaneous measurements of acceleration to two different ensembles of atoms separated spatially by a fixed distance. By using atoms whose internal states are pairwise entangled in a distributed fashion, however, it is possible to obtain the acceleration gradient information with just one phase measurement as I will show now:

Start, at time $`t=t_1`$, with atoms in the two ensembles $`A`$ and $`B`$ in the pairwise entangled initial state

$`\mathrm{\Psi }(t_1)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}(|0:\stackrel{}{p}_A_A|1:\stackrel{}{p}_B+\mathrm{}\stackrel{}{K}_B`$ (28)
$``$ $`|1:\stackrel{}{p}_A+\mathrm{}\stackrel{}{K}_A|0:\stackrel{}{p}_B_B).`$ (29)
In other words, in some subset (which my measurements will be directed at) of ensemble $`A`$, every atom is entangled with some atom in ensemble $`B`$, and, at some initial time $`t=t_1`$, each entangled pair is (upto an overall phase) in the state Eq. (17). Focus attention now on a single such entangled pair. Since the state Eq. (17) is an energy eigenstate of the total internal Hamiltonian, its explicit time evolution consists of a pure multiplicative phase only. Implicitly, it evolves in time kinematically with the inertial motion of the atoms; according to Eqs. (1)–(3), this evolution takes the form

$`\mathrm{\Psi }(t)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}[e^{i\stackrel{}{K}(\stackrel{}{x}_{Bt}\stackrel{}{x}_{B1})}|0:\stackrel{}{p}_A_A|1:\stackrel{}{p}_B+\mathrm{}\stackrel{}{K}_B`$ (30)
$``$ $`e^{i\stackrel{}{K}(\stackrel{}{x}_{At}\stackrel{}{x}_{A1})}|1:\stackrel{}{p}_A+\mathrm{}\stackrel{}{K}_A|0:\stackrel{}{p}_B_B],`$ (32)
$`t:t_1t<t_2,`$
where $`\stackrel{}{x}_{At}\stackrel{}{x}_A(t)`$, $`\stackrel{}{x}_{B1}\stackrel{}{x}_B(t_1)`$, etc., and an overall phase factor as well as the common wave packet $`\psi _A(\stackrel{}{x}_A,t)\psi _B(\stackrel{}{x}_B,t)`$ are suppressed as in Eq. (7) above. For increased readability, let me introduce the following abbreviation which I will use throughout the rest of this letter: for $`E=A,B`$,

$$|0_E|0:\stackrel{}{p}_E_E,|1_E|1:\stackrel{}{p}_E+\mathrm{}\stackrel{}{K}_E.$$
(33)
Now suppose that, at a time $`t=t_2>t_1`$, a common $`H_\pi `$ pulse is applied to both atoms $`A`$ and $`B`$ (simultaneously) in the state Eq. (18). Then the new state of the pair at $`t=t_2`$ becomes

$`\mathrm{\Psi }(t_2)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}[e^{i\stackrel{}{K}(\stackrel{}{x}_{B2}\stackrel{}{x}_{B1})}|1_A|0_B`$ (34)
$``$ $`e^{i\stackrel{}{K}(\stackrel{}{x}_{A2}\stackrel{}{x}_{A1})}|0_A|1_B].`$ (35)
At any time $`t`$ after $`t_2`$ and before the next pulse, the state Eq. (20) evolves as

$`\mathrm{\Psi }(t)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}[e^{i\stackrel{}{K}(\stackrel{}{x}_{At}\stackrel{}{x}_{A2}+\stackrel{}{x}_{B2}\stackrel{}{x}_{B1})}|1_A|0_B`$ (36)
$``$ $`e^{i\stackrel{}{K}(\stackrel{}{x}_{Bt}\stackrel{}{x}_{B2}+\stackrel{}{x}_{A2}\stackrel{}{x}_{A1})}|0_A|1_B],`$ (38)
$`t:t_2t<t_3.`$
At a time $`t=t_3>t_2`$, a final, common $`H_{\frac{\pi }{2}}`$ pulse is applied to the entangled pair $`A`$ and $`B`$, and transforms the state Eq. (21) into

$`\mathrm{\Psi }(t_3)`$ $`=`$ $`{\displaystyle \frac{1}{2\sqrt{2}}}[e^{i(\mathrm{\Phi }_{}^{A}{}_{32}{}^{}+\mathrm{\Phi }_{}^{B}{}_{21}{}^{})}(|0_A|1_A)(|0_B+|1_B)`$ (39)
$``$ $`e^{i(\mathrm{\Phi }_{}^{B}{}_{32}{}^{}+\mathrm{\Phi }_{}^{A}{}_{21}{}^{})}(|0_A+|1_A)(|0_B|1_B)],`$ (40)
where, for $`i,j=1,2,3,\mathrm{}`$,

$$\mathrm{\Phi }_{}^{A}{}_{ij}{}^{}\stackrel{}{K}(\stackrel{}{x}_{Ai}\stackrel{}{x}_{Aj}),\mathrm{\Phi }_{}^{B}{}_{ij}{}^{}\stackrel{}{K}(\stackrel{}{x}_{Bi}\stackrel{}{x}_{Bj}).$$
(41)
After absorbing a factor of $`e^{i(\mathrm{\Phi }_{}^{B}{}_{32}{}^{}+\mathrm{\Phi }_{}^{A}{}_{21}{}^{})}`$ into the overall phase multiplying $`\mathrm{\Psi }(t_3)`$, Eq. (22) takes the more manageable form

$`\mathrm{\Psi }(t_3)`$ $`=`$ $`{\displaystyle \frac{1}{2\sqrt{2}}}\text{[}e^{i\mathrm{\Phi }}(|0_A|1_A)(|0_B+|1_B)`$ (42)
$``$ $`(|0_A+|1_A)(|0_B|1_B)\text{]},`$ (43)
with the phase observable $`\mathrm{\Phi }`$ given by

$`\mathrm{\Phi }`$ $``$ $`\mathrm{\Phi }_{}^{A}{}_{32}{}^{}+\mathrm{\Phi }_{}^{B}{}_{21}{}^{}\mathrm{\Phi }_{}^{B}{}_{32}{}^{}\mathrm{\Phi }_{}^{A}{}_{21}{}^{}`$ (44)
$`=`$ $`\stackrel{}{K}(\mathrm{\Delta }\stackrel{}{x}_32\mathrm{\Delta }\stackrel{}{x}_2+\mathrm{\Delta }\stackrel{}{x}_1),`$ (45)
where $`\mathrm{\Delta }\stackrel{}{x}\stackrel{}{x}_A\stackrel{}{x}_B`$ and $`\mathrm{\Delta }\stackrel{}{x}_i\stackrel{}{x}_{Ai}\stackrel{}{x}_{Bi}`$ for $`i,j=1,2,3,\mathrm{}`$ Collecting terms in Eq. (24),

$`\mathrm{\Psi }(t_3)`$ $`=`$ $`{\displaystyle \frac{1}{2\sqrt{2}}}\text{[}(e^{i\mathrm{\Phi }}1)|0_A|0_B+(e^{i\mathrm{\Phi }}+1)|0_A|1_B`$ (46)
$``$ $`(e^{i\mathrm{\Phi }}+1)|1_A|0_B(e^{i\mathrm{\Phi }}1)|1_A|1_B\text{]}.`$ (47)
The phase $`\mathrm{\Phi }`$ can now be observed by measuring the relative abundance in the state $`\mathrm{\Psi }(t_3)`$ of those pairs in which both atoms are in their excited states. \[Even though $`\mathrm{\Psi }(t_3)`$ is not stationary, its time evolution consists of pure phases multiplying each of the coefficients in Eq. (26), so the observation time is not critical.\] Specifically, according to Eq. (26), the fraction $`P_{11}`$ of those pairs in which both the $`A`$-atom and the $`B`$-atom are in the excited state $`|1`$ is given by

$$P_{11}=\frac{1}{8}\left|e^{i\mathrm{\Phi }}1\right|^2=\frac{1}{2}\mathrm{sin}^2\left(\frac{\mathrm{\Phi }}{2}\right).$$
(48)
On the other hand, if the Raman pulses making up the sequence $`H_\pi H_{\frac{\pi }{2}}`$ are aligned in time as before so that $`t_2=t_1+T`$ and $`t_3=t_2+T`$ for some common interrogation time $`T`$, we can write, using Eq. (25),

$$\mathrm{\Phi }=\stackrel{}{K}\mathrm{\Delta }\ddot{\stackrel{}{x}}(t_2)T^2+O(T^3),$$
(49)
where the magnitude of the $`O(T^3)`$ remainder is again down by a factor of order $`T/T_a`$ relative to the first (acceleration) term, with $`T_a`$ denoting the timescale over which the atoms’ acceleration varies. As long as the time scale $`T_a`$ is much larger than the interrogation time $`T`$, measurement of the single phase observable $`\mathrm{\Phi }`$ following the pulse sequence $`{}_{\overline{T}}{}^{}H_{\pi }^{}{}_{\overline{T}}{}^{}H_{\frac{\pi }{2}}^{}`$ on the entangled pair Eq. (17) yields a direct measurement of the inertial acceleration gradient $`\stackrel{}{K}\mathrm{\Delta }\stackrel{}{a}=\stackrel{}{K}(\stackrel{}{a}_A\stackrel{}{a}_B)`$ between the ensembles $`A`$ and $`B`$. Signal sensitivity of $`\mathrm{\Phi }`$ to the acceleration gradient $`\mathrm{\Delta }\stackrel{}{a}`$ can be enhanced using tricks similar to the one discussed in the note.

4. Producing the entangled atom pairs

A number of methods for entangling two-level atoms have been discussed in the literature and some implemented in experiments. For completeness, I will sketch two such methods which can be used to produce states of the form Eq. (17).

One method, originally proposed in Cabrillo et. al. , involves pumping two spatially separated atoms $`A`$ and $`B`$ into their excited states, so that the joint system is initially in the product state

$`|1_A|1_B.`$

A single-photon detector, which cannot (even in principle) distinguish the atom from which a detected photon arrives, is placed halfway between the atoms. Such a detection scheme can be implemented either by using a detector which is inherently insensitive to the direction of the photon that triggers it, or by using two detectors situated at the output ports of a symmetric beam-splitter on which the light paths from the two atoms are symmetrically impingent, thus erasing the which-path information. When one of the atoms makes a transition to its ground state, and the detector registers the emitted photon, the result of its measurement is to put the combined two-atom system into the entangled state

$$\frac{1}{\sqrt{2}}\left(|0_A|1_B+e^{i\varphi }|1_A|0_B\right),$$
(50)
where the phase $`\varphi `$ is a measure of the relative delay between the paths reaching the detectors from $`A`$ versus from $`B`$, and can be adjusted to have any desired value \[e.g. as in Eq. (17)\] by adjusting the position of the detector relative to the atoms.

A second method, recently investigated in detail by Haroche and coworkers is (in very rough outline) the following: Start with a single-mode cavity whose excitation frequency is tuned to $`\mathrm{\Omega }`$. Send the pair of atoms $`A`$ and $`B`$ into the cavity one after the other, with atom $`B`$ first. Initially, both atoms and the cavity are in their ground states:

$$|0_A|0_B|0_{EM},$$
(51)
where $`|0_{EM}`$ denotes the vacuum state of the cavity. After atom $`B`$ is in the cavity, apply a $`\pi /2`$-pulse on it, which transforms the state Eq. (30) into

$$\frac{1}{\sqrt{2}}|0_A\left(|0_B|1_{EM}|1_B|0_{EM}\right).$$
(52)
When both atoms are in the cavity, apply a second, $`\pi `$-pulse, this time on the atom $`A`$, thereby transforming the state Eq. (31) into

$$\frac{1}{\sqrt{2}}\left(|1_A|0_B|0_A|1_B\right)|0_{EM},$$
(53)
which, for the atom-pair $`A`$ and $`B`$, is in the desired form Eq. (17) up to an overall phase factor. Considering that a number of proof-of-principle experiments confirming it have already been performed successfully, this method (or, more precisely, a suitable variant which maintains the necessary correlation between atomic energy levels and linear momentum as discussed Sect. 2 above) is likely to be the method of choice in practical implementations of the entangled-atom gradiometry algorithm.

5. Quantifying the improvement

To what extent does the entangled-interferometry technique proposed in Sect. 3 represent an improvement over the standard, “classical” atom-wave gradiometry as reviewed in Sect. 2? To answer this question quantitatively, I will now carry out a simple analysis of the improvement in the signal-to-noise ratio expected with entangled-interferometry.

I will take the differential phase-shift as my observed signal, so in a gradient measurement with both classical and entangled interferometry, the signal is the dimensionless quantity

$$S=\stackrel{}{K}\mathrm{\Delta }\stackrel{}{a}T^2,$$
(54)
where $`\stackrel{}{a}\ddot{\stackrel{}{x}}`$, and other symbols are defined as in Eqs. (16) and (28), with the same approximation that allows neglecting terms of order $`O(T^3)`$ being assumed implicitly. In this discussion, for ease of notation I will label the two atomic ensembles by $`i=1,\mathrm{\hspace{0.33em}2}`$ as opposed to by $`A`$ and $`B`$ before (in Sects. 2—3). In a “classical” gradient measurement, two ensembles containing $`M`$ atoms each are used to obtain the two acceleration measurements, which are later differenced to obtain the signal

$$S=\mathrm{\Theta }_1\mathrm{\Theta }_2,$$
(55)
where the phase observables $`\mathrm{\Theta }_i`$ are related to the observed abundance (fluorescence) of excited-state atoms in each ensemble $`i=1,\mathrm{\hspace{0.33em}2}`$ via Eq. (14):

$$\mathrm{sin}^2\left(\frac{\mathrm{\Theta }_i}{2}\right)=\frac{M_i}{M}.$$
(56)
Here $`M_i`$ denotes the number of excited-state atoms at the end of the interrogation process in ensemble $`i`$, $`i=1,\mathrm{\hspace{0.33em}2}`$.

For comparison, imagine a similar experiment using the entangled algorithm of Sect. 3, involving $`M`$ pairs of entangled atoms interrogated in the pulse sequence of Eqs. (17)—(23) with the same interrogation-time constant $`T`$. At the end of that sequence, the number of atom pairs $`M_{11}`$ in which both atoms end up in their excited states is measured, which yields the signal

$$S=\mathrm{\Phi },$$
(57)
where according to Eq. (27)

$$\frac{1}{2}\mathrm{sin}^2\left(\frac{\mathrm{\Phi }}{2}\right)=\frac{M_{11}}{M}.$$
(58)
For simplicity, I will now make the additional assumption that the only sources of noise in both kinds of signal measurements (classical and entangled) are (i) phase noise, and (ii) shot noise due to the statistical fluctuations in the observed numbers of excited atoms: in $`M_1`$ and $`M_2`$ for the classical case, and in $`M_{11}`$ for the entangled case. I will assume a binomial distribution for these random variables with $`M`$ trials and with probabilities $`P_1`$, $`P_2`$, and $`P_{11}`$, respectively. The expected RMS shot noise in each number measurement $`M_i`$ is then given by

$$\mathrm{\Delta }M_i\sqrt{<M_{i}^{}{}_{}{}^{2}><M_i>^2}=\sqrt{P_i(1P_i)M},$$
(59)
and the expected RMS noise in $`M_{11}`$ is, similarly,

$$\mathrm{\Delta }M_{11}\sqrt{<M_{11}^{}{}_{}{}^{2}><M_{11}>^2}=\sqrt{P_{11}(1P_{11})M}.$$
(60)
In the classical measurement \[Eqs. (34)—(35)\], the shot noise contribution to each $`\mathrm{\Theta }_i`$ measurement can be written as

$$N_i=\frac{\mathrm{\Theta }_i}{M_i}\mathrm{\Delta }M_ii=1,\mathrm{\hspace{0.33em}2}.$$
(61)
Since the two measurements of $`\mathrm{\Theta }_i`$ are completely independent, the total variance (noise squared) in the quantity $`S=\mathrm{\Delta }\mathrm{\Theta }`$ is the sum of the variances in each $`\mathrm{\Theta }_i`$:

$$N_c=\left(N_{1}^{}{}_{}{}^{2}+N_{2}^{}{}_{}{}^{2}\right)^{\frac{1}{2}}\sqrt{2}N_1=\sqrt{2}\frac{\mathrm{\Theta }_1}{M_1}\mathrm{\Delta }M_1,$$
(62)
where the approximate equality holds because of the assumption that the signal $`S1`$. Substituting Eqs. (35), (38) and (14) in Eq. (41) gives

$$N_c=\left(\frac{2}{M}\right)^{\frac{1}{2}}.$$
(63)
In the entangled measurement, according to Eqs. (36)—(37)

$$S=\mathrm{\Phi }=2\mathrm{arcsin}\sqrt{\frac{2M_{11}}{M}}2\sqrt{2}\sqrt{\frac{M_{11}}{M}},$$
(64)
and the shot noise contribution is

$$N_e=\frac{\mathrm{\Phi }}{M_{11}}\mathrm{\Delta }M_{11}=\left(\frac{2(1p_{11})}{M}\right)^{\frac{1}{2}}\left(\frac{2}{M}\right)^{\frac{1}{2}},$$
(65)
where the second equality follows from Eq. (39), and the third by the assumption that the signal $`S1`$. Remarkably, the shot noise contributions in the classical and entangled measurements are identical. Since a typical classical gradient measurement operates in a regime where the number of atoms $`M`$ in each ensemble $`i`$ is large enough to make shot noise negligible, the same operating conditions will make it negligible also in the entangled case, and therefore phase noise is the dominant source of error in both kinds of measurement.

To compare phase-noise performance of the two kinds of measurement, I will model the noise as simply proportional to the magnitude of the fluorescence signal:

$$<\mathrm{\Theta }_{i}^{}{}_{}{}^{2}><\mathrm{\Theta }_i>^2=\sigma ^2\mathrm{sin}^2\frac{<\mathrm{\Theta }_i>}{2},i=1,\mathrm{\hspace{0.33em}2}$$
(66)
for the classical case, and

$$<\mathrm{\Phi }^2><\mathrm{\Phi }>^2=\sigma ^2\frac{1}{2}\mathrm{sin}^2\frac{<\mathrm{\Phi }>}{2},i=1,\mathrm{\hspace{0.33em}2}$$
(67)
in the entangled case, where $`\sigma `$ is a constant. By the same arguments as above \[Eqs. (40)—(41)\], the phase noise in the differential acceleration signal is, for the classical case

$$N_c\sqrt{2}\sigma \mathrm{sin}\frac{<\mathrm{\Theta }_1>}{2},$$
(68)
and for the entangled case

$$N_e=\frac{1}{\sqrt{2}}\sigma \mathrm{sin}\frac{<\mathrm{\Phi }>}{2}.$$
(69)
Therefore the ratio of phase noise in the entangled measurement to that in the classical measurement can be written in the form \[cf. Eqs. (36), (34), and (33)\]

$$\frac{N_e}{N_c}=\frac{1}{2}\frac{\mathrm{sin}\frac{<\mathrm{\Phi }>}{2}}{\mathrm{sin}\frac{<\mathrm{\Theta }_1>}{2}}\frac{1}{2}\frac{<\mathrm{\Phi }>}{<\mathrm{\Theta }_1>}\frac{1}{2}\frac{\mathrm{\Delta }\stackrel{}{a}}{\stackrel{}{a}}.$$
(70)
For a typical measurement performed on earth with a 1m separation between the two atomic ensembles, the earth’s $`\mathrm{\hspace{0.33em}3000}`$E gradient field (1E$`=10^9\mathrm{sec}^2`$) and 1g gravity suggests suppression of phase noise via entanglement by a factor of up to $`1.5\times 10^7`$, or more than six orders of magnitude. For measurements performed in spacecraft in low earth orbit, the ambient acceleration $`\stackrel{}{a}`$ will be due to non-isolated common-mode vibrations originating from drag forces acting on the spacecraft. Therefore, for a gradiometer in low earth orbit (common-mode accelerations due to atmospheric drag being typically of order $`10^5`$g or higher) the phase-noise suppression $`N_e/N_c`$ via entangled interferometry can be expected to be a reduced but still significant two orders of magnitude or better.

6. Future work

The Entangled Quantum Interferometry algorithm as presented here can be generalized to design sensors of other “distributed” quantities (such as higher derivatives of the gravitational curvature tensor); these generalizations will be discussed in a future publication.

###### Acknowledgements.

I would like to acknowledge valuable discussions with Lute Maleki and Nan Yu of the JPL Frequency Standards Laboratory, and with John Clauser, Jonathan Dowling, Dmitry Strekalov, and Colin Williams of the JPL Quantum Computing Technologies group. The research described in this paper was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration, and was supported by a contract with the National Reconnaissance Office.


text:
# Strong “quantum” chaos in the global ballooning mode spectrum of three-dimensional plasmas

## Abstract

The spectrum of ideal magnetohydrodynamic (MHD) pressure-driven (ballooning) modes in strongly nonaxisymmetric toroidal systems is difficult to analyze numerically owing to the singular nature of ideal MHD caused by lack of an inherent scale length. In this paper, ideal MHD is regularized by using a $`k`$-space cutoff, making the ray tracing for the WKB ballooning formalism a chaotic Hamiltonian billiard problem. The minimum width of the toroidal Fourier spectrum needed for resolving toroidally localized ballooning modes with a global eigenvalue code is estimated from the Weyl formula. This phase-space-volume estimation method is applied to two stellarator cases.

In design studies for new magnetic confinement devices for fusion plasma experiments (e.g. investigations  leading to the proposed National Compact Stellarator Experiment, NCSX ), the maximum pressure that can stably be confined in any proposed magnetic field configuration is routinely estimated by treating the plasma as an ideal magnetohydrodynamic (MHD) fluid. One linearizes about a sequence of equilibrium states with increasing pressure, and studies the spectrum of normal modes (frequency $`\omega `$) to determine when there is a component with $`\mathrm{Im}\omega >0`$, signifying instability.

Even with the simplification obtained by using the ideal MHD model, the computational task of determining the theoretical stability of a three-dimensional (i.e. nonaxisymmetric) device, such as NCSX or the four currently operating helical axis stellators , remains a challenging one.

The problem can be posed as a Lagrangian field theory, with the potential term being the energy functional $`\delta W`$ . For a static equilibrium, the kinetic energy is quadratic in $`\omega `$, so that $`\omega ^2`$ is real. Thus instability occurs when $`\omega ^2<0`$. There are two main approaches to analyzing the spectrum—local and global.

In the *local* approach, which is used for analytical simplification, one orders the scale length of variation of the eigenfunction across the magnetic field lines to be short compared with equilibrium scale lengths . Both interchange and ballooning stability can be treated by solving the general *ballooning equations* , a system of ordinary differential equations defined on a given magnetic field line.

The *global* (Galerkin) approach is to expand the plasma displacement field in a finite basis set, inserting this ansatz in the Lagrangian to find a matrix eigenvalue representation of the spectral problem. This approach has been implemented for ideal MHD in three-dimensional plasmas in two codes, TERPSICHORE  and CAS3D .

Although the Galerkin approach is potentially exact, if one could use a complete, infinite basis set, it is in practice computationally challenging due to the large number of basis functions required to resolve localized instabilities. This leads to very large matrices which must be diagonalized by iterative methods. There is a need for analytical insight to determine a suitable truncated basis set and to predict the nature of the spectrum, e.g. whether it is continuous or discrete.

Such insight may be obtained by a *hybrid* local-global approach, in which one uses a Wentzel–Kramers–Brillouin (WKB) representation of the eigenfunction. In the short-wavelength limit, the same analytical simplifications as are obtained in the local approach are found to give a local dispersion relation that can be used to give information on the global spectrum by using ray tracing and semiclassical quantization.

In axisymmetric systems  or in cases where helical ripple can be averaged out, giving an adiabatic invariant, , the ray equations are integrable and hence the spectrum is characterized by “good quantum numbers”.

However, it has been known for many years  that the ray-tracing problem in strongly three-dimensional systems is singular because, in the absence of an adiabatic invariant, the phase-space motion is not bounded—the rays escape to infinity in the wavevector sector. Dewar and Glasser  argued that this gives rise to a *continuous* unstable spectrum, with correspondingly singular generalized eigenfunctions. (A more rigorous treatment involves the concept of the essential spectrum and Weyl sequences .)

Our proposed regularization of this singularity can be understood using a simple quantum analogy. Consider the one-dimensional time-independent Schrödinger equation $`H\psi =E\psi `$ in the limit as the mass of the particle goes to infinity. Then the kinetic energy disappears and the Hamiltonian becomes $`H=V(x)`$, where $`V`$ is the potential energy, assumed here to be the harmonic oscillator potential, $`\frac{1}{2}x^2`$ in suitable units. In the usual Hilbert space the energy spectrum is continuous: $`E0`$ and the (generalized) eigenfunctions singular: $`\psi (x)=\delta (xx_\mathrm{E})\pm \delta (x+x_\mathrm{E})`$, where $`V(x_\mathrm{E})E`$.

We now seek a regularization of this problem by restricting $`\psi `$ to the space of functions with a finite bandwidth in wavenumber $`k`$:

$$\psi (x)=_{k_{\mathrm{max}}}^{k_{\mathrm{max}}}\frac{dk}{2\pi }\psi _k\mathrm{exp}ikx.$$
(1)
This truncated Fourier-integral representation models what occurs when one seeks to find the spectrum numerically using a truncated Fourier-series representation.

We take as starting point a Lagrangian for the wavefunction,

$$L=_{\mathrm{}}^{\mathrm{}}\psi ^{}[EV(x)]\psi 𝑑x.$$
(2)
Inserting Eq. (1) in Eq. (2) gives

$`L`$ $`=`$ $`{\displaystyle _{k_{\mathrm{max}}0}^{k_{\mathrm{max}}+0}}[E|\psi _k|^2|{\displaystyle \frac{d\psi _k}{dk}}+\psi _k\delta (k+k_{\mathrm{max}})`$ (4)
$`\psi _k\delta (kk_{\mathrm{max}})|^2]{\displaystyle \frac{dk}{2\pi }}.`$
This is infinite unless we require the coefficients of the $`\delta `$-functions to vanish. That is, $`\psi _k=0`$ at $`k=\pm k_{\mathrm{max}}`$. The Euler–Lagrange equation is $`(d^2/2dk^2+E)\psi _k=0`$, which has the solutions $`\mathrm{exp}\pm i(2E)^{1/2}k`$. These waves would propagate to infinity if it were not for the reflecting boundary conditions at $`\pm k_{\mathrm{max}}`$ we have just derived.

That is, we have removed the continuum by box quantization in $`k`$-space. In the following we shall do the same for the ballooning mode problem.

As in  we write the magnetic field of an arbitrary three-dimensional toroidal equilibrium plasma with nested magnetic flux surfaces labeled by an arbitrary parameter $`s`$ as $`𝐁=\mathbf{}\zeta \mathbf{\times }\mathbf{}\psi q\mathbf{}\theta \mathbf{\times }\mathbf{}\psi \mathbf{}\alpha \mathbf{\times }\mathbf{}\psi `$, where $`\alpha \zeta q\theta `$. Here, $`\theta `$ and $`\zeta `$ are the poloidal and toroidal angles, respectively, $`\psi (s)`$ is the poloidal flux function, and $`q(s)`$ is the inverse of the rotational transform. Since $`𝐁\mathbf{}\mathbf{}s=𝐁\mathbf{}\mathbf{}\alpha =0`$, $`s`$ and $`\alpha `$ serve to label an individual field line.

We take the stream function  to be given by $`\phi =\widehat{\phi }\mathrm{exp}(iSi\omega t),`$ where $`\widehat{\phi }(\theta |s,\alpha )`$ is assumed to vary on the equilibrium scale. The phase variation is taken to be rapid, so $`𝐤\mathbf{}S`$ is ordered to be large. The frequency $`\omega `$ is ordered $`O(1)`$, which requires that the wave vector be perpendicular to $`𝐁`$: $`𝐤\mathbf{}𝐁0`$. (In this study we consider unstable ideal MHD modes, $`\omega ^2<0`$.)

It immediately follows that the eikonal is constant on each field line: $`S=S(\alpha ,s)`$. From the definition of the wave vector, $`𝐤=k_\alpha \mathbf{}\alpha +k_s\mathbf{}sk_\alpha [\mathbf{}\alpha +\theta _kq^{}(s)\mathbf{}s]`$ where $`k_\alpha S/\alpha `$ and $`k_sS/s`$. Here the anglelike ballooning parameter $`\theta _k`$ appears naturally as the ratio $`k_s/q^{}(s)k_\alpha `$ .

The ballooning equation emerges in the large $`|𝐤|`$ expansion  as an ordinary differential equation to be solved on each field line $`(\alpha ,s)`$ with given $`(k_\alpha ,k_s)`$ under the boundary condition $`\widehat{\phi }(\theta )0`$ at infinity to give the eigenvalue $`\lambda (\alpha ,s,k_\alpha ,k_s)`$. This constitutes a local dispersion relation $`\lambda \rho \omega ^2`$ (the mass density $`\rho `$ being assumed constant everywhere).

The ray equations are the characteristics of the eikonal equation $`\lambda (\alpha ,s,_\alpha S,_sS)=\rho \omega ^2`$. These are Hamiltonian equations of motion with $`\alpha ,s`$ the generalized coordinates, $`k_\alpha ,k_s`$ the canonically conjugate momenta, and $`\lambda `$ as the Hamiltonian.

In axi- or helically symmetric systems all field lines on a given magnetic surface are equivalent—$`\alpha `$ is ignorable and $`k_\alpha `$ is a constant of the motion. In this case the equations are integrable and semiclassical quantization can be used to predict the approximate spectrum of global ballooning instabilities . This technique can sometimes be applied successfully, even in nonsymmetric systems, if there are regions of phase space with a large measure of invariant tori . In  this was verified using the global eigenvalue code TERPSICHORE .

At the other extreme, if the ray orbits are chaotic (but still bounded) then the global spectrum is not regularly structured, but must rather be described statistically by the density of states and the probability distribution of level spacings using the techniques of quantum chaos theory (see e.g. ).

However, because of the scale invariance of the ideal MHD equations, $`\lambda `$ depends only on the direction of $`𝐤`$, *not* on its magnitude: $`\lambda =\lambda (\alpha ,s,\theta _k)`$. This has the consequence that the ray orbits are *unbounded* in phase space, so, strictly speaking, ideal MHD gives rise to a quantum *chaotic scattering*  problem rather than a straight quantum *chaos* problem. This leads to the continuous spectrum  with singular generalized eigenfunctions that cannot really be represented using the simple eikonal ansatz.

On the other hand, the absence of a natural length scale in ideal MHD is a mathematical artifact. Physically, the ion Larmor radius provides a lower cutoff in space, or an upper cutoff in $`|𝐤|`$, beyond which ideal MHD ceases to apply. The ballooning equation is also physically regularized by inclusion of diamagnetic drift .

However, since in general it leads to a complex ray tracing problem , we shall not attempt to model diamagnetic drift stabilization in this paper. Rather, we regularize the ray equations simply by adding a barrier term to the effective ray “Hamiltonian” $`H(\alpha ,s,k_\alpha ,k_s)`$,

$$H=\lambda (\alpha ,s,k_\alpha ,k_s)+U(k_\alpha ),$$
(5)
where the barrier potential we use is $`U(k_\alpha )K(|k_\alpha |k_{\mathrm{max}})^2`$ for $`|k_\alpha |>k_{\mathrm{max}}`$ and 0 for $`|k_\alpha |<k_{\mathrm{max}}`$. In the limit of the constant $`K\mathrm{}`$, this infinite box potential gives the ideal MHD ray equations for $`|k_\alpha |<k_{\mathrm{max}}`$ and reflecting boundary conditions at $`|k_\alpha |=k_{\mathrm{max}}`$. Thus we have a two-degree of freedom Hamiltonian billiard problem.

Although overly crude for modeling FLR stabilization, the cutoff at $`|k_\alpha |=k_{\mathrm{max}}`$ provides a reasonable model for representing the finite spectral bandwidth in the toroidal Fourier mode number ($`n`$) representation used in the global eigenvalue codes TERPSICHORE  and CAS3D .

Using ballooning-unstable plasma equilibria calculated for the H-1NF heliac  using the VMEC code , detailed parameter scans have been undertaken for two cases. The first case studied  was obtained by increasing the pressure gradient of a marginally stable equilibrium  uniformly across the plasma and thus was ballooning unstable at the edge of the plasma. The ray tracing problem for this case would involve consideration of the effect of the plasma boundary.

Thus a second equilibrium, ballooning stable near the edge of the plasma, was calculated for the purposes of the present paper. This case has a more peaked pressure profile than the first, but both have average $`\beta 1\%`$, where $`\beta `$ is the ratio of plasma pressure to magnetic field pressure.

The $`q`$-profiles are not monotonic—in the peaked pressure profile case studied in this paper, $`q`$ was 0.8895 on the magnetic axis, rising to a maximum value of 0.8964 quite close to the magnetic axis, then falling monotonically to 0.8675. Clearly the (global) magnetic shear is very weak. Despite this fact and the non-monotonicity, there is some formal simplification in choosing $`sq`$, and we have taken $`s=q`$ since the region of plasma studied is in a monotonic-decreasing part of the $`q`$-profile (the decreasing region outside the maximum-$`q`$ surface).

In these scans the most unstable ballooning eigenvalue was tabulated on a three-dimensional grid in $`s,\alpha ,\theta _k`$ space. The dependence on $`\alpha `$ was found to be rapid. The dependence on $`\theta _k`$ was much slower, but the variation was sufficient that the higher-growth-rate isosurfaces formed a set of distinct, topologically spherical branches. It was argued in  that this branch structure is produced by Anderson localization in bad curvature regions due to the strong breaking of both helical and axisymmetry in H-1NF.

According to the perturbation expansion in $`q^{}`$ described in , a quadratic form in $`\alpha ,\theta _k`$ should form a good approximation to $`\lambda \lambda _{\mathrm{min}}(q)`$ in the neighborhood of the central branch. Accordingly a least-squares fit on each surface was performed to provide a simple analytical description of the $`(0,0)`$  branch.

The radial dependence of the fitting coefficients was approximated by fitting to third-degree polynomials in $`q`$. Sections of the resulting approximation to the central branch are shown in Fig. 1. The isosurface spans a substantial range of magnetic surfaces within the plasma — the narrow range of variation in $`q`$ is due to the low magnetic shear in H-1NF.

In order to establish the nature of the ray dynamics described by the regularized Hamiltonian, Eq. (5), a numerical integration with cutoff at $`k_{\mathrm{max}}=50`$ was performed with initial conditions $`q=q_2`$, $`\alpha =0`$, and $`k_\alpha =5`$, where $`[q_1,q_2]=[0.8852,0.8951]`$ is the $`q`$-range spanned by the $`\lambda =6`$ isosurface as seen in Fig. 1. (A run with $`k_\alpha =10`$ was also performed, with similar results.) Choosing the value $`K=1`$ gave a good compromise between the sharp boundary potential to be modeled, and the smooth potential required for the numerical integration. The orbit remained on the “energy shell” $`\lambda =6`$ to within an accuracy of one part in $`10^6`$ over the “time” interval of the integration, $`7500`$.

The two Poincaré plots in Fig. 2 show the orbit to be strongly chaotic, filling the “energy shell” ergodically, except that the regions $`k_\alpha >0`$ and $`k_\alpha <0`$ are dynamically disjoint. The solid curve shown surrounding the outer limits of the “energetically accessible” region is calculated by solving $`\lambda (0,q,k_q/k_{\mathrm{max}})=6`$.

According to the Weyl formula \[16, pp. 257–261\], the number, $`N(\lambda _{\mathrm{max}})`$, of global eigenmodes with eigenvalues below the eigenvalue $`\lambda _{\mathrm{max}}`$ is given, asymptotically in the limit $`N\mathrm{}`$, as $`N(\lambda _{\mathrm{max}})v_{4\mathrm{D}}(\lambda _{\mathrm{max}})/(2\pi )^2`$. Here $`v_{4\mathrm{D}}(\lambda _{\mathrm{max}})`$ is the volume of the dynamically acessible 4-dimensional phase-space region $`\lambda (\alpha ,q,k_q/k_\alpha )<\lambda _{\mathrm{max}}`$, $`0<k_\alpha <k_{\mathrm{max}}`$. The $`k_\alpha `$ integration can be performed analytically, giving $`v_{4\mathrm{D}}(\lambda _{\mathrm{max}})=\frac{1}{2}k_{\mathrm{max}}^2v_{3\mathrm{D}}(\lambda _{\mathrm{max}}),`$ where $`v_{3\mathrm{D}}(\lambda _{\mathrm{max}})`$ is the volume within the isosurface $`\lambda (\alpha ,q,\theta _k)=\lambda _{\mathrm{max}}`$. Thus

$$N(\lambda _{\mathrm{max}})\frac{1}{8\pi ^2}k_{\mathrm{max}}^2v_{3\mathrm{D}}(\lambda _{\mathrm{max}}).$$
(6)
We can make a rather rough estimate of the minimum value of $`n_{\mathrm{max}}`$ required for CAS3D or TERPSICHORE to find even one eigenvalue with $`\lambda <\lambda _{\mathrm{max}}`$ by setting $`N(\lambda _{\mathrm{max}})=1`$ and calculating $`k_{\mathrm{max}}n_{\mathrm{max}}`$ from Eq. (6). This gives $`n_{\mathrm{max}}(N=1)(8\pi ^2/v_{3\mathrm{D}})^{1/2}`$.

The isosurface $`\lambda =6`$ studied above is about the largest of the disjoint topologically spherical isosurfaces corresponding to the highly toroidally localized strongly ballooning unstable regions of $`\alpha ,q,\theta _k`$ space. (For $`\lambda >6`$ the isosurfaces are no longer topologically spherical.) Using the polynomial fits described above, we calculate $`v_{3\mathrm{D}}(6)=0.02158`$. This gives $`n_{\mathrm{max}}(N=1)60`$. Assuming that the dominant contributions to the MHD energy $`\delta W`$ come from the rational surfaces intersecting the $`\lambda =6`$ isosurface, we thus predict that it would be necessary to include, as a minimum set, basis functions corresponding to one of the two “mode families”  contained in the set $`(n,m)=(9,8)`$, $`(18,16)`$, $`(19,17)`$, $`(27,24)`$, $`(28,25)`$, $`(35,31)`$, $`(36,32)`$, $`(37,33)`$, $`(38,34)`$, $`(44,39)`$, $`(45,40)`$, $`(46,41)`$, $`(47,42)`$, $`(53,47)`$, $`(54,48)`$, $`(55,49,)`$, $`(56,50)`$, and $`(57,51)`$ to resolve a toroidally localized ballooning mode. (Here $`n,m`$ are the toroidal and poloidal Fourier mode numbers, respectively.)

The large value of $`n_{\mathrm{max}}(N=1)`$ required, and the unusual spread in $`n`$ required in the basis set, will make these modes difficult to resolve using global eigenvalue codes (e.g. the simplifying phase factor method sometimes used in CAS3D studies  would not be appropriate). It is hoped that the Weyl formula estimate above will act as a guide in a future more extensive study using such a code. Physically, the large value of $`n_{\mathrm{max}}`$ suggests that toroidally localized ballooning modes in H-1NF should be subject to strong FLR stabilization.

We can also apply the same approach to the toroidally localized ballooning branches found in the Large Helical Device (LHD) study . From the plots in  we estimate $`v_{3\mathrm{D}}0.05`$, which gives $`n_{\mathrm{max}}(N=1)40`$.

The ballooning calculations were carried out on the Australian National University Supercomputer Facility’s Fujitsu VPP300 vector processor. We thank Dr. H. J. Gardner for providing the H-1 heliac VMEC input files and Dr. S. P. Hirshman for use of the VMEC equilibrium code. Some of this work was done while one of us (RLD) was a visiting scientist at Princeton University Plasma Physics Laboratory, supported under US DOE contract No. DE-AC02-76CH0-3703. Useful conversations with Drs. M. Redi and A.H. Boozer are gratefully acknowledged.


text:
# NICS–TNG infrared spectroscopy of NGC1068: the first extragalactic measurement of [P II] and a new tool to constrain the origin of [Fe II] line emission in galaxies.

## 1 Introduction

Since the first infrared spectroscopic observations of galaxies and supernova remnants in the 80’s, the emission lines of \[Fe II\] have become a popular and debated issue (Moorwood & Oliva moorwood (1988), Forbes & Ward forbes93 (1993), Simpson et al. simpson96 (1996), Veilleux et al. veilleux97 (1997), Alonso–Herrero et al. alonso97 (1997), Mouri et al. mouri00 (2000)).

From the observational point of view, \[Fe II\] is weak in HII regions and planetary nebulae while extremely strong in shock–excited filaments of supernova remnants. Since relatively bright \[Fe II\] emission is commonly found in the IR spectra of normal and active galaxies, many authors have considered the possibility of using this line as shock tracer and, even, to count the number of supernova remnants (e.g. Colina colina93 (1993), Vanzi & Rieke vanzi97 (1997), Engelbracht et al. engelbracht98 (1998)).

From the theoretical point of view, a low density region with normal abundances can become a strong source of \[Fe II\] only if the following conditions are satisfied.

* Most of iron must be in the gas phase, i.e. dust grains must have been destroyed.
* The gas electron temperature must be large enough to collisionally excite the upper levels of the lines, in practice $`T_e5000`$ K.
* Most of iron must be singly ionized, i.e. Fe+/Fe must be close to unity.

Given the low ionization potential of Fe+ and the high efficiency of the Fe++ \+ Ho charge–exchange recombination reactions, the latter two conditions are equivalent to saying that bright \[Fe II\] lines can only be formed in regions where hydrogen is partly ionized (e.g. Oliva et al. oliva\_snr1 (1989)).

The most efficient mechanisms for creating extended regions of hot, partially ionized gas are shocks and photoionization by soft X–rays. In both cases the volume emission measure of the partially ionized region could easily exceed that of the fully ionized gas. The only important difference between the two mechanisms is that photoionization is unable to destroy the toughest iron–based grains which are otherwise easily sputtered by shock fronts.

In practice, therefore, a purely photoionized region emitting \[Fe II\] can be easily recognized by measuring the abundance of Fe+ relative to any non-refractory species which forms in the same partially ionized region. Finding a very low iron abundance would unequivocally imply that the gas has not been significantly processed by shocks.

This could be in principle obtained by comparing the infrared \[Fe II\] and optical \[O I\] lines but, in practice, the large difference in critical densities, reddening and the problems of comparing data taken with different instruments makes it impossible to obtain a clear–cut conclusion (see e.g. Mouri et al. mouri93 (1993), Alonso–Herrero et al. alonso97 (1997), Larkin et al. larkin98 (1998)). A much more reliable determination of the iron relative abundance should be derived from lines close in wavelengths and with similar critical densities such as the \[P II\] and \[Fe II\] lines which are discussed in this Letter.

## 2 Observations and results

The data were collected at the Telescopio Nazionale Galileo (TNG) in November 2000 during the commissioning phase of NICS, the near infrared camera and spectrometer expressly designed and built for this telescope. This instrument is a FOSC–type cryogenic focal reducer equipped with two interchangeable cameras feeding a Rockwell Hawaii 10242 array. The camera used for the spectroscopic observations has a projected scale of 0.25”/pixel (Oliva & Gennari gennari1 (1995), Baffa et al. baffa (2000)). The spectroscopic modes are achieved by means of a series of glass–resin grisms which can be inserted in the 22 mm collimated beam (Vitali et al. vitali (1997)). The spectrum of NGC1068 was collected through a slit of 0.75” (=3 pixels) width and using the IJ grism which yields a 0.89-1.46 $`\mu `$m spectrum with a dispersion of 5.7 Å/pix. The slit was oriented N–S (i.e. at PA=0) and centered on the 1 $`\mu `$m continuum peak. The acquisition consisted of a series of four 5–minute exposures with the object set at different positions along the slit followed by halogen flats. The atmospheric spectral response and the instrumental efficiency were determined using spectra of the O6.5V star HD42088 whose intrinsic spectrum was approximated by $`F_\lambda `$=1.7$`10^9`$ ($`\lambda `$/1.25)-3.7 erg cm-2 s-1 $`\mu `$m-1.

The spectrum of the central 0.75”x2” region is displayed in Fig. 1 where we also show, for comparison, spectra of the Orion Bar (Walmsley et al. 2000) and unpublished spectra of supernova remnants collected in 1992 using IRSPEC at the ESO-NTT telescope. The relative line fluxes are summarized in Table 1. Evident is the difference between the very large \[Fe II\]/\[P II\] ratio measured in SNR’s and the much smaller values found in the Orion Bar and in NGC1068.

In principle, the emission feature peaking at 1.188 $`\mu `$m could be contaminated by \[Ni II\] 1.191 $`\mu `$m. However, this line was measured in the Crab nebula at a level of only $``$15% of \[Fe II\] (Rudy et al. rudy (1994)) and was not detected in the supernova remnants listed in Table 1. Moreover, \[Ni II\] has a critical density of only $``$500 cm-3 (Nusbbaumer & Storey nussbaumer2 (1982)), two orders of magnitude lower than the critical density of \[Fe II\] and much lower than the electron density pertinent to the region of NGC1068 under consideration here. For these reasons we believe that \[Ni II\] does not affect significantly the measurement of \[P II\].

It is also interesting to note that our data confirm the identification of \[SIX\] by Marconi et al. (marconi96 (1996)) and include the first detection of \[Fe XIII\] in an extragalactic object. The latter identification is also supported by measurements of the higher ionization green line of \[Fe XIV\] (Kraemer & Crenshaw kraemer (2000)).

## 3 Discussion

### 3.1  \[Fe II\], \[P II\] and the Fe/P abundance ratio

The near–IR lines of \[P II\] and \[Fe II\] have several interesting similarities. They lie nearby in wavelength, have similar excitation temperatures and critical densities and their parent ions have similar ionization potentials and radiative recombination coefficients. Using the collision strengths of Krueger & Czyzak (krueger (1970)), Zhang & Pradhan (zhang (1995)) and the transition probabilities of Mendoza & Zeippen (mendoza (1982)), Nussbaumer & Storey (nussbaumer (1988)) yields

$$\frac{n(\text{Fe}\text{+})}{n(\text{P}\text{+})}2\frac{I(\text{[Fe II]}1.257\mu \text{m})}{I(\text{[P II]}1.188\mu \text{m})}$$
$`(1)`$
which is accurate within a factor of 2 for all the temperatures and densities of interest. The only important difference between the two species is that charge exchange recombinations between P++ and neutral hydrogen are $``$2 orders of magnitude less efficient than the Fe++ \+ Ho reactions (Kingdon & Ferland kingdon96 (1996)). This implies

$$\frac{n(\text{Fe}\text{+})}{n(\text{P}\text{+})}\frac{n(\mathrm{Fe})}{n(\mathrm{P})}$$
$`(2)`$
which combined with Eq. (1) yields

$$\frac{n(\mathrm{Fe})}{n(\mathrm{P})}2\frac{I(\text{[Fe II]}1.257\mu \text{m})}{I(\text{[P II]}1.188\mu \text{m})}$$
$`(3)`$
i.e. the Fe/P abundance ratio is quite well constrained by the \[Fe II\]/\[P II\] ratio and, if anything, it is overestimated. Finally, it is interesting to note that, for a solar Fe/P$``$100 abundance ratio, one expects \[Fe II\]/\[P II\]$``$50 i.e. a ratio similar to that measured in supernova remnants.

### 3.2  Fe/P abundance ratio and dust destruction

Iron is a well known refractory species whose gas phase abundance in the ISM is often found to be down by many orders of magnitudes relative its cosmic value. The only regions where iron is not found to be significantly depleted are those associated with fast ($``$100 km/s) shocks which can effectively destroy the grains by sputtering. In normal photoionized regions the depletion of iron ranges between the factor of $``$0.1 measured in the Orion Bar (Baldwin et al. baldwin (1996)) to significantly lower values found in planetary nebulae (e.g. Oliva et al. oliva6302 (1996), Perinotto et al. perinotto (1999)). These differences probably reflect the fact that a variable (though small) fraction of iron is locked into relatively soft grains which can be easily destroyed without the need of fast shocks.

Phosphorus is a non-refractory species whose measured depletion in ionized gas is close to unity. Therefore, the Fe/P relative abundance should give a direct estimate of the iron depletion or, equivalently, of the presence of robust dust in a given region. This is indeed confirmed by the data presented here (Fig. 1 and Table 1) which show variations by more than one order of magnitude between the large, quasi–solar Fe/P ratio found in SNR’s and the much smaller values found in the other objects.

### 3.3  \[Fe II\]/\[P II\] and the origin of \[Fe II\] in galaxies

As already discussed in the introduction, determining the origin of \[Fe II\] line emission is of crucial importance for any program aiming at using \[Fe II\] for tracing shock fronts and/or constraining the supernovae rate in galaxies. The results obtained here indicate that the \[Fe II\]/\[P II\] ratio could provide a clear–cut answer to this problem. The line ratio is large ($``$20) in regions excited by fast shocks while low ($``$2) in normal photoionized region and in NGC1068. This indicate that shocks are not the dominant source of \[Fe II\] in the central $``$200 pc of NGC1068 where, most likely, the lines are produced by photoionization from the active nucleus, as already indicated by detailed photoionization modeling (e.g. Kraemer & Crenshaw kraemer (2000)). The copious flux of soft X–rays from the AGN creates a very extended partially ionized region which is responsible for the strong emission of low ionization species such as \[S II\], \[O I\], \[P II\] and \[Fe II\]. However, the latter is relatively weak because most of the iron is locked into dust grains. The relative intensities of the low and higher ionization lines depend on a complex combination of the ionization parameter, density and of the spectral shape of the ionizing radiation. Nevertheless, the ratio between “close relatives” such as \[Fe II\]/\[P II\] are little influenced by this and, as discussed in Sect. 3.1, almost solely depend on the Fe/P relative abundance which, in turn, is a direct measurement of iron depletion (Sect. 3.2).

## 4 Conclusions

Given the above results and considerations, we propose using the ratio between \[Fe II\] (1.257 $`\mu `$m) and \[P II\] (1.188 $`\mu `$m) as a general tool for constraining the origin of \[Fe II\] line emission in galaxies. The diagnostic works as follows

* In objects with low \[Fe II\]/\[P II\] ratios shocks do not play an important role in the lines excitation.
* Large values of \[Fe II\]/\[P II\] ($``$20) indicate that the emitting gas has recently passed through a fast shock which sputtered and destroyed most of the dust grains. It is therefore likely that the lines are also produced by shock excitation

###### Acknowledgements.

This paper is based on observations made with the Italian Telescopio Nazionale Galileo (TNG) operated on the island of La Palma by the Centro Galileo Galilei of the CNAA (Consorzio Nazionale per l’Astronomia e l’Astrofisica) at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias. We are grateful to all the technical staff and telescope operators for their assistance during the commissioning phase of NICS.