text:
# The Fermion Boson Interaction Within the Linear Sigma Model at Finite Temperature

## I Introduction

The fermion dispersion relation to one-loop order at high temperature has been studied several times in the literature . The dispersion relation of a fermion, which gives its energy $`\omega `$ as it propagate through the medium as a function of its momentum $`k`$, is very important in different kinds of physical situations  and has been focused on various approximations and limits . During the last few years, there has been some controversy if the damping rate is gauge dependent or not and other problems . It has been proposed that a proper resummation cures such problems found in most one-loop calculations . In fact, if one wants consistent results of damping rates at positive temperatures, it is crucial the use of methods like the well known hard thermal loop (HTL) resummation of Braaten and Pisarski  which employs resumed propagators and vertices and keep the correct order of the coupling constant. Whith this method, Braaten and Pisarski has given in ref.  a formal proof that resummation produces gauge-invariant results for the damping rates of quarks and gluons. In this article we reinvestigate the self-energy of massless fermions interacting with massless bosons at high temperature in the framework of the linear sigma model. This is the simpler, but instructive situation of fermions interacting with scalar and pseudoscalar fields. These are similar to Yukawa couplings which has been discussed by M. Thoma in ref. . Using the HTL, Thoma has found a kinematic restriction (not observed in the case of quarks or gluons) depending on the effective fermion and boson masses involved in the calculation of the damping rate of the fermion. More recently, in ref. , a fermion-scalar plasma has been considered in the context of real-time formulation and it was found that the Yukawa fermions acquire a width through the induced decay of the scalar in the medium. Here, as in , we use imaginary time formalism and calculate the fermion dispersion relation for interacting massless bosons and fermions in the limit $`k_0,kT`$ and compute the fermion damping rate at rest, $`k=0`$, for the case where the effective dressed (induced by the thermal medium) fermion and boson masses are consistently considered.

In a recent paper  we proposed a modified self-consistent resummation (MSCR) which resums higher-order terms in a non-perturbative way in order to cure the problem of breakdown of the perturbative expansion at finite temperature. We have applied the method up to one  and two-loop  order in the perturbative expansion. We have shown that the MSCR, when applied to the study of the chiral fermion meson model, has the essential features which lead to the satisfaction of Goldstone’s theorem and renormalization of the UV divergences, in the low and high temperature regions. We have explicitly shown that the scheme breaks down around $`T_c`$ i.e., in the region of intermediate temperatures, since quantum fluctuations are known to play a major role there. In this region higher-order terms in the perturbative expansion are required.

It is well known that at high temperature the perturbative expansion can also be broken in theories with spontaneous symmetry breaking (SSB) or in massless field theories (like QCD) because powers of the temperature can compensate for powers of the coupling constant, even if the strength of the coupling is small . Infrared divergence appears. When a set of infrared-divergent diagrams is summed up one gets an infrared-finite result. This is implemented in the MSCR by the recalculation of the self-energy. A comparison between the MSCR and the HTL resummation methods showing its similarities and differences, which is out of the scope of this paper, is under preparation and will be reported elsewhere . So, a motivation to study the fermion dispersion relation and damping rate at high temperature is the fact that the MSCR showed itself an efficient method to execute the resummation in a divergence-free way in this region. Thus, in the computation of the damping rate of the fermion, we shall use at the starting point the dressed by interaction boson mass obtained in  rather than the zero mass parameter of the Lagrangian. In this way, this calculation is interesting since it provides another simple example for the application of the MSCR method.

This paper is organized as follows. In Section II we study the fermion self-energy and obtain the dispersion relation of the fermion with some of its interesting limits. In Section III we compute the fermion damping rate at rest. Section IV is devoted to conclusions.

## II The Fermion Self-Energy

We describe the fermion-bosons vertex by the interaction Lagrangian extracted from the linear sigma model

$$^{int.}=g\overline{\psi }\left[\sigma +i\gamma ^5\stackrel{}{\pi }\stackrel{}{\tau }\right]\psi ,$$
(1)
where $`\psi `$, $`\sigma `$, and $`\pi `$ represent the quark, sigma and pion fields, respectively, and $`g`$ is a non-dimensional positive coupling constant.

The fermion self-energy is defined by

$$𝒟(\omega _n,𝐤)^1=𝒟_0(\omega _n,𝐤)^1+\mathrm{\Sigma }(\omega _n,𝐤),$$
(2)
where $`𝒟_0(\omega _n,𝐩)`$ is the tree-level fermion propagator, expressed as

$$𝒟_0(\omega _n,𝐤)^1=\gamma _\mu k^\mu m_\psi ,$$
(3)
and $`\mathrm{\Sigma }=\mathrm{\Sigma }_s+\gamma ^\mu \mathrm{\Sigma }_\mu `$, with $`\mathrm{\Sigma }_s`$, $`\mathrm{\Sigma }_0`$ and $`\stackrel{}{\mathrm{\Sigma }}`$ being the contributions proportional to the unit, $`\gamma ^0`$ and $`\stackrel{}{\gamma }`$ matrices respectively.

To one-loop order the fermion self-energy expression is given by

$`\mathrm{\Sigma }(k_0,\stackrel{}{k})=\left({\displaystyle \frac{\delta \mathrm{ln}Z_I^{2loop}}{\delta 𝒟_{0\psi }}}\right)_{1PI}=`$ (4)
$`g^2T{\displaystyle \underset{n}{}}{\displaystyle \frac{d^3p}{(2\pi )^3}𝐃_{0\sigma }(\omega _{n+l},𝐩+𝐤)𝒟_{0\psi }(\omega _n,𝐩)}+`$ (5)
$`3g^2T{\displaystyle \underset{n}{}}{\displaystyle \frac{d^3p}{(2\pi )^3}𝐃_{0\pi }(\omega _{n+l},𝐩+𝐤)𝒟_{0\psi }(\omega _n,𝐩)},`$ (6)
since the logarithm of the two-loop interaction partition function is found to be :

$$\mathrm{ln}Z_I^{2loop}=\frac{1}{2}g^2_0^\beta 𝑑\tau _1𝑑\tau _2d^3x_1d^3x_2\frac{[d\varphi ]e^{S_0}[(\overline{\psi }\sigma \psi )^2+(\overline{\psi }i\gamma ^5\stackrel{}{\pi }\stackrel{}{\tau }\psi )^2]}{[d\varphi ]e^{S_0}}.$$
(7)
In eq.(4) $`\omega _n`$ are the Matsubara frequencies, defined as $`\omega _n=2n\pi T`$ for bosons and $`\omega _n=(2n+1)\pi T`$ for fermions and $`𝐃_{0\sigma ,\pi }(\omega _n,𝐩)`$ is the tree-level boson propagator, expressed as

$$𝐃_{0\sigma ,\pi }(\omega _n,𝐤)^1=\omega _n^2+𝐤^2+m_{\sigma ,\pi }^2.$$
(8)
An evaluation of eq.(4) at zero three momentum gives $`\mathrm{\Sigma }(k_0,|𝐤|=0)=(\gamma ^0\mathrm{\Sigma }_0+\mathrm{\Sigma }_s)_\sigma +3(\gamma ^0\mathrm{\Sigma }_0(m_\sigma m_\pi )+\mathrm{\Sigma }_s(m_\sigma m_\pi ))_\sigma `$ with the “$`\gamma ^0`$” and “scalar”- parts of the sigma contribution given respectively by

$`\mathrm{\Sigma }_{0\sigma }=k_0{\displaystyle \frac{g^2}{2}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dpp^2}{\pi ^2}}{\displaystyle \frac{n_\sigma }{\omega _\sigma }}{\displaystyle \frac{k_0^2+\omega _\sigma ^2+\omega _\psi ^2}{[k_0^2(\omega _\psi \omega _\sigma )^2][k_0^2(\omega _\psi +\omega _\sigma )^2]}}+`$ (9)
$`k_0{\displaystyle \frac{g^2}{2}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dpp^2}{\pi ^2}}{\displaystyle \frac{n_\psi }{\omega _\psi }}{\displaystyle \frac{2\omega _\psi ^2}{[k_0^2(\omega _\psi \omega _\sigma )^2][k_0^2(\omega _\psi +\omega _\sigma )^2]}},`$ (10)
$`\mathrm{\Sigma }_{s\sigma }={\displaystyle \frac{g^2}{2}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dpp^2}{\pi ^2}}{\displaystyle \frac{n_\sigma }{\omega _\sigma }}{\displaystyle \frac{m_\psi (k_0^2\omega _\sigma ^2+\omega _\psi ^2)}{[k_0^2(\omega _\psi \omega _\sigma )^2][k_0^2(\omega _\psi +\omega _\sigma )^2]}}+`$ (11)
$`{\displaystyle \frac{g^2}{2}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dpp^2}{\pi ^2}}{\displaystyle \frac{n_\psi }{\omega _\psi }}{\displaystyle \frac{m_\psi (k_0^2\omega _\sigma ^2+\omega _\psi ^2)}{[k_0^2(\omega _\psi \omega _\sigma )^2][k_0^2(\omega _\psi +\omega _\sigma )^2]}},`$ (12)
where $`\omega _{\sigma ,\pi }^2𝐩^2+m_{\sigma ,\pi }^2`$ and $`\omega _\psi ^2𝐩^2+m_\psi ^2`$. The one-loop fermion self-energy is shown in Fig.1.

### A The Dispersion Relation Of Massless Fermions Interacting With Massless Bosons

As a first approximation, in this subsection by considering the interaction of massless fermions with massless bosons we get an effective thermal fermion mass and also calculate the dispersion relation of the fermions.

The poles of the massless fermion propagator ($`m_\psi =\mathrm{\Sigma }_s=0`$) gives the dispersion relation which occurs at the positive-energy root of

$$[k_0(k)+\mathrm{\Sigma }_0(k_0,k)]^2=\left|\stackrel{}{k}+\stackrel{}{\mathrm{\Sigma }}(k_0,\stackrel{}{k})\right|^2,$$
(13)
where $`k|\stackrel{}{k}|`$.

For our intentions, it is sufficient to evaluate eq.(4) in the limit $`k_0,kT`$, and consider the interaction of massless bosons and fermions in order to obtain an effective fermion thermal mass and dispersion relation. Thus,

$$\mathrm{\Sigma }_0(k_0,k)=\frac{1}{8}g^2\frac{T^2}{k}\mathrm{ln}\left|\frac{k_0+k}{k_0k}\right|,$$
(14)
$$\stackrel{}{\mathrm{\Sigma }}(k_0,k)=\frac{1}{4}g^2\frac{T^2}{k^2}\left[\frac{k_0}{2k}\mathrm{ln}\left|\frac{k_0+k}{k_0k}\right|1\right]\stackrel{}{k}\stackrel{~}{\mathrm{\Sigma }}\stackrel{}{k},$$
(15)
where we have defined $`\frac{1}{4}g^2\frac{T^2}{k^2}\left[\frac{k_0}{2k}\mathrm{ln}\left|\frac{k_0+k}{k_0k}\right|1\right]\stackrel{~}{\mathrm{\Sigma }}`$.

Some limits of expressions (14) and (15) are :

$`\mathrm{\Sigma }_0(k_0=0,k)=0,`$ (16)
$`\mathrm{\Sigma }_0(k_0,k=0)={\displaystyle \frac{g^2T^2}{4k_0}}.`$ (17)
$`\stackrel{}{\mathrm{\Sigma }}(k_0,k=0)=0,`$ (18)
$`\stackrel{}{\mathrm{\Sigma }}(k_0=0,k)={\displaystyle \frac{g^2T^2}{4k}}\widehat{k}.`$ (19)
Defining the fermion mass as the location of the pole in the limit $`k=0`$, we have $`k_0+\mathrm{\Sigma }_0(k_0,k0)=0`$, which implies

$$M_\psi ^2=\frac{g^2T^2}{4}.$$
(20)
From (13), we see that the fermion dispersion relation is given by $`k_0+\mathrm{\Sigma }_0=k(1+\stackrel{~}{\mathrm{\Sigma }})`$, that is

$$k_0\frac{1}{2}\frac{M_\psi ^2}{k}\mathrm{ln}\left|\frac{k_0+k}{k_0k}\right|=k\frac{M_\psi ^2}{k}\left[\frac{k_0}{2k}\mathrm{ln}\left|\frac{k_0+k}{k_0k}\right|1\right],$$
(21)
which has the following well known form in the low momentum expansion

$$k_0(k)k_0(k)_+=M_\psi +\frac{1}{3}k+\frac{k^2}{3M_\psi }.$$
(22)
This dispersion relation, $`k_0(k)_+M_\psi +\frac{1}{3}k`$ for $`kM_\psi `$, represents an ordinary fermion whose chirality is equal to its helicity. There is another dispersion relation, $`k_0(k)_{}M_\psi \frac{1}{3}k`$, termed a plasmino , which describes a quasiparticle with chirality opposite to its helicity .

## III The Fermion Damping Rate

Let us now proceed with the computation of the damping rate at rest ($`k=0`$). These calculations will be done considering in the first step $`m_\psi =0`$ and the dressed (effective) boson mass in the internal lines of the fermion self-energy. The effective fermion mass will be taken into account in the second step of the recalculation of the self-energy as the MSCR dictates . In the high temperature region, the bosons dressed masses (given by the MSCR) to be used in internal lines of the fermion self-energy read

$$m_\pi ^2=m_\sigma ^2M_B^2=\frac{g^2}{3}T^2.$$
(23)
So, eq.(9) may be written as

