text:
# References

Non trivial generalizations of the Schwinger pair production result
J. AVAN111Supported in part by a CNRS/Brown Accord

LPTHE, Tour 16
Université Pierre et Marie Curie
75252 Paris Cedex 05 France
H. M. FRIED

Department of Physics
Brown University
Providence R.I. 02912 USA
Y. GABELLINI

Institut Non Linéaire de Nice
1361 Route des Lucioles
06560 Valbonne France

Abstract

We present new, non trivial generalizations of the recent Tomaras, Tsamis and Woodard extension of the original Schwinger formula for charged pair production in a constant electric field. That extension generalized the Schwinger result to electric fields $`E_3(x_\pm )`$ dependent upon one $`\underset{¯}{\mathrm{or}}`$ the other light–cone coordinates, $`x_+`$ or $`x_{}`$, $`x_\pm =x_3\pm x_0`$; the present work generalizes their result to electric fields $`E_3(x_+,x_{})`$ dependent upon $`\underset{¯}{\mathrm{both}}`$ coordinates. Displayed in the form of a final, functional integral, or equivalent linkage operation, our result does not appear to be exactly calculable in the general case; and we give a simple, approximate example when $`E_3(x_+,x_{})`$ is a slowly varying function of its variables. We extend this result to the more general case where $`\stackrel{}{E}`$ can point in a varying direction, and where an arbitrary magnetic field $`\stackrel{}{B}`$ is present; both extensions can be cast into the form of Gaussian–weighted functional integrals over well defined factors, which are amenable to approximations depending on the nature and variations of the fields.

A recent, non trivial generalization of Schwinger’s 1951 calculation of the probability/vol.sec for $`e^+e^{}`$ production in a constant electric field has been given by Tomaras, Tsamis and Woodard ( TTW )  in which the electric field may depend upon either light–cone coordinates $`x_\pm =x_3\pm x_0`$, but not upon both. The form of their result is exactly the same as Schwinger’s , which raises the question as to whether further generalizations are possible; and if so, what form they would take.

Shortly after the TTW paper appeared, a second and independent calculation of pair production in an electric field depending on either $`x_+`$ or $`x_{}`$ was performed  using functional techniques, and verifying the result of ref.. Those functional methods are here employed to attempt a further generalization to the case where the electric field depends upon both light–cone variables, $`E=E_3(x_+,x_{})`$. The result is more complicated than the original Schwinger form, although the important ( and non trivial ) essential singularity in the region of small coupling is apparently preserved. The exact statement here requires the evaluation of an additional ( and non trivial ) functional integral; however, for certain situations, such as particle production in the overlap volume of a pair of high intensity lasers  reasonable kinematic approximations are surely justified. If

$$\mathrm{ln}P_0=\frac{\alpha E^2}{\pi ^2}\underset{n=1}{\overset{\mathrm{}}{}}\frac{1}{n^2}e^{n\pi m^2/gE},\alpha =\frac{g^2}{4\pi }$$

denotes the logarithmic density of the exact vacuum–persistence probability density of the Schwinger constant field calculation, the main result of this paper can be expressed as

$$\mathrm{ln}P_0^{}=\frac{\alpha E^2}{\pi ^2}\underset{n=1}{\overset{\mathrm{}}{}}\frac{1}{n^2}e^{n\pi m^2/gE}_n$$

with $`E=E(x_+,x_{})`$, $`\tau _n=n\pi /gE`$, $`J(\lambda _1,\lambda _2)=\theta (\lambda _1\lambda _2)\lambda _2+\theta (\lambda _2\lambda _1)\lambda _1\lambda _1\lambda _2`$, and

$`_n=(1)^ne^{4m^2\tau _n{\displaystyle _0^1𝑑\lambda _1𝑑\lambda _2\frac{\delta }{\delta v_+(\lambda _1)}J(\lambda _1,\lambda _2)\frac{\delta }{\delta v_{}(\lambda _2)}}}`$
$`\times \mathrm{cos}\left[g\tau _n{\displaystyle _0^1}d\lambda E(x_+m^1v_+(\lambda ),x_{}m^1v_{}(\lambda ))\right]|_{v_\pm 0}`$ (1)
where the “ linkage operation ” of $`(1)`$ may be cast into a corresponding functional integral. It doesn’t seem to be possible to evaluate $`(1)`$ exactly, although kinematic approximations are certainly possible. One of these will be used below to illustrate the formula. Note that if $`E`$ depends on only one variable $`x_+`$ or $`x_{}`$, $`_n1`$, and the TTW extension of the Schwinger result is obtained.

We hold to the notation and technique of ref., beginning with the exact statement of the vacuum persistence amplitude in the absence of radiative corrections

$$<0|S|\mathrm{\hspace{0.17em}0}>=e^{L\left[A_{ext}\right]}$$
(2)
The vacuum persistence probability of the present problem, here called $`P_0^{}`$, is then given by $`P_0^{}=\mathrm{exp}\left[2\mathrm{Re}L[A_{ext}]\right]`$, and the probability of producing ( at least ) one charged pair is $`P_1^{}=1P_0^{}`$. Note that $`L[A]`$ is really $`L[F]`$, since $`L`$ is rigorously gauge invariant; but we shall continue to use the vector–potential description, beginning with the exact Fradkin representation

$`L[A]={\displaystyle \frac{1}{2}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{ds}{s}}e^{ism^2}\mathrm{exp}\left[i{\displaystyle _0^s}𝑑s^{}{\displaystyle \underset{\mu }{}}{\displaystyle \frac{\delta ^2}{\delta v_\mu ^2(s^{})}}\right]`$
$`\times {\displaystyle }{\displaystyle \frac{d^4p}{(2\pi )^4}}e^{ip!{\displaystyle _0^s}ds^{}v\left(s^{}\right)}{\displaystyle }d^4xe^{ig{\displaystyle _0^s}𝑑s^{}v_\mu \left(s^{}\right)A_\mu \left(x{\displaystyle _0^s^{}}v\right)}`$ (3)
$`\times \mathrm{tr}\left(e^{g{\displaystyle _0^s}𝑑s^{}\sigma F\left(x{\displaystyle _0^s^{}}v\right)}\right)_+|_{v_{mu}0}(g0)`$
The functional linkage operation of $`(3)`$ may be recast into that of a Gaussian–weighted functional integral, since

$$e^{i{\displaystyle _0^s}𝑑s{\displaystyle \frac{\delta ^2}{\delta v^2(s^{})}}}[v]|_{v0}=Nd[v]e^{{\displaystyle \frac{i}{4}}{\displaystyle _0^s}𝑑sv^2\left(s^{}\right)}[v]$$

for arbitrary $`[v]`$, with the normalization $`N`$ given by

$$N^1=d[v]\mathrm{exp}\left[\frac{i}{4}_0^s𝑑sv^2(s^{})\right]$$

The linkage formulation, which dispenses with normalization constants, is somewhat more convenient, and will be followed whenever possible.

To represent a given, external field in the $`\widehat{3}`$ direction, we may choose $`A_\mu =(\stackrel{}{A}_{},A_3,A_0)(0,A_3,A_0)(0,A_+,A_{})`$ with $`\stackrel{}{A}_{}=0`$, $`A_\pm (x)=_\mu n_\pm ^\mu A_\mu (x_+,x_{})`$. Here, $`n_\pm ^\mu =(0,0,1,1)`$, so that $`n_\pm a=a_\pm =a_3\pm a_0`$. Note that $`n_+^2=n_{}^2=0`$, $`n_+n_{}=2`$, and that $`pAp_3A_3p_0A_0=\frac{1}{2}(p_+A_{}+p_{}A_+)`$. We shall choose the gauge specified by $`A_{}=0`$, so that $`A^2=A_{}^2+A_+A_{}=0`$, and $`E(x_+,x_{})={\displaystyle \frac{A_+}{x_+}}(x_+,x_{})`$.

As in ref., we extract the $`v_\mu (s^{\prime \prime })`$ dependence inside the arguments of $`A_\mu `$ and $`F_{\mu \nu }`$, or of $`A_+`$ and $`E`$, by introducing a functional form of unity into $`(3)`$, replacing $`[{\displaystyle _0^s^{}}v_+(s^{\prime \prime }),{\displaystyle _0^s^{}}v_{}(s^{\prime \prime })]`$ by

$$d[u_+]d[u_{}][u_+(s^{}),u_{}(s^{})]\delta \left[u_+(s^{})_0^s^{}v_+(s^{\prime \prime })\right]\delta \left[u_{}(s^{})_0^s^{}v_{}(s^{\prime \prime })\right]$$
(4)
where the $`\delta `$–functional notation means that when the region $`0`$$`s`$ is broken up into many small intervals labeled by the discrete indices $`s_i`$, $`i=1,\mathrm{},N`$ ( with $`N\mathrm{}`$, subsequently ), $`\delta \left[Q(s^{})\right]{\displaystyle \underset{i=1}{\overset{N}{}}}\delta \left(Q(s_i)\right)`$.

Note that since we are treating each $`v_\pm (s^{})`$ as a continuous ( although not necessarily differentiable ) function of $`s^{}`$, the $`u_\pm (s^{})`$ introduced in $`(4)`$ have, at least, a continuous first derivative.

Each $`\delta `$ function of $`(4)`$, and the overall $`\delta `$ functional, may be written in terms of a standard Fourier representation; so that, in the continuous limit, $`(4)`$ may be replaced by

$`N^2{\displaystyle d[u_+]d[u_{}][u_+(s^{}),u_{}(s^{})]d[\mathrm{\Omega }_+]d[\mathrm{\Omega }_{}]}`$ (5)
$`\times e^{i{\displaystyle _0^s}ds^{}[u_+\left(s^{}\right)\mathrm{\Omega }_+\left(s^{}\right)+u_{}\left(s^{}\right)\mathrm{\Omega }_{}\left(s^{}\right]}e^{i{\displaystyle _0^s}𝑑s^{}\left[\mathrm{\Omega }_+\left(s^{}\right){\displaystyle _0^s^{}}v_+\left(s^{\prime \prime }\right)+\mathrm{\Omega }_{}\left(s^{}\right){\displaystyle _0^s^{}}v_{}\left(s^{\prime \prime }\right)\right]}`$
where $`N^{}`$ is an appropriate normalization factor ( which disappears from the final result ). With the aid of Abel’s trick, the last exponential factor of $`(5)`$ may be written as

$$\mathrm{exp}\left\{i_0^s𝑑s^{}v_\mu (s^{})\left[n_+^\mu _s^{}^s𝑑s^{\prime \prime }\mathrm{\Omega }_+(s^{\prime \prime })+n_{}^\mu _s^{}^s𝑑s^{\prime \prime }\mathrm{\Omega }_{}(s^{\prime \prime })\right]\right\}$$

so that the entire $`v`$–linkage operation is immediate :

$`e^{i{\displaystyle _0^s}𝑑s{\displaystyle \frac{\delta ^2}{\delta v^2(s^{})}}}e^{i{\displaystyle _0^s}𝑑s^{}v_\mu \left(s^{}\right)\left[p_\mu gA_\mu \left(s^{}\right)n_+^\mu {\displaystyle _s^{}^s}𝑑s^{\prime \prime }\mathrm{\Omega }_+\left(s^{\prime \prime }\right)n_{}^\mu {\displaystyle _s^{}^s}𝑑s^{\prime \prime }\mathrm{\Omega }_{}\left(s^{\prime \prime }\right)\right]}|_{v0}`$
$`=\mathrm{exp}\{i{\displaystyle _0^s}ds^{}[p_\mu gA_\mu (s^{})n_+^\mu {\displaystyle _s^{}^s}\mathrm{\Omega }_+n_{}^\mu {\displaystyle _s^{}^s}\mathrm{\Omega }_{}]^2\}`$ (6)
where we have used the notation $`A_\mu (s^{})A_\mu (x_+u_+(s^{}),x_{}u_{}(s^{}))`$. As in ref., the explicitly quadratic $`\mathrm{\Omega }_+`$ and $`\mathrm{\Omega }_{}`$ terms of the exponential factor of $`(6)`$ are removed because $`n_+^2=n_{}^2=0`$, but there remains a non zero $`\mathrm{\Omega }_+\mathrm{\Omega }_{}`$ cross term, so that $`(6)`$ becomes