$`\mathrm{\Sigma }_0(k_0,\left|𝐤\right|=0)={\displaystyle \frac{2g^2\gamma ^0k_0}{\pi ^2}}{\displaystyle _0^{\mathrm{}}}𝑑pp^2{\displaystyle \frac{n_B}{\omega _B}}{\displaystyle \frac{k_0^2M_B^22p^2}{[k_0^2M_B^2]^24k_0^2p^2}}+`$ (24)
$`{\displaystyle \frac{4g^2\gamma ^0k_0}{\pi ^2}}{\displaystyle _0^{\mathrm{}}}𝑑pp^4{\displaystyle \frac{n_\psi }{\omega _\psi }}{\displaystyle \frac{1}{[k_0^2M_B^2]^24k_0^2p^2}},`$ (25)
where $`n_B`$ and $`n_\psi `$ are the usual distribution functions for bosons and fermions given respectively by

$$n_B(\omega _B;T)=\frac{1}{e^{\beta \omega _B}1},$$
(26)
$$n_\psi (\omega _\psi ;T)=\frac{1}{e^{\beta \omega _\psi }+1}.$$
(27)
with $`\omega _B\sqrt{𝐩^2+M_B^2}`$ and $`\omega _\psi \left|𝐩\right|`$.

It is worth to note that for $`k_0M_B`$ in eq.(24), it is easy to see that $`\mathrm{\Sigma }_0(k_0,\left|𝐤\right|=0)`$ reduces to the second line of eq.(16) which is a stable state without singularities and we have that there is no decay.

An explicit evaluation of eq.(24) in the pole of the correctec propagator (where the square of the $`\gamma `$ matrices is equal the unit matrix) furnishes

$`\mathrm{\Sigma }_0={\displaystyle \frac{g^2}{4\pi ^2}}\alpha {\displaystyle _0^{\mathrm{}}}𝑑x{\displaystyle \frac{n_B(\stackrel{~}{\omega }_B)}{\stackrel{~}{\omega }_B}}\left[1+{\displaystyle \frac{\beta \alpha }{4}}\left({\displaystyle \frac{1}{x\beta \frac{\alpha }{2}}}{\displaystyle \frac{1}{x+\beta \frac{\alpha }{2}}}\right)\right]`$ (28)
$`{\displaystyle \frac{3g^2}{8\pi ^2}}\alpha {\displaystyle _0^{\mathrm{}}}𝑑xn_\psi (\stackrel{~}{\omega }_\psi )\left({\displaystyle \frac{1}{x\beta \frac{\alpha }{2}}}+{\displaystyle \frac{1}{x+\beta \frac{\alpha }{2}}}\right){\displaystyle \frac{3g^2T^2}{\pi ^2\alpha }}{\displaystyle _0^{\mathrm{}}}𝑑xx^2\left({\displaystyle \frac{n_B(\stackrel{~}{\omega }_B)}{\stackrel{~}{\omega }_B}}+{\displaystyle \frac{n_\psi (\stackrel{~}{\omega }_\psi )}{\stackrel{~}{\omega }_\psi }}\right),`$ (29)
with the definitions $`\alpha k_0\frac{M_B^2}{k_0}`$, $`\stackrel{~}{\omega }_B\sqrt{\beta ^2𝐩^2+\beta ^2M_B^2}`$, $`\stackrel{~}{\omega }_\psi \left|\beta 𝐩\right|`$ and $`\beta px`$. The interesting physics happens when $`k_0>M_B`$. Otherwise (if $`k_0<M_B`$) one would get imaginary (forbidden) frequency.

The expression for $`\mathrm{\Sigma }_0`$ in (28) has singularities, and now we adopt the prescription $`\alpha =\omega i\gamma `$, since in general $`k_0`$ is complex, where $`\omega `$ is the real frequency and $`\gamma `$ is a real constant. With this assumption for $`\alpha `$, eq.(28) is expressed as

$`\mathrm{\Sigma }_0={\displaystyle \frac{g^2}{4\pi ^2}}\alpha {\displaystyle _0^{\mathrm{}}}𝑑x{\displaystyle \frac{n_B(\stackrel{~}{\omega }_B)}{\stackrel{~}{\omega }_B}}\left[1+{\displaystyle \frac{\alpha }{2}}\left({\displaystyle \frac{\frac{2x}{\beta }\omega }{\left(\frac{2x}{\beta }\omega \right)^2+\gamma ^2}}{\displaystyle \frac{\frac{2x}{\beta }+\omega }{\left(\frac{2x}{\beta }+\omega \right)^2+\gamma ^2}}\right)\right]+`$ (30)
$`i{\displaystyle \frac{g^2}{8\pi ^2}}\alpha ^2{\displaystyle _0^{\mathrm{}}}𝑑x{\displaystyle \frac{n_B(\stackrel{~}{\omega }_B)}{\stackrel{~}{\omega }_B}}\left[{\displaystyle \frac{\gamma }{\left(\frac{2x}{\beta }\omega \right)^2+\gamma ^2}}+{\displaystyle \frac{\gamma }{\left(\frac{2x}{\beta }+\omega \right)^2+\gamma ^2}}\right]+`$ (31)
$`{\displaystyle \frac{3g^2}{4\pi ^2}}{\displaystyle \frac{\alpha }{\beta }}{\displaystyle _0^{\mathrm{}}}𝑑xn_\psi (\stackrel{~}{\omega }_\psi )\left[{\displaystyle \frac{\frac{2x}{\beta }\omega }{\left(\frac{2x}{\beta }\omega \right)^2+\gamma ^2}}+{\displaystyle \frac{\frac{2x}{\beta }+\omega }{\left(\frac{2x}{\beta }+\omega \right)^2+\gamma ^2}}\right]+`$ (32)
$`i{\displaystyle \frac{3g^2}{4\pi ^2}}{\displaystyle \frac{\alpha }{\beta }}{\displaystyle _0^{\mathrm{}}}𝑑xn_\psi (\stackrel{~}{\omega }_\psi )\left[{\displaystyle \frac{\gamma }{\left(\frac{2x}{\beta }\omega \right)^2+\gamma ^2}}+{\displaystyle \frac{\gamma }{\left(\frac{2x}{\beta }+\omega \right)^2+\gamma ^2}}\right]+`$ (33)
$`{\displaystyle \frac{3g^2T^2}{\pi ^2\alpha }}{\displaystyle _0^{\mathrm{}}}𝑑xx^2{\displaystyle \frac{n_B(\stackrel{~}{\omega }_B)}{\stackrel{~}{\omega }_B}}{\displaystyle \frac{3g^2T^2}{\pi ^2\alpha }}{\displaystyle _0^{\mathrm{}}}𝑑xx^2{\displaystyle \frac{n_\psi (\stackrel{~}{\omega }_\psi )}{\stackrel{~}{\omega }_\psi }}.`$ (34)
Now making use of the definition of the delta function

$$\delta (y)=\stackrel{lim}{ϵ0}\frac{1}{\pi }\frac{ϵ}{y^2+ϵ^2},$$
(35)
and the definitions

$$F(x,\omega )\frac{\frac{2x}{\beta }\omega }{\left(\frac{2x}{\beta }\omega \right)^2+\gamma ^2}\frac{\frac{2x}{\beta }+\omega }{\left(\frac{2x}{\beta }+\omega \right)^2+\gamma ^2},$$
(36)
and

$$G(x,\omega )\frac{\frac{2x}{\beta }\omega }{\left(\frac{2x}{\beta }\omega \right)^2+\gamma ^2}+\frac{\frac{2x}{\beta }+\omega }{\left(\frac{2x}{\beta }+\omega \right)^2+\gamma ^2},$$
(37)
we get

$`\mathrm{\Sigma }_0={\displaystyle \frac{g^2}{4\pi ^2}}\alpha {\displaystyle _0^{\mathrm{}}}𝑑x{\displaystyle \frac{n_B(\stackrel{~}{\omega }_B)}{\stackrel{~}{\omega }_B}}\left[1+{\displaystyle \frac{\alpha }{2}}F(x,\omega )\right]+i{\displaystyle \frac{g^2}{8\pi }}\left[{\displaystyle \frac{\omega ^2}{\sqrt{\frac{\omega ^2}{4}+M_B^2}}}{\displaystyle \frac{1}{e^{\beta \sqrt{\frac{\omega ^2}{4}+M_B^2}}1}}\right]+`$ (38)
$`{\displaystyle \frac{3g^2}{4\pi ^2}}{\displaystyle \frac{\alpha }{\beta }}{\displaystyle _0^{\mathrm{}}}𝑑xn_\psi (\stackrel{~}{\omega }_\psi )G(x,\omega )+i{\displaystyle \frac{3g^2}{4\pi }}{\displaystyle \frac{|\omega |}{e^{\frac{\beta |\omega |}{2}}+1}}+`$ (39)
$`{\displaystyle \frac{3g^2T^2}{\pi ^2\alpha }}{\displaystyle _0^{\mathrm{}}}𝑑xx^2{\displaystyle \frac{n_B(\stackrel{~}{\omega }_B)}{\stackrel{~}{\omega }_B}}{\displaystyle \frac{3g^2T^2}{\pi ^2\alpha }}{\displaystyle _0^{\mathrm{}}}𝑑xx^2{\displaystyle \frac{n_\psi (\stackrel{~}{\omega }_\psi )}{\stackrel{~}{\omega }_\psi }}.`$ (40)
Here we use some results from high temperature expansion of one-loop integrals derived by Dolan and Jackiw in :

$$_0^{\mathrm{}}𝑑x\frac{n_B(\stackrel{~}{\omega }_B)}{\stackrel{~}{\omega }_B}=\frac{\pi T}{2M_B}+\frac{1}{2}\mathrm{ln}\left(\frac{M_B}{4\pi T}\right)+O\left(\frac{M_B^2}{T^2}\right)\frac{\pi \sqrt{3}}{2g},$$
(41)
where in eq.(41) we have used from (23) that $`M_B=\frac{gT}{\sqrt{3}}`$. This means that the first term in the r.h.s. of eq.(38) is $`\frac{g\sqrt{3}}{8\pi }1`$. On the other hand, for the last two terms in the r.h.s. of this same expression, we have

$$\frac{T^2}{\pi ^2}_0^{\mathrm{}}𝑑xx^2\frac{n_B(\stackrel{~}{\omega }_B)}{\stackrel{~}{\omega }_B}+\frac{T^2}{\pi ^2}_0^{\mathrm{}}𝑑xx^2\frac{n_\psi (\stackrel{~}{\omega }_\psi )}{\stackrel{~}{\omega }_\psi }=\left(\frac{T^2}{6}\frac{M_BT}{2\pi }+O(g^2)\right)+\left(\frac{T^2}{12}\right)\frac{T^2}{4}.$$
(42)
The parts involving $`F(x,\omega )`$ and $`G(x,\omega )`$ are less important contributions in comparison to the dominant term that is proportional to $`\frac{g^2T^2}{\omega }`$, mainly for small $`\omega `$. So, in a first glance one can neglect them. Putting these results in eq.(38) and assuming ($`\gamma \omega `$), the leading term of the real frequency can be written as

$$\omega _1^2=\frac{9g^2T^2}{4}M_{\psi ,1}^2.$$
(43)
The next step is the recalculation of the self-energy of the fermion to get the second order corrected fermion mass from the pole location

$$k_{0,2}M_{\psi ,1}+\mathrm{\Sigma }(M_B,M_{\psi ,1},k=0)=0,$$
(44)
where from this equation we will identify $`(Rek_{0,2})^2=\omega _2^2M_{\psi ,2}^2`$ and $`Imk_{0,2}=\gamma `$. With this resummation procedure it is possible to take into account the induced by the thermal medium fermion mass. At each recalculation we take into account more infinity subsets of graphs which result is shown in the mass. Now we are able to find out the damping rate of the fermion, which is defined by the imaginary part of the self-energy on-shell .

From eq.(4) the self-energy $`\mathrm{\Sigma }(M_B,M_{\psi ,1},k=0)=\mathrm{\Sigma }_s(M_B,M_{\psi ,1},k=0)+\mathrm{\Sigma }_0(M_B,M_{\psi ,1},k=0)`$ in eq.(44) to be evaluated now is

$`\mathrm{\Sigma }=2k_0g^2{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dpp^2}{\pi ^2}}{\displaystyle \frac{n_B}{\omega _B}}{\displaystyle \frac{k_0^2(M_B^2+M_{\psi ,1}^2+2p^2)}{[k_0^2(M_B^2+M_{\psi ,1}^2)]^24k_0^2p^24M_B^2M_{\psi ,1}^2}}+`$ (45)
$`4k_0g^2{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dpp^2}{\pi ^2}}{\displaystyle \frac{n_{\psi ,1}}{\omega _{\psi ,1}}}{\displaystyle \frac{M_{\psi ,1}^2+p^2}{[k_0^2(M_B^2+M_{\psi ,1}^2)]^24k_0^2p^24M_B^2M_{\psi ,1}^2}}+`$ (46)
$`2g^2{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dpp^2}{\pi ^2}}{\displaystyle \frac{n_B}{\omega _B}}{\displaystyle \frac{M_{\psi ,1}(k_0^2M_B^2+M_{\psi ,1}^2)}{[k_0^2(M_B^2+M_{\psi ,1}^2)]^24k_0^2p^24M_B^2M_{\psi ,1}^2}}+`$ (47)
$`2g^2{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dpp^2}{\pi ^2}}{\displaystyle \frac{n_{\psi ,1}}{\omega _{\psi ,1}}}{\displaystyle \frac{M_{\psi ,1}(k_0^2M_B^2+M_{\psi ,1}^2)}{[k_0^2(M_B^2+M_{\psi ,1}^2)]^24k_0^2p^24M_B^2M_{\psi ,1}^2}}`$ (48)
where the last term ($`4M_B^2M_{\psi ,1}^2`$) in the denominators in the r.h.s. of eq.(45) will be neglected, since it is of $`O(g^4)`$.

Doing again the calculations necessary to get the real and imaginary parts of the poles, where now $`\alpha k_0\frac{M_B^2+M_{\psi ,1}^2}{k_0}`$, $`\stackrel{~}{\omega }_B\sqrt{\beta ^2𝐩^2+\beta ^2M_B^2}`$ and $`\stackrel{~}{\omega }_{\psi ,1}\sqrt{\beta ^2𝐩^2+\beta ^2M_{\psi ,1}^2}`$, one finds

$$\omega _2^2=M_{\psi ,1}^2\left(1\frac{2g}{\pi \sqrt{3}}\right)M_{\psi ,2}^2,$$
(49)
$$\gamma =\frac{3g^2}{8\pi }\left[\frac{\omega _2^2}{\sqrt{\frac{\omega _2^2}{4}+M_B^2}}\frac{1}{e^{\beta \sqrt{\frac{\omega _2^2}{4}+M_B^2}}1}\right]+\frac{9g^2}{8\pi }\left[\frac{\omega _2^2}{\sqrt{\frac{\omega _2^2}{4}+M_{\psi ,1}^2}}\frac{1}{e^{\beta \sqrt{\frac{\omega _2^2}{4}+M_{\psi ,1}^2}}+1}\right].$$
(50)
It is important to note that eq.(49) allows us to define a critical value of $`g`$ , that is $`g_{cr}2.7`$. This value agrees with the one obtained in ref . As pointed out by M. Le Bellac in ref. , collective excitations should disappear at the critical temperature, althoug we do not expect pertubation theory is still valid for such large values of $`g`$. Equations (49) and (50) has the following interpretation: The leading order of the real frequency $`\omega _2`$ is of order $`gT`$, in concordance with (20). The damping rate is proportional to the probability $`n_\psi (\frac{1}{2}\omega )`$ of having a fermion with energy $`\frac{1}{2}\omega `$ and a probability $`n_B(\frac{1}{2}\omega )`$ of having a boson with energy $`\frac{1}{2}\omega `$. These probabilities are weighted by numerical factors and the available phase space $`\omega `$ . The distribution functions can be expanded for low energies, $`n_B(\frac{1}{2}\omega )T/\omega `$ and $`n_\psi (\frac{1}{2}\omega )1/2`$ and the damping rate (up to order $`g^3`$) reduces to

$$\gamma \frac{9g^2T}{16\pi }+\frac{27g^3T}{32\pi }.$$
(51)
## IV Concluding Remarks

In this paper, we have considered the fermion boson interaction at finite temperature. First, we have calculated the self-energy of the fermions due the interaction with scalar-bosons and pseudoscalar-bosons in the framework of the linear sigma model.

Next, we have calculated the fermion dispersion relation in the limit $`k_0,kT`$ of massless fermions interacting with massless bosons and some of its limits. Also, we have obtained the thermal fermion mass which is of order $`gT`$.

Finally, we have computed the frequency and the damping rate of the fermion at rest, considering the dressed fermion and boson masses in the internal lines of the fermion self-energy rather than the zero mass parameter of the Lagrangian. The damping rate of the fermion was found to be of order $`g^2T`$ from the boson internal line of the fermion self-energy plus a part which is of order $`g^3T`$ from the fermion internal line of the self-energy. This is a remarkable result since it shows the signature of the alternative MSCR method (up to second order of the non-perturbative correction) in the desired effective mass and damping rate. One important feature of this calculation (up to this order) is that we have gotten this result algebraically in a clear manner differently from the cases where the results are reached numerically.

The calculation of the fermion damping rate at rest constitutes another simple but instructive application of the MSCR method.

## Acknowledgements

The author thanks the hospitality given by the Nuclear Theory group during his visit at the University of Minnesota were this work was initiated. He is gratefully indebted to Professors Joe Kapusta and Paul Ellis for various helpful discussions. Also, the author would like to thank Professor Paul Ellis for a critical reading of the manuscript.


format: text

text:
# A new class of semiclassical wave function uniformizations

## I Introduction

Semiclassical methods are based on classical mechanics; the relevant classical manifolds form the “skeleton” to which the wave function is attached. However, associating a semiclassical wave function with a classical manifold is problematic when enclosed loops in phase space have area less than $`\mathrm{}`$. If that happens, typically two or more stationary-phase points corresponding to distinct contributions to the semiclassical wave function are not well separated in phase. Consequently, the stationary-phase (SP) method, on which semiclassical approximation is based, breaks down. Worse, in some important situations, such as the universal homoclinic oscillations associated with chaotic regions of phase space, the bad enclosed loops may be repeated many times in a small region.

Although the semiclassical wave function diverges near caustics lurking near the small area loops, the exact quantum wave function is of course well-behaved. This is one manifestation of quantum smoothing over classical detail. We shall quantify below how this smoothing occurs. We suggest and test a new way of looking at the smoothing process in terms of “replacement manifolds”, in which new, well behaved classical manifolds are substituted for the original badly behaved ones. One application of this approach is uniformizing a semiclassical wave function in the vicinity of the ubiquitous homoclinic oscillations of a chaotic system.

It is widely understood that very small changes in potentials or walls of billiards can have little effect on quantum eigenstates, but can drastically affect the classical manifold structure. This immediately implies a many-to-one relationship between the classical manifolds and quantum wave functions. That is, many different underlying classical manifold patterns correspond to the same wave function. Now, whatever the choice of manifold (amongst these equivalent forms) we can suppose an appropriate uniformization exists which gives the correct wave function. However, the non-uniformized, simple semiclassical limit for the various choices of manifolds can differ drastically in accuracy. It is this fact which we exploit in the present paper.

## II Standard uniformization: The Airy function

The paradigm example of uniformization is the Airy function, i.e. the energy eigenfunction of the linear ramp potential. We take the potential to be of the form

$$V(x)=\beta x$$
(1)
The Schrödinger equation is

$$\frac{\mathrm{}^2}{2m}\psi _E^{^{\prime \prime }}(x)(\beta x+E)\psi _E(x)=0$$
(2)
After the shift

$$\psi _E(x)=\psi _0(x+E/\beta )$$
(3)
and scaling transformation

$$u=\frac{m^{1/3}\beta ^{1/3}}{\mathrm{}^{2/3}}x\gamma x$$
(4)
the Schrödinger equation reads

$$\frac{1}{2}\mathrm{\Psi }^{^{\prime \prime }}(u)+u\mathrm{\Psi }(u)=0;$$
(5)
with $`\mathrm{\Psi }(u)=\psi _0(x)`$. Eq. (5) is exactly solved by going to the momentum representation. We have

$$\left(\frac{1}{2}p^2i\frac{}{p}\right)\mathrm{\Phi }(p)=0$$
(6)
with solution

$$\mathrm{\Phi }(p)=e^{ip^3/6}$$
(7)
Thus

$`\mathrm{\Psi }(u)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi }}}{\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}e^{ipu}\mathrm{\Phi }(p)𝑑p`$ (8)
$`=`$ $`{\displaystyle \frac{2}{\sqrt{2\pi }}}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}\mathrm{cos}(p^3/6pu)𝑑p`$ (9)
$`=`$ $`2^{1/2}\mathrm{Ai}\left(2^{1/3}u\right)`$ (10)
and

$`\psi (x)=\eta \mathrm{Ai}(2^{1/3}\gamma x)=\eta \mathrm{Ai}\left({\displaystyle \frac{2^{1/3}m^{1/3}\beta ^{1/3}}{\mathrm{}^{2/3}}}x\right)`$

SP approximation for the momentum integral gives

$`\mathrm{\Psi }(u)`$ $`=`$ $`\left[{\displaystyle \frac{2\pi }{i(\sqrt{2u})}}\right]^{1/2}e^{(\sqrt{2u})^3/6+i\sqrt{2u}u}+\mathrm{c}.\mathrm{c}.`$ (11)
$``$ $`{\displaystyle \frac{1}{u^{1/4}}}\mathrm{cos}\left({\displaystyle \frac{2}{3}}\sqrt{2}u^{3/2}{\displaystyle \frac{\pi }{4}}\right)`$ (12)
Noting that

$$^up(u^{})𝑑u^{}=\frac{(2u)^{3/2}}{3}$$
(13)
we get finally

$`\mathrm{\Psi }(u){\displaystyle \frac{1}{\sqrt{p(u)}}}\mathrm{cos}\left({\displaystyle ^u}p(u^{})𝑑u^{}{\displaystyle \frac{\pi }{4}}\right)`$

This is the standard WKB form for the linear ramp potential.

Figure 1 illustrates the situation for the semiclassical position wave function, which becomes inaccurate when the area enclosed by the vertical $`x`$-line and the energy contour corresponding to the state $`|E`$ falls below Planck’s constant. Upon entering the “bad” region, the original momentum integral may be substituted for the WKB result, giving a uniformization. We have an integral over “initial momentum,” making this example a type of “initial value representation” (IVR) as first discussed by Miller. Of course in this case we know the result of the integral is an Airy function, but in general there may be integrals which are not known analytically or even close to ones that are known.

In this paper we encounter more complex situations which are induced by nonlinear interactions, including chaos, reflection from corrugated surfaces, etc. We will show that intricate manifolds can sometimes be replaced by much simpler ones, which may themselves be evaluated semiclassically in any basis, using a much simplified uniformization if required.

## III Replacement manifolds

We now introduce an example of such a problematic Lagrangian manifold, defined by a position-momentum relation

$$p(x)=p_0\frac{3}{2}\mathrm{}\lambda \alpha x^2\mathrm{sin}\alpha x^3.$$
(14)
This model was first introduced in Ref. but was not explored there. The function (14) has oscillations of increasing amplitude and frequency, which however have the same area $`p(x)𝑑x=\lambda \mathrm{}`$ between successive zeros (see Fig. 2, small value of $`\lambda `$ corresponds to small area of loops). This is also a characteristic property of homoclinic oscillations near a periodic orbit. The momentum $`p(x)`$ spans an ever larger range as $`x`$ increases, but we now show that for $`\lambda <1`$ the manifold’s semiclassical behavior can be understood by replacing it with three smooth manifolds. The rules for replacement are simple to derive.

The action function $`S(x)`$ is the integral of $`p(x)`$:

$$S(x)=_0^xp(x^{})𝑑x^{}=p_0x+\frac{1}{2}\mathrm{}\lambda \left(\mathrm{cos}\alpha x^31\right).$$
(15)
The “wave function”

$$\psi (x)=A(x)e^{iS(x)/\mathrm{}}$$
(16)
can be approximated by

$$\psi (x)A(x)e^{i\lambda /2}e^{ip_0x/\mathrm{}}\left(1+i\frac{\lambda }{2}\mathrm{cos}(\alpha x^3)\right)$$
(17)
$`\psi (x)`$ $``$ $`A(x)e^{i\lambda /2}e^{ip_0x/\mathrm{}}`$ (19)
$`\times \left(1+i{\displaystyle \frac{\lambda }{4}}e^{i\alpha x^3}+i{\displaystyle \frac{\lambda }{4}}e^{i\alpha x^3}\right)`$
if $`\lambda <1`$. For small $`\lambda `$ this expression is nearly identical with (16) with the original action (15), because $`\mathrm{cos}(\alpha x^3)`$ is bounded by $`\pm 1`$. However semiclassically it has the interpretation of the sum of three smooth classical manifolds, i.e. $`p=p_0`$, $`p=p_0+3\mathrm{}\alpha x^2`$ and $`p=p_03\mathrm{}\alpha x^2`$. They have the weights $`A(x)`$, $`i\frac{\lambda }{4}e^{i\lambda /2}A(x)`$, $`i\frac{\lambda }{4}e^{i\lambda /2}A(x)`$, respectively. The situation is depicted in Fig. 3.

Almost every discussion of the relation of classical and quantum mechanics for chaotic systems alludes to quantum smoothing, but here we have seen explicitly one way this smoothing comes about. For all reasonable purposes the three smooth classical manifolds accurately replace the rapid oscillations of the original manifold. Note too that depending on the parameters the outlying manifolds $`p=p_0+3\mathrm{}\alpha x^2,`$ and $`p=p_03\mathrm{}\alpha x^2`$ lie far beyond the limits of the original distributions. We shall discuss interesting consequences of this in the next section.

So far, we have considered the wave function only in position representation where caustics are absent and the semiclassical form (16) was accurate. The situation deteriorates drastically in momentum space. In order to find the wave function in momentum representation, we have to sum over all contributions from intersections of a horizontal line (corresponding to a momentum eigenstate) with the oscillating manifold (14). When $`\lambda <1`$, the adjacent intersections will be separated by a phase $`\mathrm{\Delta }S/\mathrm{}`$ that is smaller than 1 for any classically allowed momentum, and therefore the standard semiclassical approximation will break down for all classical momenta! Usual Airy-type uniformization methods of Section II (also see ) do not serve under these circumstances, because the manifold is dense with caustics (see the filled-in areas $`\mathrm{\Delta }S_1`$ and $`\mathrm{\Delta }S_2`$ in Fig. 2). Apparently, we have no other choice but to part with semiclassical method and use a numerical Fourier transform (IVR) over the position wave function to obtain $`\varphi (p)`$. However using the replacement manifolds (RMs) is a much simpler and more intuitive approach which can rescue the situation.

To find the momentum wave function $`\varphi (p)`$ associated with the original manifold we add (with appropriate weights) wave functions $`\varphi _n(p)`$ corresponding to the three RMs. Each of these partial wave functions is equal to a Fourier transform of a corresponding position wave function $`\psi _n(x)`$,