$`\mathrm{exp}\{isp^2+igp_{}{\displaystyle _0^s}ds^{}A_+(s^{})+2ip_+{\displaystyle _0^s}ds^{}s^{}\mathrm{\Omega }_+(s^{})+2ip_{}{\displaystyle _0^s}ds^{}s^{}\mathrm{\Omega }_{}(s^{})`$
$`2ig{\displaystyle _0^s}ds^{}\mathrm{\Omega }_+(s^{}){\displaystyle _0^s^{}}ds^{\prime \prime }A_+(s^{\prime \prime })4i{\displaystyle _0^s}ds_1{\displaystyle _0^s}ds_2\mathrm{\Omega }_+(s_1)h(s_1,s_2)\mathrm{\Omega }_{}(s_2)\}`$ (7)
where $`h(s_1,s_2)=\theta (s_1s_2)s_2+\theta (s_2s_1)s_1=\frac{1}{2}(s_1+s_2|s_1s_2|)`$ is that function of proper time which always appears when constructing representations of $`G_c[A]`$ . In writing $`(7)`$, we have used the “ inverse ” form of Abel’s manipulation.

Performing the $`d^4p`$, one requires

$`{\displaystyle d^4pe^{isp^2+ip\left[n_{}g{\displaystyle _0^s}A\left(s^{}\right)+2n_+{\displaystyle _0^s}s^{}\mathrm{\Omega }_+\left(s^{}\right)+2n_{}{\displaystyle _0^s}s^{}\mathrm{\Omega }_{}\left(s^{}\right)\right]}}`$
$`=i\left({\displaystyle \frac{\pi ^2}{s^2}}\right)\mathrm{exp}\left\{4i{\displaystyle _0^s𝑑s_1𝑑s_2\mathrm{\Omega }_+(s_1)\left(\frac{s_1s_2}{s}\right)\mathrm{\Omega }_{}(s_2)}+{\displaystyle \frac{2i}{s}}{\displaystyle _0^s}s^{}\mathrm{\Omega }_+(s^{})g{\displaystyle _0^s}A_+(s^{})\right\}`$ (8)
where we have again used the properties : $`n_+^2=n_{}^2=0`$, $`n_+n_{}=2`$. The new $`\mathrm{\Omega }_+\mathrm{\Omega }_{}`$ term of $`(8)`$ may be combined with the previous, like term of $`(7)`$, forming the combination $`\mathrm{exp}\left\{4i{\displaystyle _0^s\mathrm{\Omega }_+J\mathrm{\Omega }_{}}\right\}`$, $`J(s_1,s_2)=h(s_1,s_2){\displaystyle \frac{s_1s_2}{s}}`$; and the functional integral over $`\mathrm{\Omega }_\pm `$ may be written as

$`e^{4i{\displaystyle \frac{\delta }{\delta u_+(s_1)}J(s_1,s_2)\frac{\delta }{\delta u_{}(s_2)}}}N^2{\displaystyle d[\mathrm{\Omega }_+]d[\mathrm{\Omega }_{}]e^{i{\displaystyle _0^s}𝑑s^{}u_{}\left(s^{}\right)\mathrm{\Omega }_{}\left(s^{}\right)}}`$
$`\times e^{i{\displaystyle _0^s}𝑑s^{}\mathrm{\Omega }_+\left(s^{}\right)\left[u_+\left(s^{}\right)+2g{\displaystyle \frac{s^{}}{s}}{\displaystyle _0^s}𝑑s^{\prime \prime }A_+\left(s^{\prime \prime }\right)2g{\displaystyle _0^s^{}}𝑑s^{\prime \prime }A_+\left(s^{\prime \prime }\right)\right]}`$ (9)
$`=e^{4i{\displaystyle \frac{\delta }{\delta u_+}J\frac{\delta }{\delta u_{}}}}\delta \left[u_{}(s^{})\right]\delta \left[u_+(s^{})+2g{\displaystyle \frac{s^{}}{s}}{\displaystyle _0^s}A_+(s^{\prime \prime })2g{\displaystyle _0^s^{}}A_+(s^{\prime \prime })\right]`$
Eq.$`(9)`$ is multiplied by the $`u_\pm `$–dependent factor

$$\mathrm{tr}\left(e^{g{\displaystyle _0^s}𝑑s^{}\sigma F(x_+u_+\left(s^{}\right),x_{}u_{}\left(s^{}\right))}\right)_+=4\mathrm{cosh}\left(g_0^s𝑑s^{}E(s^{})\right)$$
(10)
where $`E(s^{})E(x_+u_+(s^{}),x_{}u_{}(s^{}))`$. The $`{\displaystyle d[u_+]d[u_{}]}`$ may be evaluated by transferring ( by an integration by parts over each $`u`$ variable ) the differential $`\delta /\delta u_\pm `$ operators which act upon the $`\delta `$ functionals of $`(9)`$ to equivalent operations upon the $`u_\pm `$ dependence of $`(10)`$

$$4|det\frac{\delta u}{\delta f}|e^{4i{\displaystyle \frac{\delta }{\delta u_+}J\frac{\delta }{\delta u_{}}}}\mathrm{cosh}\left(g_0^sds^{}E(s^{})\right)|_{u_{}(s^{})=0,u_+(s^{})=u(s^{})}$$

where $`|det{\displaystyle \frac{\delta u}{\delta f}}|`$ is the determinant of the transformation from the variables $`u_+(s^{})`$ to the variables $`f(s^{})`$, the latter given by the argument of the $`\delta `$ functional produced by the $`{\displaystyle d[\mathrm{\Omega }_+]}`$. Here, $`u(s^{})`$ is the solution of the integral equation

$$u(s^{})=2g\left[_0^s^{}𝑑s^{\prime \prime }A_+(s^{\prime \prime })\frac{s^{}}{s}_0^s𝑑s^{\prime \prime }A_+(s^{\prime \prime })\right]$$
(11)
with $`u_{}(s^{})=0`$, corresponding to the other $`\delta `$ functional of $`(9)`$.

By exactly the same arguments as in , the only possible solution to $`(11)`$ is $`u(s^{})=u^{}(s^{})=0`$, which means that, as in ,

$$|det\frac{\delta u}{\delta f}|=\frac{gsE(x_+,x_{})}{\mathrm{sinh}(gsE(x_+,x_{}))}$$
(12)
with the difference that this $`E`$ can depend upon both $`x_+`$ and $`x_{}`$.

The sign of $`E`$ entering into $`(10)`$ and $`(12)`$ is irrelevant, and we shall suppose it is positive. Our expression for $`L[A]`$ can therefore be put into a form which closely resembles that of Schwinger, and of ref..

$`L[A]={\displaystyle \frac{i}{8\pi ^2}}{\displaystyle d^4x_0^{\mathrm{}}\frac{ds}{s^2}e^{ism^2}\frac{gE}{\mathrm{sinh}(gsE)}e^{4ism^2{\displaystyle _0^1𝑑\lambda _1𝑑\lambda _2\frac{\delta }{\delta v_+(\lambda _1)}J(\lambda _1,\lambda _2)\frac{\delta }{\delta v_{}(\lambda _2)}}}}`$
$`\times \mathrm{cosh}\left[gs{\displaystyle _0^1}𝑑\lambda E(x_+m^1v_+(\lambda ),x_{}m^1v_{}(\lambda ))\right]|_{v_\pm 0}(g0)`$ (13)
In writing $`(13)`$, we have observed that, in the interval $`0`$$`s`$, any continuous function $`u_\pm (s^{})`$ may be given by a Fourier series representation as a function of $`(s^{}/s)`$, and we therefore replace $`u_\pm (s^{})`$ by $`m^1v_\pm (s^{}/s)m^1v_\pm (\lambda )`$, where $`\lambda =s^{}/s`$, and the factor $`m^1`$ is inserted for dimensional reasons ( and will cancel away when the functional derivatives $`\delta /\delta v_\pm `$ are taken ). However, because

$$\frac{\delta u_\pm (s_1)}{\delta u_\pm (s_2)}=\delta (s_1s_2)=\frac{1}{s}\delta (\lambda _1\lambda _2),\mathrm{and}\frac{\delta v_\pm (\lambda _1)}{\delta v_\pm (\lambda _2)}=\delta (\lambda _1\lambda _2)$$

one must adopt the relation $`{\displaystyle \frac{\delta }{\delta u_\pm (s^{})}}={\displaystyle \frac{m}{s}}{\displaystyle \frac{\delta }{\delta v_\pm (\lambda ^{})}}`$. After extracting an overall factor of $`s`$, $`J(\lambda _1,\lambda _2)`$ is understood to be given by

$$J(\lambda _1,\lambda _2)=\theta (\lambda _1\lambda _2)\lambda _2+\theta (\lambda _2\lambda _1)\lambda _1\lambda _1\lambda _2$$

Extracting the $`\mathrm{Re}L[A]`$ from an expression involving an operator such as that of $`(13)`$ requires some further thought, and we therefore imagine that the linkage operator of $`(13)`$ is expanded in powers of $`J`$, so that the rotation of contour $`si\tau `$ is permitted; here, that expansion yields powers of $`4m^2\tau {\displaystyle _0^1\frac{\delta }{\delta v_+(\lambda _1)}J(\lambda _1,\lambda _2)\frac{\delta }{\delta v_{}(\lambda _2)}}`$, which generate real results because $`E`$ and all of its derivative are real. Under this contour rotation, one again finds that the contributions to $`\mathrm{Re}L[A]`$ arise from zeroes of the denominator $`\mathrm{sin}(g\tau E)`$, occuring for $`\tau \tau _ni\epsilon `$, $`\tau _n=n\pi /gE`$; and we then sum up all $`J`$–dependent terms, so that

$$2\mathrm{R}\mathrm{e}L[A]=\frac{\alpha }{\pi ^2}d^4xE^2\underset{n=1}{\overset{\mathrm{}}{}}\frac{1}{n^2}e^{n\pi m^2/gE}_n(x_+,x_{})$$

with $`_n`$ defined in $`(1)`$. This seems to be as far as the exact analysis can be performed.

Without approximation, the linkage operation of $`(13)`$ may be converted to a pair of functional integrals,

$`e^{iab{\displaystyle \frac{\delta }{\delta v_+}J\frac{\delta }{\delta v_{}}}}[v_+,v_{}]|_{v_\pm 0}`$
$`=det[J]{\displaystyle d[\chi _+]d[\chi _{}]e^{i{\displaystyle \chi _+J\chi _{}}}[aJ\chi _+,bJ\chi _{}]}`$ (14)
and mean field/stationary phase approximations derived for the right hand side of $`(14)`$. Perhaps the simplest approximation of all occurs when $`E(x_+,x_{})=E_+(x_+)+E_{}(x_{})`$, a form suggested by recent estimates of pair production in the overlap region of two, crossed, high intensity lasers . If all the derivatives of $`E`$ of order higher than the second are neglected,