$$\varphi _n(p)=\frac{1}{\sqrt{2\pi \mathrm{}}}𝑑x\psi _n(x)e^{ipx/\mathrm{}}.$$
(20)
Except for a small region close to $`p_0`$ (highlighted in Fig. 3), this integral can be evaluated by SP method! In other words, although the semiclassical approximation completely fails for the original manifold, it works almost everywhere for the replacement manifolds. Moreover, as we show in the following section, the standard Airy uniformization procedure of Section II can be exploited to correct the inaccuracy in a narrow region near $`p_0`$. We will see that the momentum wave function found by applying SP approximation to RMs is in excellent agreement with numerical solution. We emphasize that the RMs do not yield the original, badly behaved semiclassical result, but rather something much more accurate and much simpler.

## IV Uniformization with replacement manifolds in a generalized model

In order to capture other properties of homoclinic oscillations in our model, we consider a generalized manifold

$$p(x)=p_0\frac{1}{2}\mathrm{}\lambda \frac{d\xi (x)}{dx}\mathrm{sin}\xi (x),$$
(21)
in which $`\alpha x^3`$ from (14) is replaced by a generic smooth function $`\xi (x)`$. The presence of the derivative $`\frac{d\xi (x)}{dx}`$ ensures equal loop areas between successive zeros of $`\xi (x)`$. Example of a manifold with $`\xi (x)=\alpha \mathrm{log}\frac{x}{a}`$ is displayed in Fig. 5. Manifold (21) can be obtained from a horizontal manifold representing a momentum state $`|p_0`$ by a canonical transformation of coordinates

$$x=x_0,p=p_0\frac{1}{2}\mathrm{}\lambda \frac{d\xi }{dx_0}\mathrm{sin}\xi (x_0),$$
(22)
generated by a function

$$F_2(x,p_0)=p_0x+\frac{1}{2}\mathrm{}\lambda \mathrm{cos}\xi (x),$$
(23)
We find the uniform semiclassical transformation element $`p|p_0`$, and thereby fix the amplitude prefactor to be $`A(x)=(2\pi \mathrm{})^{1/2}`$; this however retains all the problematic features of the manifold. We solve the problem to all orders in $`\lambda `$ by replacing the original manifold with an infinite series of RMs. First we note that the semiclassical transformation element

$$x|p_0_{sc}=\frac{1}{\sqrt{2\pi \mathrm{}}}e^{iF_2(x,p_0)/\mathrm{}}$$
(24)
is fairly accurate since there are no caustics in this representation (to each final $`x`$ corresponds a single initial $`x_0`$). The momentum-representation element is obtained by a Fourier transform,

$`p|p_0`$ $``$ $`{\displaystyle \frac{1}{\sqrt{2\pi \mathrm{}}}}{\displaystyle 𝑑xe^{ipx/\mathrm{}}x|p_0_{sc}}`$ (25)
$`=`$ $`{\displaystyle \frac{1}{2\pi \mathrm{}}}{\displaystyle 𝑑x\mathrm{exp}\left\{\frac{i}{\mathrm{}}\left[px+F_2(x,p_0)\right]\right\}}.`$ (26)
If evaluated by SP method, the integral yields a very poor semiclassical result since there exist many coalescing SP points. An accurate uniform answer is obtained by evaluating integral (25) exactly. The accuracy is usually further improved by changing the integration variable from $`x`$ to $`x_0`$ (in our case both forms are equivalent since $`x=x_0`$). This uniform version of $`p|p_0`$ is in IVR form, but can best be evaluated by writing $`x|p_0`$ as a sum over RMs. Recognizing that the factor $`e^{i\frac{\lambda }{2}\mathrm{cos}\xi (x)}`$ in $`x|p_0_{sc}=(2\pi \mathrm{})^{1/2}\mathrm{exp}\left[ip_0x/\mathrm{}+i\frac{\lambda }{2}\mathrm{cos}\xi (x)\right]`$ is a generating function for Bessel functions, we can extend the sum from previous section beyond the first order in $`\lambda `$. In fact we obtain an infinite sum convergent for any $`\lambda `$,

$`x|p_0`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi \mathrm{}}}}e^{ip_0x/\mathrm{}}{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}J_n\left({\displaystyle \frac{\lambda }{2}}\right)i^ne^{in\xi (x)}`$ (27)
$`=`$ $`{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}J_n\left({\displaystyle \frac{\lambda }{2}}\right)i^nx|p_0_n`$ (28)
where $`x|p_0_n=(2\pi \mathrm{})^{1/2}\mathrm{exp}\left\{i\left[p_0x+n\mathrm{}\xi (x)\right]/\mathrm{}\right\}`$. We can rewrite (25) as

$`p|p_0`$ $`=`$ $`{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}J_n\left({\displaystyle \frac{\lambda }{2}}\right){\displaystyle \frac{i^n}{\sqrt{2\pi \mathrm{}}}}{\displaystyle 𝑑xe^{ipx/\mathrm{}}x|p_0_n}`$ (29)
$`=`$ $`{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}J_n\left({\displaystyle \frac{\lambda }{2}}\right)i^np|p_0_n`$ (30)
where

$$p|p_0_n=\frac{1}{2\pi \mathrm{}}𝑑x\mathrm{exp}\left\{\frac{i}{\mathrm{}}\left[(p_0p)x+n\mathrm{}\xi (x)\right]\right\}$$
(31)
may be identified semiclassically as a transformation element corresponding to the $`n`$-th RM

$$p_n(x)=p_0+n\mathrm{}\frac{d\xi (x)}{dx},$$
(32)
This manifold for most functions $`\xi (x)`$ of interest contains no caustics, and consequently $`p|p_0_n`$ allows evaluation by SP method! As promised, we have expressed the uniform version of $`p|p_0`$ as a weighted sum over semiclassical replacement manifolds,

$`p|p_0_{unif}={\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}J_n\left({\displaystyle \frac{\lambda }{2}}\right)i^np|p_0_{n,sc}.`$

If applied to a well-behaved (i.e. decaying fast enough at $`\pm \mathrm{}`$) function of $`p_0`$, the resulting sum converges for any $`\lambda `$. The replacement-manifold method is therefore not restricted to the regime where loop areas are smaller than $`\mathrm{}`$. It also works in the strongly chaotic regime where $`\lambda >1`$ (and where the standard semiclassical approximation holds) if all manifolds up to $`|n|\lambda `$ are included in the sum. This follows since $`J_n\left(\lambda /2\right)`$ considered as a function of $`n`$ decays exponentially fast for $`|n|>\lambda /2`$. In other words, we only need to include several RMs beyond those intersecting the original manifold. Physically, the manifolds outside of the range of the original one lie in a classically forbidden region in which the wave function is exponentially suppressed.

In opposite case ($`\lambda <1`$) which interests us the most, considering RMs up to $`|n|=1`$ will suffice while the “primitive” semiclassical approximation fails completely even for $`p=p_0`$. In this regime $`\lambda `$ can be thought of as a parameter defining the strength of a perturbation which causes the manifold to oscillate around a simple manifold that describes the system if perturbation is absent.

To complete our solution for a specific function $`\xi (x)`$, we must evaluate $`p|p_0_n`$. For $`n=0`$, we cannot use SP approximation (because the action is linear), but integral (31) is trivial and we obtain in general

$$p|p_0_0=\delta (pp_0).$$
(33)
For $`n0`$, SP method applied to (31) yields a result in the form of a sum over SP points $`x_{sp}`$,

$`p|p_0_{n,sc}={\displaystyle \frac{1}{\sqrt{2\pi \mathrm{}}}}`$ (34)
$`\times {\displaystyle \underset{x=x_{sp}}{}}|{\displaystyle \frac{^2f_n}{x^2}}|^{1/2}\mathrm{exp}\left[i({\displaystyle \frac{f_n}{\mathrm{}}}+{\displaystyle \frac{\pi }{4}}sgn{\displaystyle \frac{^2f_n}{x^2}})\right]`$ (35)
where $`f_n(x,p,p_0)=(p_0p)x+n\mathrm{}\xi (x)`$ and $`\frac{f_n}{x}=0`$ for $`x=x_{sp}`$. For manifold (14), each RM $`p_n(x)=p_0+3n\mathrm{}\alpha x^2`$ has two SP points

$$x_{sp}=\pm \left(\frac{pp_0}{3n\mathrm{}\alpha }\right)^{1/2}$$
(36)
whose contributions add to give

$`p|p_0_{n,sc}`$ $`=`$ $`\pi ^{1/2}[3n\mathrm{}^3\alpha (pp_0)]^{1/4}`$ (38)
$`\times \mathrm{cos}\left[{\displaystyle \frac{2}{(n\alpha )^{1/2}}}\left({\displaystyle \frac{pp_0}{3\mathrm{}}}\right)^{3/2}{\displaystyle \frac{\pi }{4}}\right].`$
RMs still contain one caustic at $`p=p_0`$, which can be easily uniformized by evaluating integral (31) exactly,

$`p|p_0_{n,unif}={\displaystyle \frac{1}{(3|n|\alpha )^{1/3}\mathrm{}}}Ai\left({\displaystyle \frac{(sgnn)(pp_0)}{(3|n|\alpha )^{1/3}\mathrm{}}}\right).`$

Having found $`p|p_0_n`$, we can calculate $`p|p_0`$ from (29) to all orders in $`\lambda `$. $`p|p_0_{1,sc}`$ and $`p|p_0_{1,unif}`$ agree very well except for a small region near $`p_0`$ (similarly as in Fig. 1). Fig. 4 demonstrates the excellent agreement of RM method with the direct numerical computation of $`p|p_0`$ using fast Fourier transform (IVR).

Here the replacement-manifold approach should be compared to the semiclassical perturbation (SCP) approximation of Miller and Smith  who used perturbative classical dynamics to calculate the semiclassical $`S`$ matrix. They first applied this method to the collinear collision of an atom with a diatomic molecule . The starting point is the initial value representation of the classical scattering matrix,

$`S_{n_f,n_i}`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{2\pi }}𝑑q_i\left[{\displaystyle \frac{q_f(q_i,n_i)}{q_i}}\right]^{1/2}`$ (39)
$`\times `$ $`\mathrm{exp}\left(i\left\{\mathrm{\Phi }(q_i,n_i)+q_f(q_i,n_i)\left[n_f(q_i,n_i)n_f\right]\right\}\right)`$ (40)
where $`\mathrm{}=1`$, $`n_i`$, $`n_f`$, $`q_i`$, $`q_f`$ denote, respectively, the initial and final values of the action and angle variables describing the internal degree of freedom, and $`\mathrm{\Phi }`$ is the action integral. For details, see Ref.  and . If the classical quantities listed above are calculated using the first-order perturbation dynamics, the $`S`$ matrix takes the form

$`S_{n_f,n_i}={\displaystyle \frac{e^{i\mathrm{\Phi }_0}}{2\pi }}{\displaystyle _0^{2\pi }}𝑑q_i\mathrm{exp}\left\{i\left[\left(n_fn_i\right)q_i+A(q_i,n_i)\right]\right\}`$

where $`\mathrm{\Phi }_0`$ is twice the phase-shift for the unperturbed potential and $`A(q_i,n_i)`$ is the time integral of the perturbation potential along the unperturbed trajectory . For various physical systems (such as ion-dipole collisions  or atom-surface scattering ), $`A(q_i,n_i)`$ has a sinusoidal dependence on $`q_i`$ as long as perturbation remains small. In this regime we expect the SCP and RM approximations to give the same answer. We shall verify that for the scattering by a corrugated wall in Section VII. The SCP approximation of Miller and Smith presently appears to be more widely applicable since it does not require any special property of the manifold. On the other hand, advantage of the RM method lies in the fact that it is not subjected to the smallness of perturbation. Provided that the full action falls into one of the classes discussed above, we can apply the RM method without approximating classical dynamics.

## V Model of homoclinic oscillations

Homoclinic oscillations in chaotic systems have another characteristic property: the amplitude of oscillations (i.e. “height” of loops) increases exponentially as we approach an unstable periodic orbit (see Fig. 5).

To keep area of loops constant, the width of loops has to decrease accordingly. Because systems in nature are usually bounded, the exponentially growing loops must eventually twist in a complicated way to fit into available phase space. For sake of simplicity, we consider a model that satisfies the requirement of exponential growth, but is unbounded. Such a model can be found among generalized manifolds of previous section. Specifically, we want to find a manifold whose homoclinic points (i.e. intersections with line $`p=p_0`$ in our case) form a geometric series,

$$x_n=ae^{n\lambda _+}$$
(41)
where $`\lambda _+`$ is the Lyapunov exponent. In other words, we are looking for a function $`\xi (x)`$ such that $`\xi (x_n)=2\pi n`$. The logarithmic function comes to mind first since $`\mathrm{log}\frac{x_n}{a}=n\lambda _+`$. With the correct prefactor $`\alpha =\frac{2\pi }{\lambda _+}`$, we find the desired function

$$\xi (x)=\alpha \mathrm{log}\frac{x}{a}$$
(42)
because $`\xi (x_n)=\alpha \mathrm{log}\frac{x_n}{a}=\frac{2\pi }{\lambda _+}n\lambda _+=2\pi n`$. Substituting $`\xi (x)`$ from (42) into the general form (21), we obtain a manifold described by

$$p=p_0\frac{1}{2}\frac{\mathrm{}\lambda \alpha }{x}\mathrm{cos}(\alpha \mathrm{log}\frac{x}{a})$$
(43)
and displayed in Fig. 5. Replacement manifolds (19) become

$$p_n(x)=p_0+\frac{n\mathrm{}\alpha }{x}.$$
(44)
In the present case, function $`f_n`$ from (34) becomes