$`{\displaystyle _0^1}𝑑\lambda E_\pm (x_\pm m^1v_\pm (\lambda ))E_\pm (x_\pm )m^1{\displaystyle _0^1}𝑑\lambda v_\pm (\lambda ){\displaystyle \frac{E_\pm (x_\pm )}{x_\pm }}`$
$`+{\displaystyle \frac{m^2}{2}}{\displaystyle _0^1}𝑑\lambda v_\pm ^2(\lambda ){\displaystyle \frac{^2E_\pm (x_\pm )}{x_\pm ^2}}+\mathrm{}`$
which approximation assumes that the fractional variations of the electric field are very small over distances on the order of the Compton wavelenght.

Then, it is easy to see that the linkage operation yields

$$_n=\mathrm{cos}\left[\frac{a^2}{2}Q(\alpha _+^2\gamma _{}+\alpha _{}^2\gamma _+)\right]\mathrm{exp}\left[4a^2\alpha _+\alpha _{}R\frac{1}{2}\mathrm{Tr}\mathrm{ln}(1+a^2\gamma _+\gamma _{}I)\right]$$
(15)
where : $`a=4\pi n/gE`$, $`\alpha _\pm =n\pi {\displaystyle \frac{}{x_\pm }}\mathrm{ln}\left(E/m^2\right)`$, $`\gamma _\pm ={\displaystyle \frac{n\pi }{E}}{\displaystyle \frac{^2}{x_\pm ^2}}E`$, $`E=E_++E_{}`$ and :

$$Q=_0^1𝑑\lambda _1𝑑\lambda _2<\lambda _1|I\left(1+a^2\gamma _+\gamma _{}I\right)^1|\lambda _2>$$

$$R=_0^1𝑑\lambda _1𝑑\lambda _2<\lambda _1|J\left(1+a^2\gamma _+\gamma _{}I\right)^1|\lambda _2>$$

with

$$I(\lambda _1,\lambda _2)=<\lambda _1|I|\lambda _2>=_0^1𝑑\lambda J(\lambda _1,\lambda )J(\lambda ,\lambda _2)$$

Note that, formally, for large $`n`$, $`Rn^4`$ while $`a^2x_+x_{}n^4`$, so that the exponential terms of $`(15)`$ cannot grow rapidly with $`n`$. One therefore infers that the essential singularity structure of the original Schwinger result is preserved.

For an electric field $`E_3`$, which depends upon $`x_{}`$, one can retain the condition $`A_{}=A_{}=0`$, but require $`A_+`$ to depend on $`x_{}`$, as well as $`x_+`$ and $`x_{}`$. One then introduces variables $`u_{}(s^{})`$ and $`\mathrm{\Omega }_{}(s^{})`$, and following the above analysis finds

$`L[A]={\displaystyle \frac{i}{32\pi ^2}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{ds}{s^3}}e^{ism^2}{\displaystyle d^4xd[u_+]d[u_{}]d[u_{}]\mathrm{tr}\left(e^{g{\displaystyle _0^s}\sigma F}\right)_+}`$
$`\times e^{4i{\displaystyle \frac{\delta }{\delta u_+}J\frac{\delta }{\delta u_{}}}+i{\displaystyle \frac{\delta }{\delta u_{}}J\frac{\delta }{\delta u_{}}}}\delta \left[u_+(s^{})+2g{\displaystyle \frac{s^{}}{s}}{\displaystyle _0^s}A_+(s^{\prime \prime })2g{\displaystyle _0^s^{}}A_+(s^{\prime \prime })\right]`$ (16)
$`\times \delta \left[u_{}(s^{})\right]\delta \left[u_{}(s^{})\right]`$
An integration by part on the $`u_\pm `$, $`u_{}`$ dependence then allows the linkage operator to act upon the $`\mathrm{tr}\left(e^{g{\displaystyle _0^s}\sigma F}\right)_+`$ term only; and one sees that the jacobian of the transformation between the $`u_\pm (s^{})`$ and $`f(s^{})`$ variables is exactly the same as in $`(12)`$, except that here $`E=E_3(x_+,x_{},x_{})`$. If we restrict $`\stackrel{}{E}`$ to lie in the $`3`$ direction, then again, $`\mathrm{tr}\left(e^{g{\displaystyle _0^s}\sigma F}\right)_+4\mathrm{cosh}(gsE)`$ and the limits $`u_\pm =u_{}=0`$ are to be taken after the linkage operator of $`(16)`$ acts upon the $`u_\pm `$, $`u_{}`$ dependence inside $`4\mathrm{cosh}(gsE)`$, with $`E=E_3(x_+u_+(s^{}),x_{}u_{}(s^{}),x_{}u_{}(s^{}))`$.

For arbitrary $`\stackrel{}{E}`$ and $`\stackrel{}{B}`$ fields, one must include at the very beginning $`A_{}(x_+,x_{},x_{})0`$, although it is still convenient to retain the gauge condition $`A_{}=0`$. Now there will appear in this final generalization of $`(16)`$ the extra factor

$$[u]=\mathrm{exp}\left[ig^2_0^s𝑑s^{\prime \prime }A_{}^2(s^{\prime \prime })+i\frac{g^2}{s}\left(_0^s𝑑s^{\prime \prime }A_{}(s^{\prime \prime })\right)^2\right]$$

where, following the notation of $`(6)`$, $`A_{}(s^{})=A_{}(x_+u_+(s^{}),x_{}u_{}(s^{}),x_{}u_{}(s^{}))`$, and $`[u]`$ multiplies the trace of the now more complicated ordered exponential ( OE ) $`\left(e^{g{\displaystyle _0^s}\sigma F}\right)_+`$. Following techniques developed for unitary OEs , the present OE can be approximately evaluated analytically when $`\stackrel{}{E}`$ and $`\stackrel{}{B}`$ directions change in ways that can be characterized as slowly varying ( “ adiabatic ” ) or rapidly fluctuating ( “ stochastic ” ); and the same linkage operator as in $`(16)`$ is then to act upon the product $`[u]\left(e^{g{\displaystyle _0^s}\sigma F}\right)_+`$. Note that there will appear another jacobian, involved in the variable change from $`u_{}(s^{})`$ to $`f_{}(s^{})`$, as in $`(12)`$, and both jacobians will simultaneously involve $`\stackrel{}{E}`$ and $`\stackrel{}{B}`$. Alternatively, if one wishes to avoid the direct evaluation of these jacobians, one can rewrite this generalization of $`(16)`$ in terms of functional integrals, and attempt to use mean field/stationary phase methods for their approximate evaluation. This extension of $`(16)`$ provides a functional representation of $`L[A]`$ appropriate for the most general choice of $`\stackrel{}{E}`$ and $`\stackrel{}{B}`$ fields. When the latter are constants, and in particular for $`\stackrel{}{B}=0`$, it reduces immediately to Schwinger’s 1951 solution, while it generates a definition of $`L[A]`$ at any stage inbetween.

Finally, a more symmetric formulation, explicitly depending only on the $`F_{\mu \nu }`$ may be obtained by starting from formula (5.22) of the book of ref.

$`L[F]={\displaystyle \frac{1}{2}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{ds}{s}}e^{ism^2}{\displaystyle d^4x\frac{d^4p}{(2\pi )^4}e^{i{\displaystyle _0^s}{\displaystyle \frac{\delta ^2}{\delta v^2}}}e^{ip{\displaystyle _0^s}v}}`$
$`\times N^{}{\displaystyle d[\mathrm{\Omega }]d[u]e^{i{\displaystyle _0^s}u_\mu \mathrm{\Omega }_\mu }e^{i{\displaystyle _0^s}\mathrm{\Omega }_\mu {\displaystyle _0^s^{}}v_\mu }}`$ (17)
$`\times e^{ig{\displaystyle _0^s}𝑑s^{}v_\mu \left(s^{}\right)u_\nu (s^{}){\displaystyle _0^1}\lambda 𝑑\lambda F_{\mu \nu }\left(x\lambda u\left(s^{}\right)\right)}\mathrm{tr}\left(e^{g{\displaystyle _0^s}\sigma F\left(xu\left(s^{}\right)\right)}\right)_+(g0)`$
where every factor of $`{\displaystyle _0^s^{}}𝑑s^{\prime \prime }v_\mu (s^{\prime \prime })`$ has been replaced by $`u_\mu (s^{})`$. Performing the linkage operation, $`{\displaystyle d^4p}`$ and $`{\displaystyle d[\mathrm{\Omega }]}`$, then leads to the result

$`L[F]={\displaystyle \frac{i}{32\pi ^2}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{ds}{s^3}}e^{ism^2}{\displaystyle d^4x\underset{\mu =1}{\overset{4}{}}d[u]\mathrm{\Delta }[u_\mu ]}`$
$`\times e^{i{\displaystyle \frac{\delta }{\delta u_\mu }J\frac{\delta }{\delta u_\mu }}}\mathrm{tr}\left(e^{g{\displaystyle _0^s}\sigma F}\right)_+e^Q,`$ (18)
$`\mathrm{\Delta }[u_\mu ]=\delta \left[u_\mu (s^{})2g{\displaystyle _0^s}𝑑s^{\prime \prime }(\theta (s^{}s){\displaystyle \frac{s^{}}{s}}){\displaystyle _0^1}\lambda 𝑑\lambda F_{\mu \nu }(x\lambda u(s^{\prime \prime }))u_\nu (s^{\prime \prime })\right]`$
$`Q=ig^2{\displaystyle _0^s}𝑑s^{}{\displaystyle _0^1}\lambda 𝑑\lambda {\displaystyle _0^1}\lambda ^{}𝑑\lambda ^{}{\displaystyle _0^s}𝑑s^{\prime \prime }u_\nu (s^{})u_\sigma (s^{\prime \prime })F_{\mu \nu }(x\lambda u(s^{}))F_{\mu \sigma }(x\lambda ^{}u(s^{\prime \prime }))\left[\delta (s^{}s^{\prime \prime }){\displaystyle \frac{1}{s}}\right]`$
The only solutions for each $`u_\mu (s^{})`$ allowed by the product of the four $`\delta `$ functional of $`(18)`$ are $`u_\mu (s^{})0`$, although the jacobians obtained from each variable change $`u_\mu (s^{})f_\mu (s^{})`$ will each be unity only for the Schwinger case of constant $`F_{\mu \nu }`$; in this latter situation, the $`\mathrm{exp}\left[i{\displaystyle }{\displaystyle \frac{\delta }{\delta u_\mu }}J{\displaystyle \frac{\delta }{\delta u_\mu }}\right]`$ operation may be recast into that of soluble gaussian functional integrals over the $`u_\mu (s^{})`$, leading back to Schwinger’s original result.


format: text

text:
# New Way to Compute Excited States and Thermodynamics: Monte Carlo Hamiltonian

## 1 MONTE CARLO HAMILTONIAN

The idea of the Monte Carlo Hamiltonian has been introduced in Refs.. Here we want to recall the basic idea, outline some of its characteristic features and present some selected results. Comparing the Lagrangian approach with the Hamiltonian approach in lattice field theory, the Hamiltonian approach has not flourished in producing a wealth of numerical data like the Lagrangian approach. On the other hand, there are a few topics, where the Lagrangian approach has only given slow progress, e.g. on excited states, wave functions, particle scattering. In our opinion, the lattice Hamiltonian approach might have done a lot better, if it were not for the lack of effective many-body methods to solve such Hamiltonian and compute physical observables. Recall that the Kogut-Susskind Hamiltonian of lattice $`QED_{3+1}`$ and $`QCD_{3+1}`$ is a quantum mechanical operator in a high-dimensional Hilbert space and a numerical many-body solution is quite difficult.