$`f_n(x,p,p_0)=(p_0p)x+n\mathrm{}\alpha \mathrm{log}{\displaystyle \frac{x}{a}},`$

and its single SP point,

$$x_{sp}=\frac{n\mathrm{}\alpha }{pp_0}.$$
(45)
For $`p>p_0`$, only manifolds with $`n>0`$ contribute, and for $`p<p_0`$ only $`n<0`$ is allowed. Defining $`j=|n|`$, for $`pp_0`$ we replace the sum over $`n`$ by sum over all positive $`j`$. In both cases, the $`j`$-th manifold gives a semiclassical contribution

$`p|p_0_{j,sc}`$ $`=`$ $`\left({\displaystyle \frac{j\alpha }{2\pi }}\right)^{1/2}{\displaystyle \frac{1}{|pp_0|}}`$ (46)
$`\times `$ $`\mathrm{exp}\left\{isgn(pp_0)\left[j\alpha \mathrm{log}{\displaystyle \frac{j\mathrm{}\alpha }{ae|pp_0|}}{\displaystyle \frac{\pi }{4}}\right]\right\}`$ (47)
with $`j=|n|>0`$. In this case caustics are missing because there is a single SP point $`x_{sp}`$. We have a finite integration limit at zero, but it is separated from $`x_{sp}`$ by an infinite action (equal to the area delimited by lines $`x=0`$, $`p=p_0`$, and the $`n`$-th RM, see Fig. 5). In fact, it can be shown that further terms in asymptotic expansion of $`p|p_0_{j,sc}`$ have the same dependence on $`\mathrm{}`$ and $`pp_0`$ as (46), and only differ in their dependence on $`j\alpha `$. Luckily, these claims can be easily verified by evaluating integral (31) analytically. This oscillatory integral is made convergent by displacing momentum $`p`$ with an imaginary infinitesimal term $`iϵ`$, $`ϵ>0`$,

$`p|p_0_{n,unif}={\displaystyle \frac{1}{2\pi \mathrm{}}}`$ (48)
$`\times `$ $`{\displaystyle _0^{\mathrm{}}}𝑑x\mathrm{exp}\left\{{\displaystyle \frac{i}{\mathrm{}}}\left[(piϵp_0)x+n\mathrm{}\alpha \mathrm{log}{\displaystyle \frac{x}{a}}\right]\right\}.`$ (49)
In this form, the answer is found by rotating the contour about the origin of complex plane by $`\frac{\pi }{2}sgn(pp_0)`$.

Let us explicitly solve the case $`p>p_0`$. Fig. 6 shows a contour in complex plane enclosing a region with no singularities, and as a result the sum of integrals along appropriate parts of the contour is zero, $`I_1+I_R+I_2+I_r=0`$. In the limit $`R\mathrm{}`$, $`r0`$, $`I_1`$ becomes our desired integral (48) and both $`I_R`$ and $`I_r`$ vanish, implying $`I_1=I_2`$. Contour for integral $`I_2`$ is the set of complex points $`z=ix=e^{i\pi /2}x`$, $`x>0`$. Consequently,

$`I_2={\displaystyle \frac{1}{2\pi \mathrm{}}}e^{i\pi /2}\times `$ (50)
$`{\displaystyle _0^{\mathrm{}}}𝑑x\mathrm{exp}\left\{{\displaystyle \frac{i}{\mathrm{}}}[(pp_0iϵ)(ix)+n\mathrm{}\alpha \mathrm{log}{\displaystyle \frac{xe^{i\pi /2}}{a}}]\right\}.`$ (51)
Transforming to a dimensionless variable $`y=\frac{(pp_0iϵ)x}{\mathrm{}}`$, the integral becomes

$`I_2`$ $`=`$ $`{\displaystyle \frac{1}{2\pi \mathrm{}}}\mathrm{exp}\left\{(n\alpha i){\displaystyle \frac{\pi }{2}}in\alpha \mathrm{log}\left[(pp_0iϵ)a/\mathrm{}\right]\right\}`$ (52)
$`\times `$ $`{\displaystyle \frac{\mathrm{}}{pp_0iϵ}}{\displaystyle _0^{\mathrm{}}}𝑑ye^yy^{in\alpha }.`$ (53)
Recognizing the remaining integral as a gamma function of complex argument $`1+in\alpha `$,

$`I_2={\displaystyle \frac{\mathrm{\Gamma }(1+in\alpha )}{2\pi (pp_0iϵ)}}`$ (54)
$`\times \mathrm{exp}\left\{(n\alpha i){\displaystyle \frac{\pi }{2}}in\alpha \mathrm{log}\left[(pp_0iϵ)a/\mathrm{}\right]\right\}.`$ (55)
For $`p<p_0`$, we need to close the contour in the upper half-plane, but otherwise the procedure is analogous. Combining the two cases, we obtain a single analytic expression for the uniform transformation element,

$`p|p_0`$ $`_{n,unif}={\displaystyle \frac{\mathrm{\Gamma }(1+in\alpha )}{2\pi (|pp_0|isϵ)}}`$ (56)
$`\times `$ $`\mathrm{exp}\left[(n\alpha i){\displaystyle \frac{\pi }{2}}sin\alpha \mathrm{log}{\displaystyle \frac{(|pp_0|isϵ)a}{\mathrm{}}}\right]`$ (57)
where $`s=sgn(pp_0)`$. As promised, the uniform and semiclassical forms have identical dependence on $`\mathrm{}`$ and $`pp_0`$. The only difference lies in the dependence on $`n\alpha `$, in particular the semiclassical form is obtained if we keep only the first term in the asymptotic expansion of $`\mathrm{\Gamma }(1+in\alpha )`$ in $`p|p_0_{n,unif}`$. However, for our conservative choice of $`\alpha 33`$ the agreement between $`p|p_0_{n,sc}`$ and $`p|p_0_{n,unif}`$ is such that they may be used interchangeably for any practical purposes. Finally, in Fig. 7, $`p|p_0`$ is applied to an initial Gaussian wave packet centered around $`X_0`$, $`P_0`$, namely

$`\varphi ^i(p_0)=\left({\displaystyle \frac{\pi \mathrm{}^2}{\sigma ^2}}\right)^{1/4}\mathrm{exp}\left[{\displaystyle \frac{i}{\mathrm{}}}(P_0p_0)X_0{\displaystyle \frac{\sigma ^2(p_0P_0)^2}{2\mathrm{}^2}}\right].`$

The final position wave function is

$`\psi _{sc}^f(x)=(\pi \sigma ^2)^{1/4}\mathrm{exp}\left[{\displaystyle \frac{i}{\mathrm{}}}F_2(x,P_0){\displaystyle \frac{(xX_0)^2}{2\sigma ^2}}\right]`$

and the uniform momentum wave function,

$`\varphi _{unif}^f(p)=J_0\left({\displaystyle \frac{\lambda }{2}}\right)\varphi ^i(p)`$ (58)
$`+`$ $`{\displaystyle \underset{j=1}{\overset{\mathrm{}}{}}}J_j\left({\displaystyle \frac{\lambda }{2}}\right)i^\lambda p|p_0_{j,sc}{\displaystyle \frac{(2\pi \mathrm{})^{1/2}}{(\pi \sigma ^2)^{1/4}}}\mathrm{exp}\left[{\displaystyle \frac{(x_{sp,j}x_0)^2}{2\sigma ^2}}\right]`$ (59)
with $`x_{sp,j}=j\mathrm{}\alpha /|pp_0|`$ and $`p|p_0_{j,sc}`$ given by (46). The large oscillation near $`p=p_0`$ in Fig. 7 corresponds to the zeroth RM contribution $`\varphi ^i(p)`$, the smaller wavelet to the right corresponds to the remaining RMs. The figure confirms the excellent accord between RM method and numerical IVR evaluation by FFT. We have thus succeeded in uniformizing something with many of the properties of homoclinic oscillations near an unstable periodic orbit using replacement-manifold “technology”.

## VI Comparison with the stationary-phase method

It is instructive to check explicitly if the expansion in terms of RMs agrees with semiclassical method applied to the original manifold when areas of loops are larger than $`\mathrm{}`$ ($`\lambda >1`$). Let us choose an analytically solvable example with $`\xi (x)=\alpha x`$. Corresponding manifold

$$p(x)=p_0\frac{1}{2}\mathrm{}\lambda \alpha \mathrm{sin}\alpha x$$
(60)
has another advantage compared to manifolds in Figures 3 and 5. Unlike those, for $`\lambda 1`$ and small enough $`pp_0`$, manifold (60) has all caustics in a safe distance. The RMs are horizontal lines $`p_n=p_0+n\mathrm{}\alpha `$, independent of $`x`$ and corresponding transformation elements can be evaluated exactly as $`p|p_0_n=\delta (pp_0n\mathrm{}\alpha )`$. Using (29), the uniform expression for $`p|p_0`$ is

$$p|p_0_{unif}=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}J_n\left(\frac{\lambda }{2}\right)i^n\delta (pp_0n\mathrm{}\alpha ),$$
(61)
so the wave function is determined by its Fourier coefficients $`c_n(\lambda )=J_n(\lambda /2)i^n`$, which is natural since manifold (60) is exactly periodic. Moreover, each RM contributes only to a single momentum ($`p_n`$). Put differently, we have calculated the wave function for any specific momentum to all orders in $`\lambda `$. Now let us find the semiclassical form of $`p|p_0`$. For a given momentum $`p`$ there is an infinite number of SP points which occur in pairs $`(x_n,y_n)`$,

$$x_n=x_0+\frac{2\pi }{\alpha }n,y_n=\frac{\pi }{\alpha }x_0+\frac{2\pi }{\alpha }n,$$
(62)
where $`x_0=\frac{1}{\alpha }\mathrm{arcsin}\frac{2(pp_0)}{\lambda \mathrm{}\alpha }`$. Analogously to (34), the semiclassical transformation element for manifold (60) is

$`p|p_0_{sc}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi \mathrm{}}}}{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}\left|{\displaystyle \frac{^2S}{x^2}}\right|_{x=x_n}^{1/2}`$ (63)
$`\times `$ $`\left\{\mathrm{exp}\left[{\displaystyle \frac{i}{\mathrm{}}}S(x_n)i{\displaystyle \frac{\pi }{4}}\right]+\mathrm{exp}\left[{\displaystyle \frac{i}{\mathrm{}}}S(y_n)+i{\displaystyle \frac{\pi }{4}}\right]\right\}`$ (64)
with

$`S(x,p,p_0)=px+F_2(x,p_0)=(p_0p)x+{\displaystyle \frac{1}{2}}\mathrm{}\lambda \mathrm{cos}\alpha x.`$

Using (62), (63) and defining $`r=\frac{pp_0}{\mathrm{}\alpha }`$, we get

$$p|p_0_{sc}=\frac{1}{\mathrm{}\alpha }F(\lambda ,r)\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}e^{2\pi inr}$$
(65)
where

$`F(\lambda ,r)`$ $`=\sqrt{{\displaystyle \frac{2}{\pi }}}e^{ir\pi /2}`$ (66)
$`\times `$ $`{\displaystyle \frac{\mathrm{cos}\left\{r\left(\frac{\pi }{2}+\mathrm{arcsin}\frac{2r}{\lambda }\right)+\left[\left(\frac{\lambda }{2}\right)^2r^2\right]^{1/2}\frac{\pi }{4}\right\}}{\left[(\frac{\lambda }{2})^2r^2\right]^{1/4}}}.`$ (67)
Employing the Poisson summation formula

$`{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}e^{2\pi inr}={\displaystyle \underset{m=\mathrm{}}{\overset{\mathrm{}}{}}}\delta (mr)`$

and reverting to $`p=p_0+r\mathrm{}\alpha `$, we obtain

$`p|p_0_{sc}`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{}\alpha }}F(\lambda ,r){\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}\delta (nr)`$ (68)
$`=`$ $`{\displaystyle \frac{1}{\mathrm{}\alpha }}{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}F(\lambda ,n)\delta (nr)`$ (69)
$`=`$ $`{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}F(\lambda ,n)\delta (Pp_0n\mathrm{}\alpha ).`$ (70)
Comparing (68) with (61), we see that the semiclassical and uniform versions of $`p|p_0`$ will be asymptotically equal if $`F(\lambda ,n)J_n(\lambda /2)i^n`$ for large $`\lambda `$. That is indeed true since for $`\lambda n^2`$,

$`F(\lambda ,n)`$ $``$ $`\left({\displaystyle \frac{2}{\pi \frac{\lambda }{2}}}\right)^{1/2}e^{in\pi /2}\mathrm{cos}\left({\displaystyle \frac{\lambda }{2}}{\displaystyle \frac{n\pi }{2}}{\displaystyle \frac{\pi }{4}}\right)`$ (71)
$``$ $`J_n\left({\displaystyle \frac{\lambda }{2}}\right)i^n`$ (72)
Table I shows that for $`\lambda =10`$ semiclassical and RM values of $`a_n`$ differ by less than $`0.02`$ up to $`n=3`$. In the opposite limit, for $`\lambda <1`$, even $`a_0`$ evaluated semiclassically is completely off.

## VII Scattering from a corrugated wall

Let us apply the replacement-manifold method to a physical problem generating loop structure in phase space. Scattering of a plane wave by a corrugated wall has exactly such a property. This system has been used to model elastic scattering of atoms by solid surfaces (for review see ). Consider a two-dimensional coordinate space divided into two parts by a periodically curved boundary consisting of the set of points with coordinates related by

$$\stackrel{~}{y}(x)=\frac{\epsilon }{\beta }\mathrm{sin}\beta x$$
(73)
where $`\beta `$ gives the spatial frequency of corrugation and $`\epsilon `$ the maximum slope $`\frac{d\stackrel{~}{y}}{dx}`$ of the wall. A plane wave with momentum $`\mathrm{}𝐤_1`$ is incident from the upper half-plane at an angle $`\alpha _1`$ from the $`y`$ axis, so the incident wave vector is $`𝐤_1=k(\mathrm{sin}\alpha _1,\mathrm{cos}\alpha _1).`$ Classically, the wave reflects specularly from the curved boundary. The reflected rays are shown in Fig. 8.

Notice also the clearly visible caustics. The Poincaré surface of section in Fig. 9 displays dependence of momentum component $`p_{2x}`$ of reflected rays on coordinate $`x`$ at a given distance $`y_2`$ from $`x`$ axis. We can see familiar loops with constant area and predict the failure of semiclassical approximation when this area gets smaller than $`\mathrm{}`$. Below we uniformize the semiclassical solution using RMs and compare the outcome with the exact and primitive semiclassical results.

There exist several approaches to solve this scattering problem exactly. Garcia and Cabrera  have thoroughly compared merits of various methods. The main issue is the solution’s convergence, which improves with decreasing corrugation parameter $`ϵ`$. Luckily, this regime overlaps with the small-loop limit in which we are interested. We can therefore use the most straightforward approach based on the so-called Rayleigh hypothesis  which has been shown to be valid for $`ϵ<0.448`$ by Millar . Details of the solution are provided in the Appendix. A probability density plot of the wave function is shown in Fig. 10.

The whole scattering problem can be formulated using an analogue to the smooth-potential Lippmann-Schwinger equation. Invoking the Green’s theorem, it can be shown that the total wave function satisfies

$$|\psi _{tot}=|\psi _{inc}+\widehat{G}_0_{boundary}𝑑x|𝐫\widehat{n}(x)\psi _{tot}(𝐫)$$
(74)
where $`\widehat{G}_0`$ is the free-space Green’s operator and $`\widehat{n}(x)`$ is a normal unit vector at the boundary pointing into free space. Now we could use periodicity of the wall, expand $`\psi _{tot}(𝐫)`$ in its Fourier modes,

$`\psi _{tot}(𝐫)`$ $`=`$ $`\psi _{inc}(𝐫)+\psi _{scat}(𝐫)`$ (75)
$`=`$ $`{\displaystyle \frac{1}{2\pi }}e^{i𝐤_1𝐫}+{\displaystyle \frac{1}{2\pi }}{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}a_ne^{i𝐤_{2,n}𝐫},`$ (76)
and recover the linear system (A8) for coefficients $`a_n`$. Our goal, however, is to make a connection with a semiclassical picture in the Poincaré surface of section (Fig. 9). That figure implicitly assumes that classically only a single scattering event takes place before a ray leaves the wall permanently. If corrugation is deep enough ($`\epsilon >\epsilon _{\mathrm{max}}`$, where $`\epsilon _{\mathrm{max}}<1`$), multiple scattering will occur for any incident angle $`\alpha _1`$. The larger $`\epsilon `$ gets (within the range $`0<\epsilon <\epsilon _{\mathrm{max}}`$), the smaller the maximum allowed incident angle $`\alpha _1`$.

Let us consider a wall with very shallow corrugation ($`\epsilon =0.2)`$ and a small angle of incidence ($`\alpha _1=0.2)`$, a situation which classically allows only a single reflection. We also choose the incident-wave length $`\lambda =\frac{2\pi }{k}`$ to be $`0.3`$ times the period of corrugation. The surprising result we shall obtain below is that although there exist seven real Bragg peaks in the scattered wave, the semiclassics break down while the single-scattering approximation using replacement manifolds works with excellent accuracy.

The simple semiclassical wave function may be evaluated e.g. by tracing of individual rays and employing the Van Vleck propagator . The action must be adjusted by correct Maslov indices, corresponding to reflection from a hard wall and to passage through caustics . A probability density plot of the semiclassical wave function is shown in Fig. 11. Note the caustics, which clearly separate regions with one, three, and five contributions to the scattered wave. Between caustics, solution looks qualitatively the same as the exact quantum analogue in Fig. 10. By looking at a surface of section at $`y_2=4\pi /\beta `$ (where there are already at least three contributions for any $`x`$) in Fig. 9, we expect that semiclassical and quantum wave functions should disagree everywhere. This guess is confirmed in Fig. 12.

Now let us find a uniform expression for the scattered wave. The Lippmann-Schwinger equation (74) is an integral equation, and it is not immediately obvious how it could be used to find $`\psi _{scat}`$ from $`\psi _{inc}`$. We circumvent this problem by proceeding in four steps: first we canonically transform the incident wave function to a new coordinate system $`𝐫^{}`$ in which the wall becomes straight. We solve the scattering problem in these coordinates, since there the semiclassical approach is exact. Then we transform back to the original coordinate system in which we propagate the wave using the free-space Green’s function to obtain our final answer. Semiclassical approach is used for each separate step but their combination is evaluated exactly.

The wall becomes flat if we apply a canonical transformation

$`x^{}=x,y^{}=y{\displaystyle \frac{\epsilon }{\beta }}\mathrm{sin}\beta x`$

generated by

$$F_2(𝐫^{},𝐤)=k_xx^{}+k_y(y^{}+\frac{\epsilon }{\beta }\mathrm{sin}\beta x^{}).$$
(77)
The incident wave function $`\psi _{inc}(𝐤_1)`$ can be transformed into new coordinates using

$$𝐫^{}|𝐤_{sc}=\frac{1}{2\pi }e^{iF_2(𝐫^{},𝐤)}$$
(78)
In primed coordinates, the semiclassical solution of scattering is exact and “angle of incidence” equals “angle of reflection”. Therefore $`\widehat{n}^{}\psi _{tot}(𝐫^{})`$ may be replaced by $`2\widehat{n}^{}\psi _{inc}(𝐫^{})`$ in equation (74). The non-diagonal part of the semiclassical scattering matrix becomes

$$\widehat{T}=\widehat{S}\widehat{1}=2𝑑x^{}|𝐫^{}\widehat{n}^{}^{}𝐫^{}||_{y^{}=0}$$
(79)
where $`\widehat{n}^{}=\widehat{y}^{}`$. Transformation back to original coordinates is achieved by a complex conjugate of (78). The free-space propagation is accomplished by a semiclassical Green’s function, exact in mixed representation

$$k_{2x},y_2|\widehat{G}_0^{sc}|𝐤=\frac{\delta (k_{2x}k_x)e^{ik_yy_2}}{\sqrt{2\pi }(k^2k_2^2+iϵ)}$$
(80)
Putting all four pieces together, we obtain an expression for the scattered wave on a surface of section,

$``$ $`k_{2x},y_2|\widehat{G}_0^{sc}\widehat{T}|𝐤_1`$ (81)
$`=`$ $`2{\displaystyle d^2k𝑑x^{}k_{2x},y_2|\widehat{G}_0^{sc}|𝐤𝐤|𝐫^{}_y^{}𝐫^{}|𝐤_{sc}}|_{y^{}=0}`$ (82)
$`=`$ $`{\displaystyle \frac{2ik_{1y}}{(2\pi )^{5/2}}}{\displaystyle d^2k\frac{\delta (k_{2x}k_x)e^{ik_yy_2}}{\sqrt{2\pi }(k^2k_2^2+iϵ)}}`$ (83)
$`\times `$ $`{\displaystyle 𝑑x^{}\mathrm{exp}[i(k_{1x}k_x)x^{}+i(k_{1y}k_y)\frac{\epsilon }{\beta }\mathrm{sin}\beta x^{}]}`$ (84)
Before we proceed, it should be noted that the outcome of this four-step process, an integral representation of the scattered wave, can be viewed as a continuous superposition of waves emanating from sources along the boundary, with strength proportional to the normal derivative of the incident wave. It can be shown, in fact, that this process is equivalent to a generalized Kirchhoff diffraction method .

Returning to expression (81), we recognize that the exponent in integral over $`x^{}`$ corresponds to a classical manifold

$$k_x=k_{1x}+(k_{1y}k_y)\epsilon \mathrm{cos}\beta x^{}$$
(85)
which is “begging” to be replaced by partial manifolds $`k_{x,n}=k_{1x}+n\beta x^{}`$ with weights $`J_n(\lambda )`$ because

$`|\lambda |=|k_{1y}k_y|{\displaystyle \frac{\epsilon }{\beta }}{\displaystyle \frac{2k\epsilon }{\beta }}<1`$

for the classically allowed momenta. The integral over $`x^{}`$ is then simple to evaluate and is equal to

$`{\displaystyle 𝑑x^{}}`$ $`{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}J_n\left(\lambda \right)\mathrm{exp}[i(k_{1x}+n\beta k_x)x^{}]`$ (86)
$`=`$ $`2\pi {\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}J_n(\lambda )\delta (k_{1x}+n\beta k_x)`$ (87)
Using this result and evaluating the trivial integral over $`k_x`$, expression (81) becomes

$``$ $`k_{2x},y_2|\widehat{G}_0^{sc}\widehat{T}|𝐤_1={\displaystyle \frac{2ik_{1y}}{(2\pi )^{3/2}}}`$ (88)
$`\times `$ $`{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle 𝑑k_y\frac{e^{ik_yy_2}}{\sqrt{2\pi }(k^2k_2^2+iϵ)}J_n(\lambda )\delta (k_{2x}k_{2x,n})}`$ (89)
where $`k_{2x,n}=k_{1x}+n\beta `$. Integral over $`k_y`$ picks up a pole at $`k_y=k_{2y,n}=\sqrt{k^2k_{2x,n}^2}`$ and the final answer is

$`k_{2x},y_2`$ $`|`$ $`\widehat{G}_0^{sc}\widehat{T}|𝐤_1`$ (90)
$`=`$ $`{\displaystyle \frac{k_{1y}}{\sqrt{2\pi }}}{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}e^{ik_{2y,n}y_2}{\displaystyle \frac{J_n(\lambda _n)\delta (k_{2x}k_{2x,n})}{k_{2y,n}}}`$ (91)
where $`\lambda _n=(k_{1y}k_{2y,n})\frac{\epsilon }{\beta }`$. The scattered wave is thus a superposition of traveling and evanescent waves with wavevectors $`𝐤_{2,n}`$ with coefficients

$$a_{n,unif}=\frac{k_{1y}}{k_{2y,n}}J_n(\lambda _n).$$
(92)
These turn out to be very close to the exact quantum coefficients obtained in the Appendix, but here we have avoided the solution of a linear system (A8).

The uniform expression (92) for coefficients of expansion (76) is surprisingly accurate. A two-dimensional density plot in Fig. 10 does not reveal any difference from the exact quantum solution while an analogous semiclassical plot clearly shows the caustics (see Fig. 11). Even if we look at the Poincaré surface of section in a region with many classical caustics, the quantum and uniform solutions agree while semiclassical solution fails miserably (see Fig. 12). The uniform solution agrees with one found by Garibaldi et al. and is somewhat more accurate than the same result without the $`k_{1y}/k_{2y,n}`$ prefactor obtained by Hubbard and Miller  by SCP approximation and by Garibaldi et al. . As Garibaldi et al. (who obtain three different solutions differing only by prefactors) point out, these prefactors “should not be taken too seriously.” We agree: for traveling modes, they do not cause large errors, and since all these solutions neglect multiple scattering, we cannot expect high accuracy of the already small coefficients of evanescent modes.

As promised in Section IV, we have obtained almost identical result using RMs as Hubbard and Miller with SCP approximation . They found their answer for a hard wall as a limit of the solution for a smooth potential due to Lennard-Jones and Devonshire. Although we were lucky in this case to obtain exactly solvable integrals for RMs, we can quite generally expect that the replacement manifolds will be much smoother than the original one, enabling their simple evaluation by SP method.

We discuss the specific choice of parameters. Selecting $`\epsilon 1`$, $`\alpha 1`$ was necessary for the validity of the single-scattering approximation. At first, it appears that the simple semiclassical approximation should be accurate for this regime of very shallow corrugation, since scattering from a flat wall has an exact semiclassical answer. The correct criterion, however, is based on the area of loops in phase space (see Fig. 9). In this case, for small $`\epsilon `$ and $`\alpha `$, the height of loops is equal to the momentum kick in region of maximum slope, approximately $`p\epsilon `$ and the width is half the period of corrugation, $`\frac{\pi }{\beta }`$. The loop area is smaller than $`p\epsilon \frac{\pi }{\beta }=`$ $`\pi kd\mathrm{}`$ ($`d`$ is the depth of corrugation) and has to be compared to $`\mathrm{}`$. We arrive at a surprising result that the validity of semiclassical approximation has nothing to do with the periodicity of corrugation, but only depends on the product $`kd`$ of wave vector and corrugation depth. If this product is small, semiclassical approximation breaks down. In the opposite case of deep corrugation, semiclassical approximation works, although we might have to take into account classical trajectories corresponding to multiple scattering.

## VIII Conclusion

We have successfully uniformized several intimidating situations which can arise when repeated small areas or loops arise from enclosing classical manifolds using the replacement-manifold idea. This notion can rejuvenate physical intuition about the quantum wave function and also be a very convenient approach to an accurate result. Two physically motivated cases were considered (the homoclinic-like oscillation and the corrugated wall) along with several other contrived cases to test the method. The present method should work in more complicated time-dependent problems if used with Miller’s IVR applied to Van Vleck’s semiclassical propagator. After replacement manifolds are identified, the perturbative expansion contributions to the wave function could be found by SP approximation applied to these manifolds.

Future work along these lines includes the attempt to uniformize a full homoclinic tangle. It remains to be seen how convenient and generally applicable the replacement manifold idea is. However, it seems clear that we now have a new method which simplifies many problems while reinstating an intuitive foundation.

## IX Acknowledgements

This research was supported by the National Science Foundation under Grant No. CHE9610501 and by the Institute for Theoretical Atomic and Molecular Physics. One of us (J. V.) would like to acknowledge helpful discussions with A. Mody and A. Barnett.

## A Exact solution of scattering from a hard sinusoidal wall

From conservation of energy, the scattered and incident waves must have the same magnitude of momentum. Due to the periodicity of the wall in the $`x`$ direction, the $`x`$-component of momentum is constrained by the Bragg condition

$$k_{2x,n}=k_{1x}+n\beta ,nZ.$$
(A1)
Corresponding $`y`$-component of momentum $`k_{2y,n}=\sqrt{k^2k_{2x,n}^2}`$ will be real only if $`|k_{1x}+n\beta |<k`$. Otherwise the reflected wave becomes evanescent: it decays exponentially with distance from $`x`$-axis because $`k_{2y}`$ is imaginary. To complete the solution, coefficients $`a_n`$ in the momentum Green’s function

$`𝐤_{2x}|𝐤_{1x}={\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}a_n\delta (k_{2x}k_{2x,n})\delta (k_{2y}k_{2y,n})`$

must be found. This is accomplished by imposing the Dirichlet boundary condition in coordinate space. Namely, the total wave function $`\psi _{tot}(𝐫)=\psi _{inc}(𝐫)+\psi _{scat}(𝐫)`$ has to vanish along the corrugated wall. Using

$`\psi _{inc}(𝐫)`$ $`=`$ $`𝐫|𝐤_1={\displaystyle \frac{1}{2\pi }}e^{i𝐤_1𝐫},`$ (A2)
$`\psi _{scat}(𝐫)`$ $`=`$ $`{\displaystyle d^2k_2𝐫|𝐤_2𝐤_2|𝐤_1}={\displaystyle \frac{1}{2\pi }}{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}a_ne^{i𝐤_{2,n}𝐫},`$ (A3)
definition of boundary (73) and (A1) for $`k_{2x,n}`$, the boundary condition becomes

$`\mathrm{exp}\left[i(k_{1x}x+k_{1y}{\displaystyle \frac{\epsilon }{\beta }}\mathrm{sin}\beta x)\right]`$ (A4)
$`+`$ $`{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}a_n\mathrm{exp}\left[i(k_{1x}+n\beta )x+ik_{2y,n}{\displaystyle \frac{\epsilon }{\beta }}\mathrm{sin}\beta x\right]=0.`$ (A5)
Recalling the identity $`e^{i\lambda \mathrm{sin}\xi (x)}=_{m=\mathrm{}}^{\mathrm{}}J_m(\lambda )e^{imx}`$ and denoting $`\lambda _1=\frac{\epsilon }{\beta }k_{1y}`$, $`\lambda _{2,n}=\frac{\epsilon }{\beta }k_{2y,n}`$, $`\theta =\beta x`$, the boundary condition may be written as

$$\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}\left[J_m(\lambda _1)e^{im\theta }+\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}a_nJ_m(\lambda _{2,n})e^{i(n+m)\theta }\right]=0.$$
(A6)
Replacing $`m`$ by $`pn+m`$ in the second term and by $`pm`$ in the first term in the brackets in (A6), we get

$$\underset{p=\mathrm{}}{\overset{\mathrm{}}{}}e^{ip\theta }\left[J_p(\lambda _1)+\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}a_nJ_{pn}(\lambda _{2,n})\right]=0.$$
(A7)
This equation will be satisfied for all $`x=\frac{\theta }{\beta }`$ only if coefficients of $`e^{ip\theta }`$ on the left-hand side vanish,

$$J_p(\lambda _1)+\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}a_nJ_{pn}(\lambda _{2,n})=0,pZ.$$
(A8)
We recognize (A8) as an infinite system of linear equations for coefficients $`a_n`$. In matrix form, $`𝐀𝐚=𝐛`$, with $`A_{ij}J_{ij}\left(\frac{\epsilon }{\beta }\sqrt{k^2(k_{1x}+j\beta )^2}\right)`$and $`b_iJ_i(\frac{\epsilon }{\beta }k_{1y})`$. Although this system cannot be solved analytically, expressing the solution in this form is a great improvement compared to what can be done for a general Dirichlet problem satisfying time-independent Schrödinger equation. The system (A8) can be approximately evaluated numerically by restricting $`n`$ and $`p`$ to lie between limiting values $`\pm n_{\mathrm{max}}`$.


text:
# SPATIAL DIFFUSION OF X-RAY EMISSION LINES IN THE M87 COOLING FLOW; EVIDENCE FOR ABSORPTION

## 1 DIFFUSION OF X-RAY LINE PHOTONS IN M87

This estimate of the diffusion of X-ray line photons in M87 has been motivated by the exceptionally low abundances of iron, silicon and oxygen within $`5`$ kpc of the center of M87 as recently observed by Böhringer et al (2000) with XMM-Newton. We explore the possibility that these abundance variations are artifacts of the outward diffusion of optically thick X-ray lines as suggested by Gil’franov et al. (1987), Tawara et al. (1997), Shigeyama (1998) and Böhringer et al. (2000). We find that the expected outward migration of line photons does occur, but is insufficient to account for the central photon depletion observed. Instead, an additional source of continuous absorption is indicated, as suggested by Buote (2000), and we speculate about its physical nature. Before discussing our results, we begin with a review of the properties of the M87 cooling flow and describe our line transfer calculation.

## 2 VARIATION OF GAS TEMPERATURE AND DENSITY

The electron density distribution in the M87 cooling flow is taken from Nulsen and Böhringer (1995), corrected to distance 17 Mpc (1’ = 4.945 kpc):

$$n_{e,g}=n_0[1+(r/r_o)^p]^1$$
(1)
with $`n_0=0.120`$ cm-3, $`p=1.18`$ and $`r_o=0.5r_e`$ where $`r_e=8.625`$ kpc is the effective radius. According to Böhringer et al. (2000) these observations underestimate the central density which should instead be represented by a $`\beta `$ model fit:

$$n_{e,\beta }=n_1[1+(r/r_\beta )^2]^{3\beta /2}$$
(2)
with $`n_1=0.35`$ cm-3, $`\beta =0.47`$ and $`r_\beta =1.97`$ kpc. This inner density distribution is based on ROSAT HRI observations reported in Böhringer (1999). Since the central X-ray emission in M87 is known to be partly non-thermal (and time variable), this latter density determination may be uncertain. Nevertheless, in the following we assume that $`n_e(r)`$ is given by the maximum of $`n_{e,g}(r)`$ and $`n_{e,\beta }(r)`$ which intersect at 11.9 kpc.

The (single-phase) gas temperature distribution is found by combining the results of Böhringer et al. (2000) within 65 kpc with more distant temperatures measured by Nulsen & Böhringer (1995). These results can be fit with the fitting function:

$$T=T_o\left(\frac{r_m}{r_{ot}}\right)\left[\frac{r_m}{(r+r_{ot})}+\left(\frac{r}{r_m}\right)^q\right]^1$$
(3)
with $`T_o=1.5\times 10^7`$ K, $`r_m=5.695`$, $`r_{ot}=1.53`$ and $`q=0.22`$; all radii are in units of $`r_e`$. This temperature profile rises from the center toward a broad maximum, similar to thermal profiles in other well-observed bright elliptical galaxy cooling flows (Brighenti & Mathews 1997). The location of the temperature maximum in M87 is at an unusually large radius presumably due to the relatively high ambient temperature (2 - 3 keV) in the Virgo cluster. The central cooling of the gas can be understood as the mixing of hot inflowing gas from the Virgo cluster with gas lost from evolving stars in M87 which orbit with a lower virial temperature than that of the surrounding cluster (Brighenti & Mathews 1999a). The characteristic deep central temperature minimum in M87 indicates that most of the gas in M87 has not been substantially heated by the jet and radio source.

## 3 ABUNDANCES, LINE EMISSION AND ABSORPTION

We assume meteoritic solar abundances with $`\mathrm{log}(a_E)\mathrm{log}(n_E/n_H)=3.07`$ and $`4.45`$ for O and Si respectively. The radial abundance variation of each element in the hot gas relative to solar meteoritic is represented with $`A_E(r)`$. The line emissivity is assumed to vary with local temperature according to the emission coefficients $`\mathrm{\Lambda }_{\mathrm{}}`$ determined by Landini & Monsignori Fossi (L-MF; 1990):

$$\epsilon _{\mathrm{}}=A_E(r)n_en_H\alpha \mathrm{\Lambda }_{\mathrm{}}(T)(\mathrm{erg}\mathrm{cm}^3\mathrm{s}^1)$$

where $`\alpha `$ is a factor of order unity that corrects L-MF abundances to solar meteoritic and $`n_H=n_e(43\mu )/(2+\mu )`$ is the local proton density in a plasma with molecular weight $`\mu =0.61`$.

The emission lines that we consider arise from the ground state and are Doppler broadened by thermal and turbulent velocities. The absorption coefficient is

$$\kappa (x)=n_Hx_ia_EA_E\left(\frac{\pi e^2}{m_ec}f_{gu}\frac{1}{\pi ^{1/2}}\frac{1}{\mathrm{\Delta }\nu _D}\right)e^{x^2}(\mathrm{cm}^1).$$

Here $`a_EA_E(r)`$ is the local abundance of the line-emitting element and $`x_i(r)`$ is the fractional abundance of the line-emitting ion interpolated from the results of Sutherland & Dopita (1993) for collisional ionization equilibrium. Photons are assumed to be emitted in a Gaussian profile centered on the line energy $`E_o=h\nu _o`$ and

$$x=\frac{\nu \nu _o}{\mathrm{\Delta }\nu _D}$$

is the dimensionless line frequency with

$$\mathrm{\Delta }\nu _D=\frac{\nu _o}{c}\left(\frac{2kT}{A_{Ewt}m_p}+v_{turb}^2\right)^{1/2}$$
(4)
where $`A_{Ewt}`$ is the atomic weight of element E. The line center optical depth in the radial direction is found by integrating $`d\tau _o=\kappa (0)dr`$.

We are taking $`x_i(T)`$ and $`\mathrm{\Lambda }_{\mathrm{}}(T)`$ from two different sources. However, the error from this possible inconsistency is unlikely to be large since collisional ionization fractions $`x_i(T)`$ found by most authors, including L-MF, are in good agreement with those of Arnaud & Rothenflug (1985).

## 4 MONTE CARLO PROCEDURE

To estimate the diffusion of X-ray line photons, we represent the inner cooling flow of M87 with $`N=50`$ radial zones, each having different but individually uniform temperature and density. The zone spacing $`r_n`$ is chosen so that the line luminosity is the same from each zone. The Monte Carlo process begins by choosing the zone of origin of each photon based on equal weighting among the zones. The exact birth radius $`r_b`$ of each photon within the chosen zone is then found from

$$r_b=[r_n^3+(r_{n+1}^3r_n^3)]^{1/3}$$

where $``$ is a random number between 0 and 1. The initial direction of the photon is $`w=\mathrm{cos}\theta =21`$ where $``$ is another random number and $`\theta `$ is the angle between the photon velocity and the direction to the center of the cooling flow. The initial frequency $`x`$ is chosen by random selection from a Gaussian probability density. The distance that the photon moves from its origin $`r_o`$ to the local zone boundary $`r_z`$ is

$$\mathrm{}=wr_o+\sigma [r_z^2r_o^2(1w^2)]^{1/2}.$$
(5)
Here $`\sigma =1`$ and $`r_z=r_{n+1}`$ for outward moving photons ($`w<w_{cr}`$) and $`\sigma =1`$ and $`r_z=r_n`$ for inward moving photons ($`w>w_{cr}`$) where $`w_{cr}=[1(r_n/r_o)^2]^{1/2}`$. In the zone of photon birth $`r_o=r_b`$. If the photon continues to move through several zones without scattering, $`w`$ and $`r_o`$ remain constant and Equation (5) can be used to find the total path length moved $`L`$. The distance $`\mathrm{}`$ traversed by the photon across each zone must be compared to the free path in the zone:

$$\mathrm{}_s=\mathrm{ln}()/\kappa _n(x_n)$$

where $`\kappa _n(x_n)`$ is the absorption coefficient of the photon in zone $`n`$. Since $`\mathrm{\Delta }\nu _D`$ varies from zone to zone it is necessary to rescale the frequency to the local zone $`n`$ ($`r_n<r<r_{n+1}`$): $`x_n(r)=x_n(r_b)[\mathrm{\Delta }\nu _D(r_b)/\mathrm{\Delta }\nu _D(r)]^{1/2}`$. If $`\mathrm{}_s>\mathrm{}`$, the photon moves into the next zone, but if $`\mathrm{}_s<\mathrm{}`$ the photon is scattered at radius

$$r_s=(r_o^2+L_s^22r_oL_sw)^{1/2}$$

where $`L_s=\mathrm{}_n+\mathrm{}_s`$. The photon is then scattered into a new random direction beginning a new trajectory with $`r_o=r_s`$, assuming complete redistribution in the Gaussian line profile. The scattering is conservative so all photons must ultimately flow out beyond the outermost zone. When the photon reaches the outer boundary of the gas $`r_t=r_N`$, where the optical depths are very small, it is assumed to escape freely.

## 5 RESULTS AND DISCUSSION

We describe transfer calculations for two lines: SiXIV 2.003 keV (L$`\alpha `$) and OVIII 0.652 keV (L$`\alpha `$), both having f-values $`f_{gu}=0.420`$. Since these resonance lines are among the most optically thick lines of Si and O at M87 gas temperatures, the apparent abundance distribution inferred from each of these lines – subsequent to spatial diffusion – must correspond to the largest possible central abundance depletion. For each line-emitting element the apparent radial abundance profiles observed with XMM-Newton can be fit with

$$A_E^{(ob)}(r)=A_{ob}\left(\frac{r_{mob}}{r_{otob}}\right)\left[\frac{r_{mob}}{(r+r_{otob})}+\left(\frac{r}{r_{mob}}\right)^{qob}\right]^1.$$
(6)
The apparent (meteoritic) abundance dependence of Si and O are approximated with: ($`A_{ob}`$, $`r_{mob}`$, $`r_{otob}`$, $`qob`$) = (0.1 , 7.4, 0.34, 0.7) and (0.05, 8.0, 0.5, 0.94) respectively; $`r_{mob}`$ and $`r_{otob}`$ are in kpc. The line center optical depths in the three lines evaluated with $`A_E^{(ob)}(r)`$ (Fig. 1) are smaller than those suggested by Böhringer et al (2000) in part because we use the higher temperatures observed instead of their isothermal approximation. If the observed abundance maxima are due to outward photon diffusion by conservative scattering before escape, the total X-ray line luminosity within some larger radius (where $`\tau 1`$) is unaffected by the radiative transfer.

According to this line diffusion hypothesis, the true radial abundance profiles are expected to decrease monotonically as we represent with

$$A_E(r)=\frac{A_o}{1+(r/r_E)^c}+A_{cf}.$$
(7)
The first term accounts for the local enrichment by stellar mass loss and SNIa explosions, and the second term represents the abundance of cluster cooling flow gas flowing into M87, but the parameterized fit is not intended to exactly quantify these two enrichment sources. The coefficients in Equation (7) are chosen so that the total line luminosities within 50 kpc are very close to luminosities based on $`A_E^{(ob)}`$ (see below): ($`A_o`$, $`r_E`$(kpc), $`c`$, $`A_{cf}`$) = (1.0, 35, 1.4, 0.15) and (0.4, 12, 2, 0.14) for Si and O respectively. The $`A_{cf}`$ are chosen for agreement with the XMM-Newton abundances at r
>30
>𝑟30r\mathrel{\mathchoice{\lower 3.0pt\vbox{\halign{$\mathsurround 0pt\displaystyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 3.0pt\vbox{\halign{$\mathsurround 0pt\textstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 3.0pt\vbox{\halign{$\mathsurround 0pt\scriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 3.0pt\vbox{\halign{$\mathsurround 0pt\scriptscriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}}30 kpc; $`r_E`$ is assumed to be similar to $`r_e=8.625`$ kpc appropriate for a stellar contribution; $`c`$ is chosen for a reasonable fit to the XMM-Newton data in 10
<r
<30
<10𝑟
<3010\mathrel{\mathchoice{\lower 2.5pt\vbox{\halign{$\mathsurround 0pt\displaystyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.5pt\vbox{\halign{$\mathsurround 0pt\textstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.5pt\vbox{\halign{$\mathsurround 0pt\scriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.5pt\vbox{\halign{$\mathsurround 0pt\scriptscriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}}r\mathrel{\mathchoice{\lower 2.5pt\vbox{\halign{$\mathsurround 0pt\displaystyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.5pt\vbox{\halign{$\mathsurround 0pt\textstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.5pt\vbox{\halign{$\mathsurround 0pt\scriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.5pt\vbox{\halign{$\mathsurround 0pt\scriptscriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}}30 kpc; $`A_o`$ is adjusted until the total luminosity within $`r=50`$ kpc (where all lines are optically thin) matches that of each line computed with $`A_E^{(ob)}(r)`$. With these assumed intrinsic abundance profiles $`A_E(r)`$, the total line luminosities within 50 kpc exceed the approximate observed luminosities of the Si and O lines by about 2 and 13 percent respectively. We exceed the oxygen line luminosity based on $`A_O^{(ob)}(r)`$ since this abundance is less accurately determined by the XMM-Newton data. The cumulative line center optical depths are somewhat greater when evaluated with $`A_E(r)`$ than with $`A_E^{(ob)}(r)`$ as shown in Figure 1.

Monte Carlo radiative transfer for each line has been calculated with $`10^4`$ photons, 50 spatial zones and abundances $`A_E(r)`$. Results are illustrated in Figure 2. The light dotted histograms show the abundance of each element that would be inferred by observing each line in the optically thin limit, $`A_E(birth)=A_E(r)N_{ph,b}/N`$, where $`N_{ph,b}`$ is the number of photons created in each zone and $`N=200`$ is the mean number of photons per zone. The adjacent light-line histograms show the variation of the apparent abundance indicated by each line after radiative diffusion, $`A_E(escape)=A_E(r)N_{ph,e}/N`$, where $`N_{ph,e}`$ is the number of photons that escape from each zone with no further scattering. $`A_E(escape)`$ can be compared directly with the observed abundances.

The light solid-line histograms clearly indicate central photon depletions or flattening within $`5`$ kpc. The apparent Si abundance based on the 2.003 keV line is clearly self-absorbed and matches the observations quite well at r
>3
>𝑟3r\mathrel{\mathchoice{\lower 3.0pt\vbox{\halign{$\mathsurround 0pt\displaystyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 3.0pt\vbox{\halign{$\mathsurround 0pt\textstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 3.0pt\vbox{\halign{$\mathsurround 0pt\scriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 3.0pt\vbox{\halign{$\mathsurround 0pt\scriptscriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}}3 kpc. But the amount of photon depletion with conservative scattering is clearly insufficient to explain the observed O abundance. It is not possible to deepen the central minima by increasing the O abundance $`A_O(r)`$ (and line opacity) in r
<5
<𝑟5r\mathrel{\mathchoice{\lower 2.5pt\vbox{\halign{$\mathsurround 0pt\displaystyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.5pt\vbox{\halign{$\mathsurround 0pt\textstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.5pt\vbox{\halign{$\mathsurround 0pt\scriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.5pt\vbox{\halign{$\mathsurround 0pt\scriptscriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}}5 kpc since the global line luminosities would exceed those observed. The OVIII line luminosity may already be too high by about 13 percent. The central oxygen line luminosity (abundance) only levels off and does not decrease toward the origin as observed. Evidently line photons are being lost by absorption. The central abundances of Si and O – $`A_E(0)=1.15`$ and 0.54 respectively – cannot be trusted if the line scattering is non-conservative. Since the continuous opacity in a warm plasma with Tw
<106
<subscript𝑇𝑤superscript106T_{w}\mathrel{\mathchoice{\lower 2.5pt\vbox{\halign{$\mathsurround 0pt\displaystyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.5pt\vbox{\halign{$\mathsurround 0pt\textstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.5pt\vbox{\halign{$\mathsurround 0pt\scriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.5pt\vbox{\halign{$\mathsurround 0pt\scriptscriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}}10^{6}K decreases rapidly with photon energy, $`\sigma E^3`$, the OVIII line should suffer $`30`$ times more absorption than the SiXIV line, in qualitative agreement with the profiles in Figure 2.

Shigeyama (1998) noted that if the outward diffusion of strong X-ray lines is not considered in the data reduction, the central gas density can be underestimated by $`2050`$ percent. To explore this, we repeated the calculations with $`n_{e,\beta }`$ increased by 1.5. Our conclusions are strengthened by this adjustment in the gas density. With higher $`n_{e,\beta }`$ the luminosity of the lines within 50 kpc can be matched with a slightly lower central abundance, resulting in central minima in $`A_E(escape)(r)`$ even less pronounced than those in Figure 2.

Studies of optical line emission in the central regions of galactic and galaxy group cooling flows reveal widespread random motions typically about $`0.20.4`$ of the sound speed in the hot gas (Heckman et al. 1989; Caon et al. 2000). Since small, optically visible line-emitting regions at $`T10^4`$ K are likely to be strongly coupled to the ambient hot gas, we infer that the hot gas must share the same turbulent velocities. Turbulence is expected to strongly influence the spatial distribution of X-ray resonance lines, particularly for more massive line-emitting ions (Equation 4). To demonstrate this, we repeated the transfer calculations for both lines with $`v_{turb}^2=_{turb}^2(5A_E/6\mu )`$ assuming a turbulent Mach number $`_{turb}=0.3`$. In the presence of turbulence, which we regard as more realistic, the line center optical depths are greatly reduced (Fig. 1), there is less outward diffusion (Fig. 2), and the apparent central depletion of silicon photons is significantly reduced.

We have also computed line profiles for the FeXXV 6.69 keV (He 4) resonance line discussed by Böhringer et al. but have not discussed its profile here because emission in this line is unlikely to be representative of the intrinsic iron abundance. The FeXXV 6.69 keV line is very sensitive to the (low) cooling flow temperatures $`T(r)`$ in r
<10
<𝑟10r\mathrel{\mathchoice{\lower 2.5pt\vbox{\halign{$\mathsurround 0pt\displaystyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.5pt\vbox{\halign{$\mathsurround 0pt\textstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.5pt\vbox{\halign{$\mathsurround 0pt\scriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.5pt\vbox{\halign{$\mathsurround 0pt\scriptscriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}}10 kpc where the line transfer is most relevant (e.g. Fig. 2 of Buote, Canizares & Fabian 1999) and relatively few FeXXV line photons are created in this low temperature region. Because of this, the $`\chi ^2`$ iron abundance determination in M87 by Böhringer et al. is dominated by the multitude of strong Fe L lines near $`1`$ keV. The many lines that constitute the Fe L feature have different atomic parameters, so no single L line can faithfully represent the iron abundance, as we have assumed for the O and Si Lyman-$`\alpha `$ lines. Since the continuous plasma opacity varies as $`\sigma E^3`$, the opacity at the FeXXV 6.69 keV is $`1000`$ times less than the opacity at the OVIII 0.652 keV line. An observation of the equivalent width of the essentially unabsorbed FeXXV 6.69 keV line with projected radius in M87 would allow a calculation of the intrinsic iron abundance distribution, provided $`T(r)`$ is sufficiently well determined.

We conclude that the central minima in apparent abundances and line luminosities in M87 cannot be understood solely in terms of conservative photon scattering. An additional source of continuous opacity seems to be required. We speculate that the source of this opacity is the warm absorbing gas at $`T10^510^6`$ K as suggested for M87 by Buote (2000). Buote’s conclusion is derived from modeling the ROSAT PSPC data of M87 with coronal plasma emission modified by a single absorption edge and Galactic absorption. The rest energy of the edge is determined to lie at $`0.6`$ keV, near the position of the OVIII L$`\alpha `$ line. The absorption derived from the ROSAT data ($`\tau 0.51`$) is spatially more extended than that required to account for central minima in the X-ray resonance lines. Because of the limited energy resolution of ROSAT, the effect of using a simple absorption edge cannot be distinguished from a model in which the O abundance is allowed to take low values as done by Böhringer et al (2000).

The central declines in Fe and Si abundances in the Böhringer et al. data are difficult to understand since we expect hot gas closer to the galactic core to be more enriched by centrally peaked stellar mass loss and Type Ia supernovae (Brighenti & Mathews 1999b). Since the oxygen abundance found by Böhringer et al. has a generally negative gradient beyond $`r_e`$, it is apparent that the stellar and Type Ia ejecta are more oxygen-rich than cooling flow gas entering M87 from larger radii, so a central peak in the O abundance would also be expected. The need for an additional source of continuous opacity in our line transfer calculations lends new support to the absorption interpretation. If this absorption arises from warm gas, then its mass is similar to that of the hot gas (Buote 2000).

The energy source that heats the warm gas is of great interest. While some of the cooling flow gas may be heated by the M87 jet, this region is avoided in the Böhringer et al. observations; the deep central minimum in their $`T(r)`$ shows no evidence of AGN heating. The most likely energy source for heating the warm gas is the turbulent motion of the hot gas as indicated by the dispersion velocity of optical emission lines $`v_t200`$ km s-1 in M87 (Heckman et al. 1989). These velocities are subsonic in the hot gas, but would be supersonic in the warm, X-ray absorbing gas. Colliding regions of warm gas would be shock-heated at at a rate $`H=fn_cm_pv_t^3/r_c`$ ergs cm-3 s-1 where $`f`$ is the filling factor of warm gas clouds of mean radius $`r_c`$ and electron density $`n_c`$. Although we expect that some net cooling occurs, if shock heating balances radiative cooling $`C=n_c^2\mathrm{\Lambda }`$ ergs cm-3 s-1 then the cloud size can be estimated: $`r_cfm_pv_t^3/n_c\mathrm{\Lambda }400f`$ pc for $`v_t=200`$ km s-1, $`n_c=0.1`$ cm-3 and cooling coefficient $`\mathrm{\Lambda }=10^{22}`$ erg cm3 s-1. With $`f0.5`$, this cloud size is consistent with the scale of irregularities observed in the optical line emission and the relative uniformity of the warm plasma absorption.

Finally, we speculate further that such extended absorbing regions are generally present in cooling flows and account for the apparent deficiency of emission lines characteristic of cooling ions in the centers of these flows (e.g. Peterson et al. 2000; McNamara et al. 2000).

We thank referee Paul Nulsen for his helpful suggestions. Studies of the evolution of hot gas in elliptical galaxies at UC Santa Cruz are supported by NASA grant NAG 5-3060 and NSF grant AST-9802994 for which we are very grateful. DAB was supported by NASA Chandra Fellowship grant PF8-10001 awarded by the Smithsonian Astrophysical Observatory under contract NAS8-39073.