The Monte Carlo Hamiltonian using a stochastic basis is an attempt to address this problem. Let us outline the idea. First, the Monte Carlo Hamiltonian is not an operator of canonical form. It is rather defined in a subspace of Hilbert space and is given in terms of matrix elements. The physical idea is that this object plays the role of an effective Hamiltonian. It has a window of validity. This can be visualized as a window in the energy spectrum. Beyond this window, the effective Hamiltonian is unphysical. When computing thermodynamical functions, the window of the energy spectrum translates into a temperature window. Such a window can be observed in the numerical data (see below). The notion of a window is familar in physics from the concept a scaling window and also reminds of Wilson-Kadanoff’s renormalisation group, where a Hamiltonian is constructed which gives the correct physics of critical phenomena (at phase transitions), but is unphysical otherwise.

What is the guiding principle in constructing such effective Hamiltonian? In contrast to the transfer matrix formalism which defines a Hamilton operator in the limit $`a_t0`$, we start by considering physical amplitudes corresponding to finite transition time $`T`$. The finiteness of $`T`$ is in part responsable for the finite energy window. However, in contrast to Lagrangian lattice field theory, where one considers n-point functions corresponding to vacuum-to-vacuum transitions, we go beyond the vacuum state and consider transition amplitudes from other states which probe in part the wave functions of excited states. To give an example in quantum mechanics, we consider transition amplitudes $`G(x_f,t=T;x_i,t=0)`$, where $`|x_i>`$ and $`|x_f>`$ denote eigen states of initial and final position. We compute those quantities via standard path integrals on the lattice using Monte Carlo with importance sampling (the trick is to write the path integral as a ratio). In order to get a reasonable representative subspace of Hilbert space (which also influences the size of the window), we took in quantum mechanics $`x_i`$ and $`x_f`$ from a finite grid of position states and computed all possible combinations of transitions. Then an algebraic diagonalisation of such a matrix yields a finite set of eigenstates and eigen energies of the effective Hamiltonian, i.e. the Monte Carlo Hamiltonian. For test cases from quantum mechanics this has been shown to work .

## 2 MANY-BODY SYSTEMS AND STOCHASTIC BASIS

It is simple to convince oneself that a straight forward generalisation to many-body systems and field theory is not feasible. The solution is to select physically important degrees of freedom. Here we use as guiding principle the enormous success of Lagrangian lattice field theory, and in particular the successful working of Monte Carlo with importance sampling. In Lagrangian lattice field theory it is possible to estimate an observable by taking the mean of the observable over a relatively small number of ”equilibrium’ field configurations, which have been constructed using a weightfactor built from the action.

We suggest that ’equilibrium’ field configurations in Lagrangian lattice field theory have a close analogy in Hamiltonian lattice field theory in terms of states of a stochastic basis. In other words, taking ’equilibrium’ configurations at a fixed time slice yields Hilbert states (Bargman states) which represent physically important degrees of freedom. This constitutes the so-called stochastic basis. Thus, via Monte Carlo with importance sampling, involving the weightfactor with the action, we generate a subspace of Hilbert space of a relatively small dimension (in the order of 100 to 1000), which represents the important degrees of freedom for the physics in a certain window.

The demonstration that such repesentative basis exists and allows to compute excited states and thermodynamical functions has been first made for the chain of coupled harmonic oscillators (Klein-Gordon model). In Ref. the excitation spectrum has been compared with the analytic answer, establishing the existence of an energy window and thermodynamic functions (specific heat) have been computed.

The method has been extended to the scalar model in Ref.. This model is not analytically soluble, hence we have chosen to compare our results with standard Langian lattice calculations. Fig.1 presents the the first 50 eigenvalues of the excitation spectrum of the Monte Carlo Hamiltonian. A sudden increase of errors indicates the upper limit of the energy window. The thermodynamical functions free energy, average energy, entropy and specific heat are displayed in Figs.2-5. The temperature window reaches from temperature zero to an upper temperature given by $`\beta 11.5`$. The size of the window depends on transition time $`T`$, the size of the stochastic basis, statistical errors, etc. In Figs.2 \- 5 one observes that the data from the Monte Carlo Hamiltonian show lesser fluctuations than those from the standard Lagrangian approach.


text:
# NLTE Radiative Transfer and Model Atmospheres of Hot Stars

## 1. Introduction

Stellar atmospheres are open systems and thus cannot be in thermodynamic equilibrium (TE). The “Local Thermodynamic Equilibrium” (LTE) is a working hypothesis which assumes TE not for the atmosphere as a whole but for small volume elements. As a consequence, the atomic population numbers are depending only on the local (electron) temperature and electron density via the Saha-Boltzmann equations. Computing models by replacing the Saha-Boltzmann equations by the rate equations (statistical equilibrium) are called non-LTE (or NLTE) models. NLTE calculations are more costly than LTE calculations, however, it is hard to predict if NLTE effects are important in a specific problem. Generally, NLTE effects are large at high temperatures and low densities, which implies intense radiation fields hence frequent radiative processes and less frequent particle collisions which tend to enforce LTE conditions.

We will restrict ourselves here to the classical model atmosphere problem, i.e. the solution of the radiation transfer equations assuming hydrostatic, radiative and statistical equilibrium. The numerical problem going from LTE to realistic NLTE models has been solved only recently and is the topic of this paper. This was achieved by the development of new numerical techniques for model construction and on the availability of atomic data for many species. The replacement of the Saha-Boltzmann equations by the atomic rate equations requires a different numerical solution technique, otherwise metal opacities cannot be accounted for at all. Such techniques were developed with big success during the last decade, triggered by important papers by Cannon (1973) and Scharmer (1981). The Accelerated Lambda Iteration (ALI) is at the heart of of this development. Combined with statistical methods we are finally able to compute metal line blanketed NLTE models with a very high level of sophistication.

## 2. Overview of problem and solution method

We assume plane-parallel geometry, which is well justified for most stars because the atmospheres are thin compared to the stellar radius. The only parameters which characterize uniquely such an atmosphere are the effective temperature ($`T_{\mathrm{eff}}`$), which is a measure for the amount of energy transported through the atmosphere per unit area and time, the surface gravity ($`g`$), and the chemical composition. To construct model atmospheres we have to solve simultaneously a set of equations that is highly coupled and non-linear. Because of the coupling, no equation is determining uniquely a single quantity – all equations determine a number of state parameters. However, each of them is usually thought of as determining a particular quantity. These equations are:

* The radiation transfer equations which are solved for the (angular) mean intensities $`J_i,i=1,\mathrm{},NF`$, on a pre-chosen frequency grid comprising $`NF`$ points. The formal solution is given by $`J=\mathrm{\Lambda }S`$, where $`S`$ is the source function (Eq. 10). Although $`\mathrm{\Lambda }`$ is written as an operator, one may think of $`\mathrm{\Lambda }`$ as a process of obtaining the mean intensity from the source function.
* The hydrostatic equation which determines the total particle density $`N`$.
* The radiative equilibrium equation from which the temperature $`T`$ follows.
* The particle conservation equation, determining the electron density $`n_e`$.
* The statistical equilibrium equations which are solved for the population densities $`n_i,i=1,\mathrm{},NL`$ of the atomic NLTE levels.
* The definition equation for a fictitious massive particle density $`n_H`$ which is introduced for a convenient representation of the solution procedure.

This set of equations has to be solved at each point $`d`$ of a grid comprising $`ND`$ depth points. Thus we are looking for solution vectors

$$𝝍_d^{}=(n_1,\mathrm{},n_{NL},n_e,T,n_H,N,J_1,\mathrm{},J_{NF}),d=1,\mathrm{},ND.$$
(1)
The Complete Linearization (CL) method (Auer & Mihalas 1969) solves this set by linearizing the equations with respect to all variables. The basic advantage of the ALI (or “operator splitting”) method is that it allows to eliminate at the outset the explicit occurrence of the mean intensities $`J_i`$ from the solution scheme by expressing these variables by the current, yet to be determined, occupation densities and temperature. This is accomplished by an iteration procedure which may be written as (suppressing indices indicating depth and frequency dependency of variables):

$$J^n=\mathrm{\Lambda }^{}S^n+(\mathrm{\Lambda }\mathrm{\Lambda }^{})S^{n1}.$$
(2)
This means that the actual mean intensity at any iteration step $`n`$ is computed by applying an approximate lambda operator (ALO) $`\mathrm{\Lambda }^{}`$ on the actual (thermal) source function $`S^n`$ plus a correction term that is computed from quantities known from the previous iteration step. This correction term includes the exact lambda operator $`\mathrm{\Lambda }`$ which guarantees the exact solution of the radiation transfer problem in the limit of convergence: $`J=\mathrm{\Lambda }S`$. The use of $`\mathrm{\Lambda }`$ in Eq. 2 only indicates that a formal solution of the transfer equation is performed but in our application the operator is not constructed explicitly. Instead we employ the Feautrier solution scheme (Mihalas 1978) or a short characteristic method (Olson & Kunasz 1987) to solve the transfer equation that is set up as a differential equation. The resulting set of equations for the reduced solution vectors

$$𝝍_d=(n_1,\mathrm{},n_{NL},n_e,T,n_H,N),d=1,\mathrm{},ND$$
(3)
is of course still non-linear. The solution is obtained by linearization and iteration which is performed either with a usual Newton-Raphson iteration or other methods (Sect. 5.2.). First model atmosphere calculations with the ALI method were performed by Werner (1986).

Another advantage of the ALI method is that explicit depth coupling of the solution vectors Eq. 1 through the transfer equation can be avoided if one restricts to diagonal (i.e. local) ALOs. Then the solution vectors Eq. 3 are independent from each other and the solution procedure within one iteration step of Eq. 2 is much more straightforward. Depth coupling is provided by the correction term that involves the exact solution of the transfer equation. The hydrostatic equation which also gives an explicit depth coupling, may be taken out of the set of equations and can – as experience shows – be solved in between two iteration steps of Eq. 2. Then full advantage of a local ALO can be taken. The linearized system may be written as:

$$𝝍_d=𝝍_d^0+\delta 𝝍_d,$$
(4)
where $`𝝍_d^0`$ is the current estimate for the solution vector at depth $`d`$ and $`\delta 𝝍_d`$ is the correction vector to be computed. Using a tri-diagonal $`\mathrm{\Lambda }^{}`$ operator the resulting system for $`\delta 𝝍_d`$ is – like in the classical CL scheme – of block tri-diagonal form coupling each depth point $`d`$ to its nearest neighbors $`d\pm 1`$:

$$𝜸_d\delta 𝝍_{d1}+𝜷_d\delta 𝝍_d+𝜶_d\delta 𝝍_{d+1}=𝐜_d.$$
(5)
The quantities $`𝜶\mathbf{,}𝜷\mathbf{,}𝜸`$ are ($`NN\times NN`$) matrices where $`NN`$ is the total number of physical variables, i.e., $`NN=NL+4`$, and $`𝐜_d`$ is the residual error in the equations. The solution is obtained by the standard Feautrier scheme. As already mentioned, the system Eq. 5 breaks into $`ND`$ independent equations $`\delta 𝝍_d=𝜷_d^1𝐜_d`$ ($`d=1,\mathrm{},ND`$) when a local $`\mathrm{\Lambda }^{}`$ operator is used. The additional numerical effort to set up the subdiagonal matrices and matrix multiplications in the tri-diagonal case is outweighed by the faster global convergence of the ALI cycle, accomplished by the explicit depth coupling in the linearization procedure (Werner 1989).

The principal advantage of the ALI over the CL method becomes clear at this point. Each matrix inversion necessary to solve Eq. 5 requires $`(NL+4)^3`$ operations whereas in the CL method $`(NL+NF+4)^3`$ operations are needed. Since the number of frequency points $`NF`$ is much larger than the number of levels $`NL`$, the matrix inversion in the CL approach is dominated by $`NF`$.

Recent developments concern the problem that the total number of atomic levels tractable in NLTE with the ALI method described so far is restricted to the order of 250, from experience with our model atmosphere code PRO2. This limit is a consequence of the non-linearity of the equations, and in order to overcome it, measures must be taken in order to achieve a linear system whose numerical solution is much more stable. Such a pre-conditioning procedure has been first applied in the ALI context by Werner & Husfeld (1985) using the core saturation method (Rybicki 1972). More advanced work achieves linearity by replacing the $`\mathrm{\Lambda }`$ operator with the $`\mathrm{\Psi }`$ operator (and by judiciously considering some populations as “old” and some as “new” ones within an ALI step) which is formally defined by writing:

$$J_\nu =\mathrm{\Psi }_\nu \eta _\nu ,\mathrm{i}.\mathrm{e}.\mathrm{\Psi }_\nu \mathrm{\Lambda }_\nu /\chi _\nu ,$$
(6)
where the total opacity $`\chi _\nu `$ (as defined in Sect. 3.6.) is calculated from the previous ALI cycle. The advantage is that the emissivity $`\eta _\nu `$ (Sect. 3.6.) is linear in the populations, whereas the source function $`S_\nu `$ is not. Hence the new operator $`\mathrm{\Psi }`$ gives the solution of the transfer problem by acting on a linear function. This idea is based on Rybicki & Hummer (1991) who applied it to the line formation problem, i.e. restricting the set of equations to the transfer and rate equations and regarding the atmospheric structure as fixed. Hauschildt (1993) and Hauschildt & Baron (1999) generalized it to solve the full model atmosphere problem. In addition, splitting the set of statistical equations and solving it separately for each chemical element means that now many hundreds of levels per species are tractable in NLTE.

## 3. Basic equations

We give a short overview of the basic equations to be solved. More detailed presentations may be found in Werner & Dreizler (1999) or Werner et al. (2003).

### 3.1. Radiation transfer

Any numerical method requires a formal solution (i.e. atmospheric structure already given) of the radiation transfer problem. The radiation transfer at any particular depth point can be described by the following equation, formally written for positive and negative $`\mu `$ (which is the cosine of the angle between direction of propagation and outward directed normal to the surface) separately, i.e. for inward and for outward directional intensities $`I`$ with frequency $`\nu `$:

$$\pm \mu \frac{I_\nu (\pm \mu )}{\tau _\nu }=S_\nu I_\nu (\pm \mu ),\mu [0,1].$$
(7)
$`\tau _\nu `$ is the optical depth (which can be defined via the column mass $`m`$ that is used in the other structural equations by $`d\tau _\nu =dm\chi _\nu /\rho `$, with the mass density $`\rho `$) and $`S_\nu `$ is the local source function. Introducing the Feautrier intensity

$$u_{\nu \mu }\left(I_\nu (\mu )+I_\nu (\mu )\right)/2$$
(8)
we obtain the second-order form (Mihalas 1978):

$$\mu ^2\frac{^2u_{\nu \mu }}{\tau _\nu ^2}=u_{\nu \mu }S_\nu ,\mu [0,1].$$
(9)
We may separate the Thomson emissivity term (scattering from free electrons, assumed coherent, with cross-section $`\sigma _e`$) from the source function so that

$$S_\nu =S_\nu ^{}+n_e\sigma _eJ_\nu /\chi _\nu ,$$
(10)
where $`S_\nu ^{}`$ is the ratio of thermal emissivity to total opacity (Sect. 3.6.): $`S_\nu ^{}=\eta _\nu /\chi _\nu `$. Since the mean intensity is the angular integral over the Feautrier intensity the transfer equation becomes

$$\mu ^2\frac{^2u_{\nu \mu }}{\tau _\nu ^2}=u_{\nu \mu }S_\nu ^{}\frac{n_e\sigma _e}{\chi _\nu }_0^1u_{\nu \mu }𝑑\mu .$$
(11)
Thomson scattering complicates the situation by explicit angle coupling. Assuming complete frequency redistribution in spectral lines, no explicit frequency coupling occurs so that the parallel solution for all frequencies enables a very efficient vectorization on the computer.

### 3.2. Statistical equilibrium

#### Rate equations

For each atomic level $`i`$ the rate equation describes the equilibrium of rates into and rates out of this level:

$$n_i\underset{ij}{}P_{ij}\underset{ji}{}n_jP_{ji}=0.$$
(12)
The rate coefficients $`P_{ij}`$ have radiative and collisional components: $`P_{ij}=R_{ij}+C_{ij}`$. The radiative downward rate for example is given by:

$$R_{ji}=\left(\frac{n_i}{n_j}\right)^{}4\pi _0^{\mathrm{}}\frac{\sigma _{ij}(\nu )}{h\nu }\left(\frac{2h\nu ^3}{c^2}+J_\nu \right)e^{h\nu /kT}𝑑\nu .$$
(13)
Photon cross-sections are denoted by $`\sigma _{ij}(\nu )`$. $`(n_i/n_j)^{}`$ is the Boltzmann LTE population ratio. The computation of collisional rates is generally dependent on the specific ion or even transition.

#### Abundance definition equation

The rate equation for the highest level of a given chemical species is redundant. It is replaced by the abundance definition equation. Summation over all levels usually includes not only NLTE levels but also levels which are treated in LTE, according to the specification in the model atom. Denoting the number of ionization stages of species $`k`$ with $`NION(k)`$, the number of NLTE and LTE levels per ion with $`NL(l)`$ and $`LTE(l)`$, respectively, and the number fraction of species $`k`$ with $`y_k`$, we can write:

$$\underset{l=1}{\overset{NION(k)}{}}\left[\underset{i=1}{\overset{NL(l)}{}}n_{kli}+\underset{i=1}{\overset{LTE(l)}{}}n_{kli}^{}\right]=y_k(Nn_e).$$
(14)
#### Charge conservation

We close the system of statistical equations by invoking charge conservation. We denote the total number of chemical species with NATOM, the charge of ion $`l`$ with $`q(l)`$ (in units of the electron charge) and write:

$$\underset{k=1}{\overset{\text{NATOM}}{}}\underset{l=1}{\overset{NION(k)}{}}q(l)\left[\underset{i=1}{\overset{NL(l)}{}}n_{kli}+\underset{i=1}{\overset{LTE(l)}{}}n_{kli}^{}\right]=n_e.$$
(15)
#### Complete statistical equilibrium equations

We introduce a vector comprising the occupation numbers of all NLTE levels, $`𝐧`$ $`=(n_1,\mathrm{},n_{NL})`$. Then the statistical equilibrium equation is written as:

$$\mathrm{𝐀𝐧}=𝐛.$$
(16)
The gross structure of the rate matrix $`𝐀`$ is of block matrix form, because transitions between levels occur within one ionization stage or to the ground state of the next ion. The structure is complicated by ionizations into excited levels and by the abundance definition and charge conservation equations which give additional non-zero elements in the corresponding lines of $`𝐀`$.

### 3.3. Radiative equilibrium

Radiative equilibrium can be enforced by adjusting the temperature stratification either during the linearization procedure or in between ALI iterations. In the former case a linear combination of two different formulations can be used and in the latter case the classical Unsöld-Lucy temperature correction procedure (Lucy 1964) is utilized. The latter is particularly interesting, because it allows to exploit the blocked form of the rate coefficient matrix. This will enable an economic block-by-block solution followed by a subsequent Unsöld-Lucy temperature correction step. On the other side, however, this correction procedure may decelerate the global convergence behavior of the ALI iteration.

The two forms of expressing the radiative equilibrium condition follow from the postulation that the energy emitted by a volume element per unit time is equal to the absorbed energy per unit time (integral form):

$$_0^{\mathrm{}}\chi _\nu (S_\nu J_\nu )𝑑\nu =0.$$
(17)
This formulation is equivalent to invoking flux constancy throughout the atmosphere (differential form) involving the nominal flux $``$:

$$_0^{\mathrm{}}\frac{}{\tau _\nu }(f_\nu J_\nu )𝑑\nu =0,$$
(18)
where $`f_\nu `$ is the variable Eddington factor, defined as

$$f_\nu =_0^1\mu ^2u_{\nu \mu }𝑑\mu /_0^1u_{\nu \mu }𝑑\mu $$
(19)
and computed from the Feautrier intensity $`u_{\nu \mu }`$ (Eq. 8) after the formal solution. The differential form is more accurate at large depths, while the integral form behaves numerically better at small depths. Instead of arbitrarily selecting that depth in the atmosphere where we switch from one formulation to the other, we use a linear combination of both constraint equations which guarantees a smooth transition with depth, based on physical grounds:

$$\frac{1}{\overline{\kappa }_J}_0^{\mathrm{}}\chi _\nu (S_\nu J_\nu )𝑑\nu +\overline{\mathrm{\Lambda }}_J^{}_0^{\mathrm{}}\frac{}{\tau _\nu }(f_\nu J_\nu )𝑑\nu \overline{\mathrm{\Lambda }}_J^{}F_0=0.$$
(20)
For details on this equation and on our implementation of the Unsöld-Lucy procedure see Dreizler (2003).

### 3.4. Hydrostatic equilibrium

The total pressure $`P`$ comprises gas, radiation and turbulent pressures, so that we write the equation for hydrostatic equilibrium as:

$$\frac{d}{dm}P=\frac{d}{dm}\left(NkT+\frac{4\pi }{c}_0^{\mathrm{}}f_\nu J_\nu 𝑑\nu +\frac{1}{2}\rho v_{\mathrm{turb}}^2\right)=g$$
(21)
with Boltzmann’s constant $`k`$ and the turbulent velocity $`v_{\mathrm{turb}}`$. The hydrostatic equation may either be solved simultaneously with all other equations or separately in between iterations. The overall convergence behavior is usually the same in both cases. If taken into the linearization scheme and a local $`\mathrm{\Lambda }^{}`$ operator is used then, like in the case of the radiative equilibrium equation, explicit depth coupling enters via the depth derivative $`d/dm`$. Again, solution of the linearized equations has to proceed inwards starting at the outer boundary.

### 3.5. Particle conservation and fictitious massive particle density

The total particle density $`N`$ is the sum of electron density plus the population density of all atomic states, LTE and NLTE levels:

$$N=n_e+\underset{k=1}{\overset{\text{NATOM}}{}}\underset{l=1}{\overset{NION(k)}{}}\left[\underset{i=1}{\overset{NL(l)}{}}n_{kli}+\underset{i=1}{\overset{LTE(l)}{}}n_{kli}^{}\right].$$
(22)
A fictitious massive particle density $`n_H`$ is defined by:

$$n_H=\underset{k=1}{\overset{\text{NATOM}}{}}m_k\underset{l=1}{\overset{NION(k)}{}}\left[\underset{i=1}{\overset{NL(l)}{}}n_{kli}+\underset{i=1}{\overset{LTE(l)}{}}n_{kli}^{}\right].$$
(23)
The mass of a chemical species in AMU is denoted by $`m_k`$. Introducing the mass of a hydrogen atom $`m_H`$, we then may simply write for the material density

$$\rho =n_Hm_H.$$
(24)
### 3.6. Opacity and emissivity

Thermal opacity and emissivity are made up by radiative bound-bound, bound-free and free-free transitions. For each species we compute and sum up:

$`\kappa _\nu `$ $`=`$ $`{\displaystyle \underset{l=1}{\overset{NION}{}}}[{\displaystyle \underset{i=1}{\overset{NL(l)}{}}}{\displaystyle \underset{j>i}{\overset{NL(l)}{}}}\sigma _{lilj}(\nu )(n_{li}n_{lj}{\displaystyle \frac{g_{li}}{g_{lj}}}e^{h(\nu \nu _{ij})/kT})`$
$`+{\displaystyle \underset{i=1}{\overset{NL(l)}{}}}{\displaystyle \underset{j>i}{\overset{NL(l+1)}{}}}\sigma _{lil+1,k}(\nu )\left(n_{li}n_{li}^{}e^{h\nu /kT}\right)`$
$`+n_e\sigma _{kk}(l,\nu )(1e^{h\nu /kT})({\displaystyle \underset{i=1}{\overset{NL(l+1)}{}}}n_{l+1,i}+{\displaystyle \underset{i=1}{\overset{LTE(l+1)}{}}}n_{l+1,i}^{})]`$
where the total opacity includes Thomson scattering, i.e. $`\chi _\nu =\kappa _\nu +n_e\sigma _e`$, and

$`{\displaystyle \frac{\eta _\nu }{2h\nu ^3/c^2}}`$ $`=`$ $`{\displaystyle \underset{l=1}{\overset{NION}{}}}[{\displaystyle \underset{i=1}{\overset{NL(l)}{}}}{\displaystyle \underset{j>i}{\overset{NL(l)}{}}}\sigma _{lilj}(\nu )n_{lj}{\displaystyle \frac{g_{li}}{g_{lj}}}e^{h(\nu \nu _{ij})/kT}`$
$`+{\displaystyle \underset{i=1}{\overset{NL(l)}{}}}{\displaystyle \underset{j>i}{\overset{NL(l+1)}{}}}\sigma _{lil+1,k}(\nu )n_{li}^{}e^{h\nu /kT}`$
$`+n_e\sigma _{kk}(l,\nu )e^{h\nu /kT}({\displaystyle \underset{i=1}{\overset{NL(l+1)}{}}}n_{l+1,i}+{\displaystyle \underset{i=1}{\overset{LTE(l+1)}{}}}n_{l+1,i}^{})].`$
$`\sigma _{lil+1,k}(\nu )`$ denotes the cross-section for photoionization from level $`i`$ of ion $`l`$ into level $`k`$ of ion $`l+1`$. The double summation over the bound-free continua takes into account the possibility that a particular level may be ionized into more than one level of the next high ion. The source function used for the approximate radiation transfer is $`\eta _\nu /\kappa _\nu `$, thus, excludes Thomson scattering. For the exact formal solution of course, the total opacity $`\chi _\nu `$ in the expression Eq. 10 includes the Thomson term ($`n_e\sigma _e`$).

## 4. The Accelerated Lambda Iteration (ALI)

In all constraint equations described above the mean intensities $`J_\nu `$ are substituted by the approximate radiation field Eq. 2 in order to eliminate these variables from the solution vector Eq. 1. In principle the ALO may be of arbitrary form as long as the iteration procedure converges. In practice however an optimum choice is desired in order to achieve convergence with a minimum amount of iteration steps.

In the case of diagonal (local) lambda operators the mean intensity $`J_d`$ at a particular depth $`d`$ in the current iteration step is computed solely from the local source function $`S_d`$ and a correction term $`\mathrm{\Delta }J_d`$, the latter involving the source functions (of all depths) from the previous iteration. Dropping the iteration count and introducing indices denoting depth points we can rewrite Eq. 2:

$$J_d=\mathrm{\Lambda }_{d,d}^{}S_d+\mathrm{\Delta }J_d.$$
(27)
In the discrete form we now think of $`\mathrm{\Lambda }^{}`$ as a matrix acting on a vector whose elements comprise the source functions of all depths. Then $`\mathrm{\Lambda }_{d,d}^{}`$ is the diagonal element of the $`\mathrm{\Lambda }^{}`$ matrix corresponding to depth point $`d`$. Writing $`\mathrm{\Lambda }_{d,d}^{}B_d`$ we have a purely local expression for the mean intensity:

$$J_d=B_dS_d+\mathrm{\Delta }J_d.$$
(28)
Much better convergence is obtained if the mean intensity is computed not only from the local source function but also from the source function of the neighboring depths points. Then the matrix representation of $`\mathrm{\Lambda }^{}`$ is of tri-diagonal form and we may write:

$$J_d=C_{d1}S_{d1}+B_dS_d+A_{d+1}S_{d+1}+\mathrm{\Delta }J_d,$$
(29)
where $`C_{d1}`$ and $`A_{d+1}`$ represent the upper and lower subdiagonal elements of $`\mathrm{\Lambda }^{}`$ and $`S_{d\pm 1}`$ the source functions at the adjacent depths. We emphasize again that the actual source functions in Eq. 29 are computed from the actual population densities and temperature which are unknown. We therefore have a non-linear set of equations which is solved by either a Newton-Raphson iteration or other techniques, resulting in the solution of a tri-diagonal linear equation of the form Eq. 5. It was shown in Olson et al. (1986) that the elements of the optimum $`\mathrm{\Lambda }^{}`$ matrix are given by the corresponding elements of the exact $`\mathrm{\Lambda }`$ matrix.

## 5. Solution of the non-linear equations by iteration

The complete set of non-linear equations for a single iteration step Eq. 2 comprises at each depth the equations for statistical, radiative, and hydrostatic equilibrium and the particle conservation equation. For the numerical solution we introduce discrete depth and frequency grids. The equations are then linearized and solved by a suitable iterative scheme. Explicit angle dependency of the radiation field is not required here and consequently eliminated by the use of variable Eddington factors. Angle dependency is only considered in the formal solution of the transfer equation.

### 5.1. Discretization and Linearization

A depth grid is set up and we start from a gray approximation computing a LTE continuum model using the Unsöld-Lucy temperature correction procedure. Depth points (typical number is 90) are set equidistantly on a logarithmic (Rosseland) optical depth scale. Our NLTE code uses the column mass as an independent depth variable.

The frequency grid is established based upon the atomic data input file. Frequency points are set blue- and redward of each absorption edge and for each spectral line. Gaps are filled in by setting continuum points. Finally, the quadrature weights are computed. Frequency integrals appearing e.g. in Eq. 20 are replaced by quadrature sums and differential quotients involving depth derivatives by difference quotients.

All variables $`x`$ are replaced by $`xx+\delta x`$ where $`\delta x`$ denotes a small perturbation of $`x`$. Terms not linear in these perturbations are neglected. The perturbations are expressed by perturbations of the basic variables:

$$\delta x=\frac{x}{T}\delta T+\frac{x}{n_e}\delta n_e+\frac{x}{N}\delta N+\frac{x}{n_H}\delta n_H+\underset{l=1}{\overset{NL}{}}\frac{x}{n_l}\delta n_l.$$
(30)
### 5.2. Newton-Raphson iteration and alternative methods

As described in Sect. 2. the linearized equations have a tri-diagonal block-matrix form, see Eq. 5. Inversion of the grand matrix ($`𝐓`$ sized $`(NNND)\times (NNND)`$, i.e. about $`10^4\times 10^4`$ in typical applications) is performed with a block-Gaussian elimination scheme, which means that our iteration of the non-linear equations represents a multi-dimensional Newton-Raphson method. The problem is structurally simplified when explicit depth coupling is avoided by the use of a local ALO, however, the numerical effort is not much reduced, because in both cases the main effort lies with the inversion of matrices sized $`NN\times NN`$. The Newton-Raphson iteration involves two numerically expensive steps, first setting up the Jacobian (comprising $`𝜶\mathbf{,}𝜷\mathbf{,}𝜸`$) and then inverting it. Additionally, the matrix inversions necessary to solve Eq. 5 limit their size to about $`NN=250`$ because otherwise numerical accuracy is lost.

Two variants recently introduced in stellar atmosphere calculations are able to improve both, numerical accuracy and, most of all, computational speed. Broyden’s (1965) variant belongs to the family of quasi-Newton methods and it was first used in model atmosphere calculations in Dreizler & Werner (1991), Hamann et al. (1991) and Koesterke et al. (1992). It avoids the repeated set-up of the Jacobian by the use of an update formula. On top of this, it also gives an update formula for the inverse Jacobian. Another variant, the Kantorovich method was introduced into model atmosphere calculations by Hubeny & Lanz (1992). This method simply keeps fixed the Jacobian during the linearization cycle and it is surprisingly stable.

## 6. NLTE metal line blanketing

Despite the capacity increase for the NLTE treatment of model atmosphere problems by introducing the ALI method combined with pre-conditioning techniques, the blanketing by millions of lines from the iron group elements arising from transitions between some $`10^5`$ levels could only be attacked with the help of statistical methods. These have been introduced into NLTE model atmosphere work by Anderson (1989). At the outset, model atoms are constructed by combining many thousands of levels into a relatively small number of superlevels which can be treated with ALI (or other) methods. Then, in order to reduce the computational effort, two approaches were developed which vastly decrease the number of frequency points (and hence the number of transfer equations to be solved) to describe properly the complex frequency dependence of the opacity. These two approaches have their roots in LTE modeling techniques, where for the same reason statistical methods are applied for the opacity treatment: The Opacity Distribution Function (ODF) and Opacity Sampling (OS) approaches. Both are based on the circumstance that the opacity (in the LTE approximation) is a function of two only local thermodynamic quantities. Roughly speaking, each opacity source can be written in terms of a population density and a photon cross-section for the respective radiative transition: $`\kappa _\nu n_l\sigma _{lu}(\nu )`$. In LTE the population follows from the Saha-Boltzmann equations, hence $`n_l=n_l(n_e,T)`$. The OS and ODF methods use such pre-tabulated (on a very fine frequency mesh) $`\kappa _\nu (n_e,T)`$ during the model atmosphere calculations. The NLTE situation is more complicated, because pre-tabulation of opacities is not useful. The population densities at any depth now also depend explicitly on the radiation field (via the rate equations which substitute the TE Saha-Boltzmann statistics) and thus on the populations in each other depth of the atmosphere. As a consequence, the OS and ODF methods are not applied to opacity tabulations, but on tabulations of the photon cross-sections $`\sigma (\nu )`$. These do depend on local quantities only, e.g. line broadening by Stark and Doppler effects is calculated from $`T`$ and $`n_e`$. In the NLTE case the cross-section takes over the role which the opacity played in the LTE case. More details on this topic can be found in Werner & Dreizler (1999).

Our strategy is the following. Before any model atmosphere calculation is started, the atomic data are prepared by constructing superlevels, and the cross-sections for superlines. Then these cross-sections are either sampled on a coarse frequency grid or ODFs are constructed. These data are put into the model atom which is read by the model code.

## 7. Summary

The construction of metal line blanketed models in hydrostatic and radiative equilibrium under NLTE conditions was the last and long-standing problem of classical model atmosphere theory and it is finally solved with a high degree of sophistication. The essential milestones for this development, starting from the pioneering work of Auer & Mihalas (1969) are:

* Introduction of the Accelerated Lambda Iteration (ALI, or “operator splitting” methods), based upon early work by Cannon (1973) and Scharmer (1981). First ALI model atmospheres were constructed by Werner (1986).
* Introduction of statistical approaches to treat the iron group elements in NLTE by Anderson (1989).
* Linear formulation of the statistical equations (Rybicki & Hummer 1991, Hauschildt 1993).
* Computation of atomic data by Kurucz (1991), by the Opacity Project (Seaton et al. 1994) and the Iron Project (Hummer et al. 1993).

##### Acknowledgments.

Stellar atmosphere work in Tübingen is supported by DFG and DLR. I thank for financial support from the conference organizers.

## References

Anderson, L. S. 1989, ApJ, 339, 558

Auer, L. H., & Mihalas, D. 1969, ApJ, 158, 641

Broyden, C. G. 1965, Math. Comp., 19, 577

Cannon, C. J. 1973, ApJ, 185, 621

Dreizler, S. 2003, in Stellar Atmosphere Modeling, ed. I. Hubeny, D. Mihalas, & K. Werner, ASP Conf. Series, in press

Dreizler, S., & Werner, K. 1991, in Stellar Atmospheres: Beyond Classical Models, NATO ASI Series C, Vol. 341, ed. L. Crivellari, I. Hubeny & D.G. Hummer, (Dordrecht: Kluwer), 155

Hamann, W.-R., Koesterke, L., & Wessolowski, U. 1991, in Stellar Atmospheres: Beyond Classical Models, NATO ASI Series C, Vol. 341, ed. L. Crivellari, I. Hubeny & D.G. Hummer, (Dordrecht: Kluwer), 69

Hauschildt, P. H. 1993, JQSRT, 50, 301

Hauschildt, P. H., & Baron, E. 1999, in Computational Astrophysics, ed. H. Riffert & K. Werner, J. of Comp. and Applied Mathematics, 109, 41

Hubeny, I., & Lanz, T. 1992, A&A, 262, 501

Hummer, D. G., Berrington, K. A., Eissner, W., Pradhan, A. K., Saraph, H. E., & Tully, J. A. 1993, A&A, 279, 298

Koesterke, L., Hamann, W.-R., & Kosmol, P. 1992, A&A, 255, 490

Kurucz, R. L. 1991, in Stellar Atmospheres: Beyond Classical Models, NATO ASI Series C, Vol. 341, ed. L. Crivellari, I. Hubeny & D.G. Hummer, (Dordrecht: Kluwer), 441

Lucy, L. B. 1964, in First Harvard Smithsonian Conference on Stellar Atmospheres, SAO Special Report No. 167, Cambridge, Massachusetts, 93

Mihalas, D. 1978, Stellar Atmospheres, 2nd edition, (San Francisco: Freeman)

Olson, G. L., & Kunasz, P. B. 1987, JQSRT, 38, 325

Olson, G. L., Auer, L. H., & Buchler, J. R. 1986, JQSRT, 35, 431

Rybicki, G. B. 1972 in Line Formation in the Presence of Magnetic Fields, ed. R. G. Athay, L. House & G. Newkirk, Jr. (Boulder: High Altitude Observatory), 145

Rybicki, G. B., & Hummer, D. G. 1991, A&A, 245, 171

Scharmer, G. B. 1981, ApJ, 249, 720

Seaton, M. J., Yan, Y., Mihalas, D., & Pradhan, A. K. 1994, MNRAS, 266, 805

Werner, K. 1986, A&A, 161, 177

Werner, K. 1989, A&A, 226, 265

Werner, K., & Husfeld, D. 1985, A&A, 148, 417

Werner, K., & Dreizler, S. 1999, in Computational Astrophysics, ed. H. Riffert & K. Werner, J. of Comp. and Applied Math., 109, 65

Werner, K., Deetjen, J.L., Dreizler, S., Nagel, T., Rauch, T., & Schuh, S.L. 2003 in Stellar Atmosphere Modeling, ed. I. Hubeny, D. Mihalas, & K. Werner, ASP Conf. Series, in press


text:
# New Evidence of Discrete Scale Invariance in the Energy Dissipation of Three-Dimensional Turbulence: Correlation Approach and Direct Spectral Detection

## I  Introduction

The small-scale physical quantity that has received the most attention in the context of a phenomenological description of hydrodynamical turbulence is the dissipation rate of kinetic energy Frisch96 . Richardson’s picture of turbulent cascades Rich22 , in which large eddies break down into smaller ones receiving a certain fraction of the flux of kinetic energy from larger scales, is a multiplicative process. The hypothesis that eddies of any generation are space filling led to the famous K41 theory K41a ; K41b . A conceptually appealing view, dating back to Obukhov O62  and Kolmogorov K62 , visualizes the transfer of kinetic energy as a self-similar multiplicative process. This view is still the point of departure of many phenomenological intermittency models today. The thought of the general theory of multifractals was proposed by Mandelbrot M74  in the study of certain high-order moments with a general cascade model of random curling. A phenomenological model of intermittency, referred to as the $`\beta `$-model beta78 , analyzed the multiplicative cascade of the breakup of eddies without space filling leading to an energy transfer to a smaller and smaller fraction of the total space. Benzi, Paladin, Parisi & Vulpiani BPPV84  and Frisch & Parisi FP85  have argued for a singular structure and for a multifractal nature of the energy dissipation in three-dimensional fully-developed turbulence. Meneveau & Sreenivasan MS87  presented a simple multifractal model, the binomial model for fully developed turbulence based on the two-scale Cantor set, which is known as the $`p`$-model. This model and the specific multifractal facets of turbulent energy dissipation were investigated experimentally by MS91 , using hot wire measurement technique. These works clarified the phenomenological nature of the multifractality of turbulent energy dissipation. In particular, the view of the energy cascade as belonging to the general class of multiplicative processes opens up the possibility that the scale invariant properties reported for turbulence flows might be broken in part into a discrete scale invariance, leading to the observation of log-periodic corrections to scaling SorDSI .

Based on theoretical argument and experimental evidence, Sornette Sor98TurbDSI  conjectured that structure functions of turbulent time series may exhibit log-periodic modulations decorating their power law dependence. He stressed the need for novel methods of averaging and proposed a novel “canonical” averaging scheme for the analysis of structure factors of turbulent flows in order to provide convincing experimental evidence. The strategy was proposed to determine the sample-specific control parameter $`r_c`$ and translate it into a specific “phase” in the logarithm of the scale $`\mathrm{ln}r`$ which, when used as the origin, allows one to rephase the different measurements of a structure factor in different realizations.

A first systematic search of discrete scale invariance in turbulence was carried out in freely decaying two-dimensional turbulence Johansen2000b . The number of vortices, their radius and separation were found to display log-periodic oscillations as a function of time with an average log-frequency of $`45`$ corresponding to a preferred scaling ratio of $`1.21.3`$. Most recently, significant evidence for the presence of intermittent discretely self-similar cascades was obtained from the analysis of the moments of the energy dissipation rate in three-dimensional fully developed turbulence TurbDSI . The results of spectral analyses of the canonically averaged local scaling exponents of high moments of the energy dissipation rate obtained in TurbDSI  are summarized in Fig. 1. The fundamental log-frequency is estimated to be $`f1.44`$ which suggests a preferred scaling ratio of $`\gamma =e^{1/f}2`$. In TurbDSI , a simple multifractal model was introduced to explain the fact that the log-frequency is independent of the order of the moments.

As was explained in TurbDSI , the difficulty in detecting log-periodic oscillations (the hallmark of discrete scale invariance) in various functions such as moments of the energy dissipation rate stems from the fact that, while the physics of the cascade may be universal with an universal scaling factor $`\gamma `$, the nucleating scales (translating in phases of reference in the log-periodicity analysis) have to change from one realization to another as the cascade may be dynamically triggered from a variety of scales $`\mathrm{}_a`$ from the integral scale and deep in the inertial range. Therefore, the corresponding phases are expected to be non-universal. Averaging a signal over several such transient cascades will thus wash out the information on log-periodicity by “destructive interferences” of the log-periodic oscillations, giving them the appearance of noise.

The previous analyses of real 2-D and 3-D turbulence data addressed this problem by using two versions of the canonical averaging approach, devised precisely to rephase different realizations and thus alleviate the problem of the random phases. Using $`n`$ realizations of the data, the first version consists in performing a canonical averaging on the $`n`$ samples to get an average realization and then analyzes its spectral properties; the second version performs a spectral analysis of each of the $`n`$ samples and then averages the $`n`$ periodograms. In order to perform these averages, we need to have different realizations. For the turbulence data analyzed in TurbDSI  and revisited here, a long time series of velocity increments is divided in many samples, each lasting about $`1/5`$ of a turn-over time, so that we can avoid scrambling of the phases by a suitable canonical averaging as explained below.

The goal of the present note is to revisit these two schemes and propose two novel and simpler versions. By simplifying the analysis and by providing independent implementations, we thus confirm further the presence of log-periodicity in the moments of the energy dissipation rate, and thus the existence of discrete scale invariance in 3D fully-developed hydrodynamical turbulence.

## II  Rephasing of log-periodicity by the cross-correlation method

In the first scheme of canonical averaging, the rephasing operation is the central part. In TurbDSI , the so-called central maxima approach was applied to the detrended logarithm of the moment $`M_q(\mathrm{}/\eta )`$ of the dissipation rate as a function of the logarithm $`\mathrm{ln}(\mathrm{}/\eta )`$ of the scale $`\mathrm{}`$ (in units of the Kolmogorov microscale $`\eta =0.195`$ mm TurbDSI ): it consists in finding the maximum of the detrended logarithm of the moment closest to the middle point of a single sample and set the logarithm of the corresponding scale (its abscissa) as the origin of the logarithm of the scales. The new abscissa, called $`\mathrm{\Delta }`$, becomes sample dependent and is used to superimpose and then average the detrended moments over the different samples. This approach was found to work satisfactorily except that the canonically averaged series present decreasing amplitudes with increasing $`|\mathrm{\Delta }|`$.

The central maxima rephasing approach used in TurbDSI  determines an absolute sample-dependent phase which is specific to each single sample and then rephases all samples. The alternative idea explored here is based on the remark that the detection of phases can be performed in relative value, that is, by estimating the phase shift between all possible pairs of samples. The cross-correlation function is used to calculate the phase-shift between two samples Nikias . This relative phase-shift approach thus uses the calculation of all phase-shifts between pairs of samples to rephase all samples by choosing arbitrarily one of them as the reference.

As in TurbDSI , we use the high-Reynolds turbulence longitudinal velocity data collected at the S1 ONERA wind tunnel by the Grenoble group from LEGI Anselmet1984 . We use a set of 20 records, each containing $`320`$ samples in the inertia range, each sample having $`2^{11}`$ points. These $`2^{11}`$ points correspond approximately to one-fifth of a turn-over time TurbDSI . We calculate the local moment exponent $`\tau (q,\mathrm{ln}(\mathrm{}/\eta ))`$ as the logarithmic derivative of $`M_q(\mathrm{}/\eta )`$ with respect to $`\mathrm{}/\eta `$:

$$\tau (q,\mathrm{ln}(\mathrm{}/\eta ))=\frac{d\mathrm{ln}[M_q(\mathrm{}/\eta )]}{d\mathrm{ln}(\mathrm{}/\eta )}.$$
(1)
We shall investigate in the sequel $`\mathrm{}/\eta [230,1843]`$ which is well within the inertial range. To obtain $`\tau (q,\mathrm{ln}(\mathrm{}/\eta ))`$ from (1), we used a smoothing low-pass digital filter, namely, the Savitzky-Golay filter. This filter approximates the underlying function in a moving window by means of polynomials of a low degree. The filter coefficients are choosen so that, unlike other filters, several low-order moments of the filtered function are preserved, which is a very desirable property (see Press1996  for details). The Savitzky-Golay filters are applied to $`\mathrm{ln}[M_q(\mathrm{}/\eta )]`$ expressed as a function of $`\mathrm{ln}(\mathrm{}/\eta )`$ and the local derivative is given from the the analytic derivative of the fitted polynomial Press1996 . There are two parameters in the Savitzky-Golay filters, i.e., the width of the sliding window $`N_L+N_R+1`$ where $`N_L`$ and $`N_R`$ are the numbers of points to the left and to the right of the investigating point and the order $`M`$ of fitting polynomials. We used symmetric windows such that $`N_L=N_R`$. A total of 24 filters were applied to test the robustness of the results with respect to the choice of the filters: $`N_L=5,\mathrm{},10`$ and $`M=4,\mathrm{},7`$.

We recall that, in the central maxima approach of TurbDSI , one identifies the sample-dependent reduced scale $`\mathrm{}_c/\eta `$ closest to the central point of the scale-interval at which $`\tau (q,\mathrm{ln}(\mathrm{}/\eta ))`$ is maximum for each sample. In contrast, the correlation approach takes an arbitrary sample as the reference sample and evaluates the phase shift measured by the cross-correlation function of all other samples relative to this reference. We then rephase all functions $`\tau (q,\mathrm{ln}(\mathrm{}/\eta ))`$ of the variable $`\mathrm{ln}(\mathrm{}/\eta )`$ so that the cross-correlation function of all sample pairs becomes peaked at zero lag, after the rephasing has been performed.

We then average these rephased functions $`\tau (q,\mathrm{\Delta })`$ over the $`320`$ samples. This defines the local average exponent

$$D(\mathrm{\Delta })=\tau (q,\mathrm{\Delta }),$$
(2)
where the average is performed over the different samples at fixed $`\mathrm{\Delta }=\mathrm{ln}(\mathrm{}/\eta )\mathrm{ln}(\mathrm{}_c/\eta )`$ (where $`\mathrm{ln}(\mathrm{}_c/\eta )`$ vary from sample to sample and $`\mathrm{ln}(\mathrm{}/\eta )`$ is adjusted accordingly to ensure the re-phasing). Fig. 2 shows four examples of $`D(\mathrm{\Delta })`$ as a function of $`\mathrm{\Delta }`$ based on the correlation rephasing approach for four different filters with $`N_L=6`$ and $`M=4,5,6,7`$ respectively. Notice that the amplitude of the oscillations decreases with decreasing $`\mathrm{\Delta }`$. This is in line with the concept TurbDSI  that, the smaller is the scale $`\mathrm{}`$, the more complex is the superposition of intermittent cascades as the number of starting scales for the cascade proliferate when the nucleating scale departs from the integral scale. This should lead to a stronger self-averaging property, as we observe here when $`\mathrm{\Delta }`$ decreases. Furthermore, since the samples have finite sizes, $`D(\mathrm{\Delta })`$ is well-defined only for not too large $`|\mathrm{\Delta }|`$ because there are fewer samples involved in the averaging for larger $`|\mathrm{\Delta }|`$. In practice, we will restrict our analysis to $`|\mathrm{\Delta }|`$ smaller than a threshold $`\mathrm{\Delta }_m`$ taken here approximately equal to $`1.0`$, as suggested from a visual inspection of Fig. 2.

The next step in the analysis consists in performing a spectral analysis of $`D(\mathrm{\Delta })`$. We used a kind of spectral analysis known under the name of the Lomb spectrum (see Press1996  for details). It is designed to analyze a non-uniformly sampled functions, fitting them to a sine function of varying frequency on a non-uniform grid by means of the least square method. This provides the analog of a standard periodogram. We thus obtain 24 such Lomb periodograms for each record (of $`320`$ samples), one for each filter. We then repeat the above steps for all the other 19 records of of $`320`$ samples each and we thus obtain a total of 480 Lomb periodograms. Alternatively, we have grouped the 20 records into a single large time series with 6400 samples on which the canonical averaging is performed directly. This increased statistics improves the significance level of the Lomb peaks in the resultant periodograms. The comparison of our method on different partitions of the data allows us to test the stability and robustness of the results.

Having computed the 480 Lomb periodograms (corresponding to 24 Lomb periodograms for each of the 20 records), we perform an average over all these 480 periodograms. This analysis was carried out for three different choices of the threshold $`\mathrm{\Delta }_m=1.0,0.8,0.6`$. The averaged Lomb powers are shown in Fig. 3. All the three curves have significant peaks at $`f_2=2/\mathrm{ln}2`$ and $`f_3=3/\mathrm{ln}2`$, which suggest a fundamental log-frequency $`f=1/\mathrm{ln}2`$. The periodogram with $`\mathrm{\Delta }_m=1.0`$ also possesses many spurious weaker peaks. The periodogram with $`\mathrm{\Delta }_m=0.8`$ exhibits the highest Lomb peaks and clearest harmonics. The difference in peak heights for $`\mathrm{\Delta }_m=0.8`$ compared with $`\mathrm{\Delta }_m=0.6`$ is partially explained by the fact that the former has more data points. It is indeed known that the Lomb power is approximately proportional to the number of analyzed points LombTests . The results summarized in Figure 3 thus provide a strong confirmation of the analysis of TurbDSI  on the presence of discrete scale invariance.

## III Direct detection of log-periodicity by Fast Fourier Transform

We now turn to the second method of implementation of the canonical averaging discussed in the introduction. In a nutshell, the methods consists in performing directly a spectral analysis of each sample and then average over the obtained spectra. In order to improve the signal-to-noise, we first filter the high-frequencies of the function $`\mathrm{ln}[M_q(\mathrm{}/\eta )]`$ of the variable $`\mathrm{ln}(\mathrm{}/\eta )`$ by performing an interpolation of the slightly unevenly sampled third-order moment values $`M_3(\mathrm{}/\eta )`$ with cubic splines. This provides us with equally sampled points that are required for the use of a fast Fourier transform (FFT). We then define the first difference and subtract the mean which defines the function $`\delta \mathrm{ln}[M_3(\mathrm{}/\eta )]`$ The periodogram of $`\delta \mathrm{ln}[M_3(\mathrm{}/\eta )`$\] for each sample is determined via a fast Fourier transform. We then perform an averaging over all the periodograms.

Specifically, the calculation of the third-order moment $`M_3(\mathrm{}/\eta )`$ is done as in Sec. II with the same data sets. We use one record with 320 samples. The number of the evenly sampled points generated by the spline interpolation is 257, such that we will have 256 points at which the function $`\delta \mathrm{ln}[M_3(\mathrm{}/\eta )]`$ is sampled. Fig. 4 shows a realization of $`M_3(\mathrm{})`$ for one specific sample. For this, $`\delta _i\mathrm{ln}[M_3(\mathrm{}/\eta )]`$ is obtained as

$$\delta _i\mathrm{ln}[M_3(\mathrm{}/\eta )]=\mathrm{ln}[M_3(\mathrm{}_{i+1}/\eta )]\mathrm{ln}[M_3(\mathrm{}_i/\eta )]\mu ,$$
(3)
where $`i=1,\mathrm{},256`$ and $`\mu =\mathrm{ln}[M_3(\mathrm{}_{j+1}/\eta )]\mathrm{ln}[M_3(\mathrm{}_j/\eta )]_j`$. The resulting detrended third-order moment is shown in Fig. 5 and its periodogram is depicted in Fig. 6. The vertical lines in Fig. 6 indicate log-frequencies equal to $`k/\mathrm{ln}2`$ with $`k=1,2,\mathrm{},21`$. One can observe well-defined peaks close to the frequencies corresponding to $`k=1,2,3,4,5,6,9,12,15,16,17,21`$. However, there are also peaks in between these vertical lines. Therefore, the signature of log-periodicity in Figs. 4-6 is ambiguous.

We then average over the 320 periodograms, one for each sample, and get Fig. 7. Most peaks at log-frequencies other than $`k/\mathrm{ln}2`$ have been washed out by the averaging procedure. The vertical dashed lines correspond to log-frequencies $`f_k=k/\mathrm{ln}2`$ with $`k=1,\mathrm{},21`$, which are the harmonics of a fundamental log-frequency $`1/\mathrm{ln}2`$ corresponding to a preferred scaling ratio equal to $`\gamma =2`$. The local peaks defined preferred frequencies $`F_m`$ that are shown as open circles. Most of the peaks ($`m=3,\mathrm{},18,20,22,24,27`$) in Fig. 7 are very close to the vertical dotted lines which can be interpreted as harmonics of $`f=1/\mathrm{ln}2`$. The background noise has a bell-shape. Comparing Fig. 7 with Fig. 6, the averaging of the spectra over the 320 samples has suppressed peaks that can be interpreted as noise and has improved remarkably the significance level of the other Lomb peaks.

In order to estimate as accurately as possible the numerical value of the fundamental log-frequency, Fig. 8 plots with open circles the log-frequencies $`F_m`$ of all the peaks found in Fig. 7 as a function of the peak sequence number $`m`$. One observes a straight line with slope $`f=1.420`$ covering 16 harmonics. This value is very close to the theoretic ansatz $`f=1/\mathrm{ln}21.4427`$. An estimation of the uncertainty in the determination of $`f`$ is obtained by removing one point in each data set, and redoing the whole procedure. Then, the standard deviation calculated over these samples with one removed point provides a quantification of the error in the determination of $`f`$. In order to include more points in the fit, we collected peaks near the vertical lines in Fig. 7 with log-frequencies $`f_k`$ and plotted $`f_k`$ with open squares in Fig. 8 as well. Note that $`f_k`$ with $`k=1,\mathrm{},20`$ corresponds to peaks enumerated $`m=3,\mathrm{},18,20,22,24,27`$ mentioned above. Again, a very good linear fit is obtained with slope $`f=1.427`$, which provides a slightly better estimate of $`1/\mathrm{ln}\gamma `$. We estimated the slopes of 19 other experimental records which, together with the slope from the first record, give the estimate of the fundamental log-frequency $`f=1.434\pm 0.007`$. An alternative estimate of $`f`$ is the mean of the ratios $`f_k/k`$. This gives $`f=1.455\pm 0.023`$. In this estimation, the first point $`(k=1,f_k=1.864)`$ was discarded.

Overall, we conclude that this new analysis agrees with and strengthens strongly previous results obtained in Sec. II and in TurbDSI  on the observation of log-periodicity. It is also quite surprising and somewhat gratifying to observe so many harmonics.

## IV Conclusion

In this paper, we have provided two novel analyses for the search of log-periodic corrections to scaling in the moments of energy dissipation of three-dimensional fully developed turbulence. The first analysis introduces a pairwise cross-correlation approach for the determination of a phase-shift, which is an alternative of the central maxima rephasing approach developed in TurbDSI  and is more natural and straightforward to implement. The resulting averaged Lomb periodogram have two significant peaks near log-frequencies equal to $`2/\mathrm{ln}2`$ and $`3/\mathrm{ln}2`$ and low peaks near other $`k/\mathrm{ln}2`$ with $`k=1,4,5`$. We have checked that these results are robust when varying the range of scales over which the analysis is performed.

The second method uses a more straightforward sample-by-sample spectral analysis, with a simple cubic spline interpolation (rather than the more delicate Savitzky-Golay filters) and a standard FFT. Averaging over the hundreds of spectra provides a strong log-periodic signals with more than 20 harmonics of what appears to be a fundamental logarithmic frequency $`f=1.434\pm 0.007`$. The corresponding preferred scaling ratio is thus estimated to be $`\gamma =2.008\pm 0.006`$. This novel analysis lifts possible remaining worries that one could keep that the results of TurbDSI  could be somehow due to the use of the Savitzky-Golay filters.

The present work further confirms the presence of log-periodicity in the energy transfer of fully developed turbulence, as well as provide additional evidence of the presence of multifractality in the energy dissipation MS91 . However, the preferred scaling ratio $`\gamma 2`$ found here does not necessarily imply that the correct model for turbulence is the binomial multifractal cascade Lapidus . We ascribe our results as an imprint of intermittent cascades from discrete hierarchical dissipation in turbulence TurbDSI .

###### Acknowledgements.

The experimental turbulence data obtained at ONERA Modane were kindly provided by Y. Gagne. We are grateful to J. Delour and J.-F. Muzy for help in pre-processing these data and to A. Helmstetter for pointing out the correlation approach. This work was partially supported by NSF-DMR99-71475 and the James S. Mc Donnell Foundation 21st century scientist award/studying complex system.