text:
# STRUCTURE OF THE YANG-MILLS VACUUM IN THE ZERO MODES ENHANCEMENT QUANTUM MODEL

## I Introduction

The modern theory of strong interactions, quantum chromodynamics (QCD)  is intrinsically ultraviolet (UV) divergent theory because of asymptotic freedom  which is due to self-interaction of massless gluons only (non-Abelian character of QCD). It implies that any exact solution to QCD should become asymptotically free in the deep UV limit. This means that the perturbative UV divergences will be inevitably encountered in order to calculate any physical observable in this theory. Thus the existence of the UV divergences in QCD does not depend on the method of solution because behind these divergences is a physical phenomenon - asymptotic freedom. In other words, there is no way to escape from the UV divergences in QCD by implementing some sophisticated (other than perturbation theory) method of solution. Since divergences are physical (not artificial due to a method of solution), a mass scale parameter appears as a result of the UV renormalization program in order to render the theory finite in the UV limit. In QCD quark masses are not observable, so the UV divergences are to be absorbed into the scale parameter which precisely becomes finite in the UV limit. Instead of the dimensionless fixed coupling constant (like in QED, see a brief discussion below) we are left with the above-mentioned dimensional asymptotic scale parameter, $`\mathrm{\Lambda }_{QCD}`$ and therefore dimensionless coupling eventually becomes running. Numerically it is a few hundred $`MeV`$ only, so it can not survive in the UV (perturbative) limit, i.e., it can not be determined by the theory in the UV region. It should come from the infrared (IR) nonperturbative region. This clearly shows that this scale is completely nonperturbative by origin and this conclusion is strongly supported by the phenomenon of ”dimensional transmutation” .

Let us emphasize that the existence of the asymptotic scale parameter, $`\mathrm{\Lambda }_{QCD}`$, means that the full gluon propagator depends in general on this scale parameter, i.e., its behavior in the deep UV limit deviates from the free perturbative one (the above-mentioned running eefective coupling constant). In other words, the perturbative dynamics even in massless QCD is nontrivial (scale violation, asymptotic freedom). Moreover, let us note in advance, that due to the origin of the QCD scale parameter from the IR region, the asymptotics of the full gluon propagator in the deep IR limit should deviate from the free one as well. In other words, in the IR limit the QCD coupling is also running with its, generally speaking, own scale parameter (see below). In this connection, let us note that if QCD itself is confining theory, then some characteristic scale (different, in general, from asymptotic scale parameter, $`\mathrm{\Lambda }_{QCD}`$) is very likely to exist.

All this is in strict contradiction with quantum electrodynamics (QED), the theory of electromagnetic interactions. In QED the UV divergences are result of the method of solution, the perturbation theory. Any exact solution to QED, in principle, should not be asymptotically free. That is why neither the perturbative UV divergences nor the dimensional scale parameter necessarily should appear in QED. So the UV divergences in QED are not physical, i.e., there is no some nontrivial physics behind these divergences (unlike QCD, where asymptotic freedom is behind the UV divergences). However, mathematically the UV renormalization program (in order to render the theory finite) is, of course, the same as in QCD***Though historically the renormalization of QCD was done following the renormalization of QED. but unlike QCD, all the UV divergences in QED can be absorbed into the electron charge and mass terms. So there is no running coupling constant and the photon propagator remains free one in the deep UV limit (no deviation from the perturbative free behavior in this limit). There is no doubt, that the drastic difference between these two theories is due to Abelian character of QED (no direct interaction between massless photons).

If the structure of QCD in the UV perturbative region is well understood (because of asymptotic freedom, as it was briefly discussed above), little is known about its structure in the nonperturbative IR region. The main reason of this is, of course, the phenomenon of confinement which remains the most important unsolved yet problem in QCD. It is still remains to understand whether confinement (in particular quark confinement) is some additional external requirement in order to complete QCD to be a self-consistent theory of strong interactions or QCD itself is able to explain confinement and other nonperturbative phenomena without involving some extra degrees of freedom. At the very beginning of QCD, it was expressed an idea , that the quantum excitations of the IR degrees of freedom because of self-interaction of massless gluons only made it possible to understand confinement as well as other nonperturbative effects, for example such as dynamical chiral symmetry breakdown (DCSB). In other words, the importance of the IR structure of the true QCD vacuum was emphasized as well as its relevance to quark confinement, DCSB and other way around.

In accordance with these ideas, let us suppose that QCD is intrinsically IR divergent theory as well. Again the physical reason of these divergences is self-interaction of massless gluons only. Asymptotic freedom is behind the UV perturbative divergences, while the nonperturbative IR divergences themselves are behind confinement and other nonperturbative effects in QCD. Thus in our picture they are more fundamental than the perturbative ones. How to take them into account? Of course, by the corresponding behavior of the full gluon propagator in the deep IR domain since there is no hope for an exact solution(s). In this connection, let us remind the reader that the full dynamical information of any quantum gauge field theory such as QCD is contained in the corresponding quantum equations of motion, the so-called Schwinger-Dyson (SD) equations for lower (propagators) and higher (vertices and kernels) Green’s functions. These equations should be also complemented by the corresponding Slavnov-Taylor (ST) identities which in general relate the above mentioned lower and higher Green’s functions to each other . These identities are consequences of the exact gauge invariance and therefore $`\mathrm{"}areexactconstraintsonanysolutiontoQCD\mathrm{"}`$ . Precisely this system of equations can serve as an adequate and effective tool for the nonperturbative approach to QCD . Among the above-mentioned Green’s functions, the two-point Green’s function describing the full gluon propagator

$$iD_{\mu \nu }(q)=\left\{T_{\mu \nu }(q)d(q^2,\xi )+\xi L_{\mu \nu }(q)\right\}\frac{1}{q^2},$$
(1)
has a central place \[1,7-11\]. Here $`\xi `$ is a gauge fixing parameter ($`\xi =0`$, Landau gauge) and $`T_{\mu \nu }(q)=g_{\mu \nu }q_\mu q_\nu /q^2=g_{\mu \nu }L_{\mu \nu }(q)`$. Evidently, its free perturbative (tree level) counterpart is obtained by simply setting the full gluon form factor $`d(q^2,\xi )=1`$ in Eq. (1.1). In particular, the solutions to the above-mentioned SD equation for the full gluon propagator (1.1), are supposed to reflect the quantum structure of the QCD ground state. It is a highly nonlinear integral equation containing many different propagators, vertices and kernels \[1,7-11\]. For this reason it may have many (several) different exact solutions with different asymptotics in the deep IR limit (the UV limit because of asymptotic freedom is uniquely determined), describing thus many different types of quantum excitations of gluon field configurations in the QCD vacuum. Evidently, not all of them can reflect the real structure of the QCD vacuum. Let us emphsize that we even do not know the complete set of boundary conditions in order to attempt to uniquely fix theory.

The deep IR asymptotics of the full gluon propagator can be generally classified into the three different types: 1) the IR enhanced (IRE) or IR singular (IRS) (for references see the next section below), 2) the IR finite (IRF)  and 3) the IR vanishing (IRV) ones . As it was emphasized above, any deviation in the behavior of the full gluon propagator in the IR domain from the free one automatically assumes its dependence on a scale parameter (at least one) responsible for the nonperturbative dynamics, say, $`\mathrm{\Lambda }_{NP}`$. This is very similar to asymtotic freedom which requires asymptotic scale parameter, $`\mathrm{\Lambda }_{QCD}`$ associated with nontrivial perturbative dynamics (scale violation). Thus, any realistic Ansatz (since exact solution(s) is(are) not known) for the full gluon propagator in QCD will contain at least the two above-mentioned scale parameters. In our approach, the nonperturbative scale is more fundamental than the asymptotic one, and if the above-mentioned Ansatz is realistic one, than, in principle, it should give rise to the relation between them, namely the asymptotic scale as $`purenumber`$ times the nonperturbative scale. The main purpose of this article is precisely to establish such relation within our approach to nonperturbative QCD and its ground state.

However, a few general remarks are in order. The nonperturbative QCD vacuum is a very complicated medium and its dynamical and topological complexity \[1,14-17\] means that its structure can be organized at various levels (classical and quantum) and it can contain many different components and ingredients which contribute to the truly nonperturbative vacuum energy density , one of the main characteristics of the QCD ground state. The classical component of the vacuum energy density is determined by the topologically nontrivial, instanton-type fluctuations of gluon fields configurations in the QCD vacuum which are nonperturbative weak-coupling limit solutions to the classical equation of motion in Euclidean space. Numerically this contribution is estimated by the random and interacting instanton liquid models (RILM and IILM) of the QCD vacuum (see Ref.  and references therein).

Here we are going to discuss the quantum part of the vacuum energy density which is determined by the effective potential approach for composite operators . It allows us to investigate the nonperturbative QCD vacuum, since in the absence of external sources the effective potential is nothing but the vacuum energy density. It gives the vacuum energy density in the form of the loop expansion where the number of the vacuum loops (consisting of the confining quarks and nonperturbative gluons properly regularized with the help of ghosts) is equal to the power of the Plank constant, $`\mathrm{}`$.

It is well known, however, that the vacuum energy density is generally badly divergent in quantum field theory, in particular QCD (see, for example, the discussion given by Shifman in Ref. ). The main problem thus was how to extract the truly nonperturbative vacuum energy density which should be finite, automatically negative without imaginary part (stable vacuum). This problem was addressed and resolved in our recent publications in Ref.  for covariant gauge QCD and in Ref.  for axial gauge QCD. Though the truly nonperturbative vacuum energy density is important in its own right (it is the bag constant, apart from the sign, by definition and precisely it is related to the gluon and quark condensates, topological susceptibility, etc. (see Refs. \[23-26\] and references therein)), here we are going to investigate the general structure of the Yang-Mills (YM) vacuum within our approach.

In the next section II we will formulate our main dynamical assumption which lies in the heart of our approach to nonperturbative QCD in general and to its ground state in particular. Then in section III we will apply to the above-mentioned approach for composite operators in order to analyze the vacuum structure. In section IV we present our main result while section V is devoted to discussion. Our conclusions are given in section VI.

## II The IRE gluon propagator. ZME quantum model

Today there are no doubts left that the dynamical mechanisms of the important non-perturbative quantum phenomena such as quark confinement and dynamical (or equivalently spontaneous) chiral symmetry breaking (DCSB) are closely related to the complicated topologically nontrivial structure of the QCD vacuum \[1,14-17\]. On the other hand, it also becomes clear that the nonperturbative IR dynamical singularities, closely related to the nontrivial vacuum structure, play an important role in the large distance behaviour of QCD . For this reason, any correct nonperturbative model of quark confinement and DCSB necessarily turns out to be a model of the true QCD vacuum and the other way around.

Our model of the true QCD ground state is based on the existence and importance of such kind of the nonperturbative, quantum excitations of the gluon field configurations (due to self-interaction of massless gluons only, i.e., without explicit involving some extra degrees of freedom) which can be effectively correctly described by the $`q^4`$ behavior of the full gluon propagator in the deep IR domain (at small $`q^2`$) . Thus our main Ansatz for the full gluon form factor in Eq. (1.1) in the whole momentum range $`[0,\mathrm{})`$ is

$$d(q^2,\xi )=\frac{\mu ^2}{(q^2)}+\frac{c}{\mathrm{ln}(q^2/\mathrm{\Lambda }_{YM}^2)}.$$
(2)
Evidently, the UV piece is nothing but the running coupling constant $`\overline{g}^2(q^2)`$, which is the leading order solution to the renormalization group equations  and $`c=(16\pi ^2/11)`$ is the inverse of the first coefficient of the $`\beta `$-function in the case of YM gluon fields (pure gluodynamics). In the above-mentioned Ref.  it is denoted as $`2/b_0`$. The flavorless QCD asymptotic scale parameter is denoted as $`\mathrm{\Lambda }_{YM}`$. In the weak coupling limit ($`q^2\mathrm{}`$) precisely the second term becomes dominant (asymptotic freedom).

On the other hand, we obtain the generally accepted form of the IR singular asymptotics for the full gluon propagator (for some references see below)

$$D_{\mu \nu }(q)(q^2)^2,q^20,$$
(3)
which may be equivalently referred to as the strong coupling regime . It describes zero momentum modes enhancement (ZMME) dynamical effect in QCD at large distances. In our previous publications and in future we prefer to use simply ZME (zero modes enhancement) since we work in momentum space always. This is our primary dynamical assumption in this section. The main problem with this strong singularity is its correct treatment by the dimensional regularization method  within the distribution theory , which precisely was one of highlights of our previous publications  (see also Refs. ). There exist many arguments in favor of this behaviour. Let us remind the reader a few of them which are the most important in our opinion.

a). Precisely this strongly singular behaviour of the full gluon propagator in the IR domain leads to the area law for static quarks (indicative of confinement) within the Wilson loop approach .

b). The cluster property of the Wightman functions in QCD fails and this allows such singular behavior like (2.2) for the full gluon propagator in the deep IR domain .

c). After the pioneering papers of Mandelstam in the covariant (Landau) gauge  and Baker, Ball and Zachariasen in the axial gauge , the consistency of the singular asymptotics (2.2) with direct solution to the SD equation for the full gluon propagator in the IR domain was repeatedly confirmed (see for example Refs.  and references therein).

d). Moreover, let us underline that without this component in the decomposition of the full gluon propagator in continuum theory it is impossible to ”see” linearly rising potential between heavy quarks by lattice QCD simulations  not involving some extra (besides gluons and quarks) degrees of freedom. This should be considered as a strong lattice evidence (though not direct) of the existence and importance of $`q^4`$-type excitations of gluon field configurations in the QCD vacuum. There also exists direct lattice evidence that the zero modes are enhanced in the full gluon propagator indeed .

e). Within the distribution theory  the structure of the nonperturbative IR singularities in four-dimensional Euclidean QCD is the same as in two-dimensional QCD, which confines quarks at least in the large $`N_c`$ limit . In this connection, let us note that $`q^4`$ IR singularity is the simplest nonperturbative power singularity in four-dimensional QCD as well as $`q^2`$ IR singularity is the simplest nonperturbative power singularity in two-dimensional QCD. The QCD vacuum is much more complicated medium than its two-dimensional model, nevertheless, the above-mentioned analogy is promising even in the case of the nonperturbative dynamics of light quarks.

f). Some classical models of the QCD vacuum also invoke $`q^4`$ behavior of the gluon fields in the IR domain. For example, it appears in the QCD vacuum as a condensation of the color-magnetic monopoles (QCD vacuum is a chromomagnetic superconductor) proposed by Nambu, Mandelstam and ’t Hooft and developed by Nair and Rosenzweig (see Ref.  and references therein. For recent developments in this model see Di Giacomo’s contribution in Ref. ) as well as in the classical mechanism of the confining medium  and in effective theory for the QCD vacuum proposed in Ref. .

g). It is also required to derive the heavy quark potential within the recently proposed exact renormalization group flow equations approach .

h). It has been shown in our papers that the singular behavior (2.2) is directly related to light quarks confinement and DCSB .

i). These excitations are also topologically nontrivial since they lead to the nontrivial (truly nonperturbative) vacuum energy density which is finite, automatically negative and it has no imaginary part (stable vacuum) .

j). Moreover, they nicely saturate the phenomenological values of the topological susceptibility, the mass of the $`\eta ^{}`$ meson in the chiral limit and the gluon condensate .

Thus we consider our main Ansatz (2.1) as physically well-motivated. The first term in the decomposition (2.1) is the truly nonperturbative part since it vanishes in the perturbative limit ($`\mu ^20`$), when the perturbative phase survives only and simultaneously it correctly reproduces the deep IR asymptotics of the full gluon propagator.

## III The effective potential for composite operators

In this section we are going to analytically investigate the quantum structure of the YM vacuum which is due to $`q^4`$-type nonperturbative, topologically nontrivial excitations of the gluon field configurations there complemented by the renormalization group improvements for the UV region, (2.1). Let us start from the gluon part of the vacuum energy density which to leading order (log-loop level $`\mathrm{}`$)Next-to-leading and higher terms (two and more vacuum loops) are suppressed by one order of magnitude in powers of $`\mathrm{}`$ at least and are left for consideration elsewhere. is given by the effective potential for composite operators  as follows

$$V(D)=\frac{i}{2}\frac{d^nq}{(2\pi )^n}Tr\{\mathrm{ln}(D_0^1D)(D_0^1D)+1\},$$
(4)
where $`D(q)`$ is the full gluon propagator (1.1) and $`D_0(q)`$ is its free perturbative (tree level) counterpart. Here and below the traces over space-time and color group indices are understood. The effective potential is normalized as $`V(D_0)=0`$, i.e., free perturbative vacuum is normalized to zero. In order to evaluate the effective potential (3.1) we use the well-known expression,

$$Tr\mathrm{ln}(D_0^1D)=8\times \mathrm{ln}det(D_0^1D)=8\times 4\mathrm{ln}\left[\frac{3}{4}d(q^2)+\frac{1}{4}\right].$$
(5)
It becomes zero (in accordance with the above mentioned normalization condition) when the full gluon form factor is replaced by its free counterpart. The composition (3.2) does not explicitly depend on a gauge choice.

Going over to four ($`n=4`$) dimensional Euclidean space in (3.1-3.2), on account of (2.1), by introducing the appropriate dimensionless momentum variable as $`z=q^2/\mu ^2`$ and after doing some algebra, we finally obtain the following expression for the vacuum energy density, $`ϵ_g=V(D)`$,

$$ϵ_g=\frac{1}{\pi ^2}q_0^4z_0^2I_g(z_0;b),$$
(6)
where

$$I_g(z_0;b)=\underset{0}{\overset{z_0}{}}𝑑zz\left\{\mathrm{ln}z+\frac{3}{4z}a+\frac{3c}{4\mathrm{ln}bz}\mathrm{ln}[3+z+\frac{3cz}{\mathrm{ln}bz}]\right\}=I_g^r(z_0;b)+I_g^s(z_0;b),$$
(7)
and

$`I_g^r(z_0;b)`$ $`=`$ $`{\displaystyle \frac{1}{2}}z_0^2\mathrm{ln}\overline{b}z_0{\displaystyle \frac{1}{4}}z_0^2+{\displaystyle \frac{3}{4}}z_0I(z_0;b)`$ (8)
$`I_g^s(z_0;b)`$ $`=`$ $`{\displaystyle \frac{\overline{a}}{2}}z_0^2+{\displaystyle \frac{3c}{4b^2}}li(bz_0)^2+{\displaystyle \frac{1}{2}}z_0^2\mathrm{ln}\mathrm{ln}bz_0+{\displaystyle \frac{1}{2b^2}}Ei(1,2\mathrm{ln}bz_0),`$ (9)
with $`\overline{a}=a+\mathrm{ln}\overline{b}`$ and $`a=(3/4)2\mathrm{ln}2`$ while $`\overline{b}=\mathrm{ln}bz_0+3c`$. In these integrals $`z_0=q_0^2/\mu ^2`$ while $`b=\mu ^2/\mathrm{\Lambda }_{YM}^2`$ and $`li(x)`$ is the logarithm-integral function. Let us note in advance, that here and everywhere below the combination $`bz_0=q_0^2/\mathrm{\Lambda }_{YM}^2=x_0`$ is fixed (constant) when $`q_0`$ itself is fixed, i.e., before going to infinity at final stage. Thus, the parameter $`x_0`$ is the truly dimensionless UV cutoff since when the UV cutoff in momentum space $`q_0^2`$ goes to infinity ($`q_0^2\mathrm{}`$) then it certainly goes to infinity as well. At the same time, the variable $`z_0`$ only looks like as the dimensionless UV cutoff. In the limit $`q_0^2\mathrm{}`$ it might be finite since the nonperturbative scale, $`\mu ^2`$, can be affected by the UV divergences as well (”bare” parameter). Precisely this takes place within our approach in this case (see below). In fact, variable $`z_0`$ determines the perturbative limit ($`z_0\mathrm{}`$) at fixed $`q_0^2`$, when the perturbative phase survives only in our model ($`\mu 0`$). $`Ei(1,2\mathrm{ln}bz_0)`$ is the exponential-integral function with unphysical singularity at $`bz_0=1`$ as well as the preceding term has. In the compositions shown in (3.5), the corresponding integral is defined as follows

$$I(z_0;b)=\underset{0}{\overset{z_0}{}}𝑑zz\mathrm{ln}[(3+z)\mathrm{ln}bz+3cz].$$
(10)
The effective potential at the log-loop level for the ghost degrees of freedom is

$$V(G)=i\frac{d^np}{(2\pi )^n}Tr\{\mathrm{ln}(G_0^1G)(G_0^1G)+1\},$$
(11)
where $`G(p)`$ is the full ghost propagator and $`G_0(p)`$ is its free perturbative counterpart. The effective potential $`V(G)`$ is normalized as $`V(G_0)=0`$. Evaluating the ghost term $`ϵ_{gh}=V(G)`$ in (3.7), we obtain $`ϵ_{gh}=\pi ^2k_0^4y_0^2I_{gh}(y_0;b)`$, where the integral $`I_{gh}(y_0;b)`$ depends on the ghost propagator, which remains arbitrary (unknown) within our approach and $`y_0=k_0^2/\mathrm{\Lambda }_{NP}^2`$.

In principle, we must sum up all contributions in order to obtain total vacuum energy density (the confining quark part of the vacuum energy density is not considered here). However, approximating the full gluon propagator by the sum of its deep IR (confinement) and UV asymptotics in (2.1), the nonzero terms appear in (3.4) which precisely violates the above mentioned normalization condition of the free perturbative vacuum to be zero when the parameter $`b`$ goes to zero ($`\mu 0,z_0\mathrm{}`$). In this case the perturbative phase only remains within our model. These terms are collected in the expression (3.5) for $`I_g^s(z_0;b)`$. It also contains two terms which suffer from unphysical singularity \- Landau pole - at $`x_0=1`$. This should be traced back to (2.1) when $`q^2=\mathrm{\Lambda }_{YM}^2`$ since we use it in the whole range. These terms should be subtracted from the gluon part of the vacuum energy density. Using the arbitrariness of the ghost term (3.7), this can be done by imposing the following condition, $`\mathrm{\Delta }=I_g^s(z_0;b)+I_{gh}(y_0;b)=0`$ (recalling that $`\mathrm{\Lambda }_{YM}^4=k_0^4y_0^2=q_0^4z_0^2`$). We regularize the gluon contribution to the vacuum energy density by subtracting unwanted terms by means of the ghost contribution, i.e., it plays the role of the corresponding counter term. So we define, $`ϵ_g+ϵ_{gh}=ϵ_g^{reg}=\pi ^2q_0^4z_0^2I_g^r(z_0;b)`$, on account of $`\mathrm{\Delta }=0`$. Thus our regularization procedure is in agreement with general physical interpretation of ghosts to cancel the effects of unphysical degrees of freedom of the gauge bosons .

Let us introduce the effective potential $`\overline{\mathrm{\Omega }}_g`$ at fixed $`q_0`$ as follows

$$\overline{\mathrm{\Omega }}_g\overline{\mathrm{\Omega }}_g(z_0;b)=\frac{1}{q_0^4}ϵ_g^{reg}=\frac{1}{2\pi ^2}\times z_0^2\left\{z_0^2\mathrm{ln}\overline{b}z_0\frac{1}{2}z_0^2+\frac{3}{2}z_02I(z_0;b)\right\},$$
(12)
where the integral $`I(z_0;b)`$ is given in (3.6). First of all, it is always useful to factorize the scale dependence in the regularized and subtracted vacuum energy density by means of the above-introduced effective potential at fixed scale. On the other hand, it makes it possible to investigate the structure of the vacuum within our approach since this quantity is finite when the UV cutoff goes to infinity (see below) as well as it is free from unphysical singularities (see above). These two remarks emphasize the importance of the effective potential (3.8). The asymptotic behavior of the effective potential (3.8) depends on the asymptotic properties of the integral $`I(z_0;b)`$, (3.6). It is almost obvious that its asymptotics at $`z_00,\mathrm{}`$ (and consequently at $`b0,\mathrm{}`$) to-leading order can be easily evaluated analytically by replacement $`\mathrm{ln}bz`$ by $`\mathrm{ln}bz_0`$ since logarithm is a very slowly varing function. Omitting all intermediate analytical calculations, one obtains

$$I(z_0;b)=\frac{1}{2}z_0^2\mathrm{ln}(\overline{b}z_0)\frac{1}{4}z_0^2+\frac{3}{2\overline{c}}z_0\frac{9}{2\overline{c}^2}\mathrm{ln}\left\{1+\frac{\overline{c}z_0}{3}\right\}+\frac{1}{2}z_0^2\mathrm{ln}\left\{1+\frac{3}{\overline{c}z_0}\right\},$$
(13)
where $`\overline{c}=1+(3c/\mathrm{ln}bz_0)`$. Thus in the perturbative limit $`z_0\mathrm{}`$ one has

$$\overline{\mathrm{\Omega }}_g(z_0;b)_{z_0\mathrm{}}\frac{9}{4\pi ^2}\frac{\mathrm{ln}bz_0c}{\mathrm{ln}bz_0+3c}z_0^1,$$
(14)
which clearly shows that in order for the effective potential (3.8) as a function of $`z_0`$ to approach zero from below in the perturbative limit, one needs to set

$$\mathrm{ln}x_0=\lambda c,\lambda >1,$$
(15)
recalling that $`bz_0=x_0`$ by definition. The expression (3.11) is nothing but useful parametrization of the UV divergences in terms of the parameter $`\lambda `$, which will go to infinity at the final stage. At the same time, at zero ($`z_00`$) the effective potential (3.8) behaves as $`z_0^1`$. Thus as a function of $`z_0`$ it has a minimum if the condition (3.11) takes place which is the case indeed ($`\lambda `$ goes to infinity at final stage).

The minimization of the effective potential (3.8) with respect to $`z_0`$ (at fixed $`bz_0`$), $`\overline{\mathrm{\Omega }}_g(z_0;b)/z_0=0`$, yields

$$z_0^2+4I_2(z_0;b)=\frac{3}{2}z_0+2z_0^2\mathrm{ln}[(3+z_0)\mathrm{ln}bz_0+3cz_0].$$
(16)
It is easy to show that to-leading order in the $`\lambda \mathrm{}`$ limit, it can be analytically evaluated as follows:

$$(\lambda 1)(\lambda +3)z_0=4\lambda ^2\mathrm{ln}\left\{1+\frac{(\lambda +3)z_0}{3\lambda }\right\}.$$
(17)
In the derivation of this expression we have already used expression for $`\overline{c}`$ (see text after Eq. (3.9)) and (3.11). Its solution in the $`\lambda \mathrm{}`$ limit is

$$z_0^{min}(\lambda \mathrm{})=\left\{\frac{x_0}{b}\right\}_{min}(\lambda \mathrm{})=2.2,$$
(18)
i. e. the point $`z_0^{min}=2.2`$ is the point of accumulation of minimums in the UV limit of the effective potential (3.8). This effect shows that the parameter $`b`$ goes to infinity in the same way as the UV cut-off $`x_0`$ in order for their ratio to go to the above displayed finite limit. This means that initial mass scale parameter $`\mu `$ is a ”bar” one. It depends on the UV cut-off $`\lambda `$, i.e., $`\mu =\mu (\lambda )`$ and consequently $`b=b(\lambda )`$Let us note that it does not make any sense to introduce finite (from the very beginning) the nonperturbative mass scale parameter. The problem is that in any case it will be affected by the UV divergences since we have explicitly kept the UV term in Eq. (2.1) (did not subtract it). At the same time, the asymptotic scale parameter, $`\mathrm{\Lambda }_{YM}`$, is finite (from the very beginning) by definition.. In what follows we will introduce the UV renormalized $`\mu `$ and consequently the UV renormalized $`b`$. To-leading order in the $`\lambda \mathrm{}`$ limit the effective potential (3.8), on account of the ”stationary” condition (3.13), becomes

$$\overline{\mathrm{\Omega }}_g(\lambda ,z_0)=\frac{1}{2\pi ^2}\left\{\frac{3}{4}z_0^1\mathrm{ln}(1+\frac{3\lambda }{z_0(\lambda +3)})\right\}.$$
(19)
## IV The YM vacuum structure

Evaluating the effective potential $`\mathrm{\Omega }_g(\lambda ,z_0^{min}(\lambda ))`$ (3.15) numerically in the vicinity of the point of accumulation of the minimums, $`z_0^{min}(\lambda \mathrm{})2.2`$ (see (3.14)), one obtains $`\mathrm{\Omega }_g(\lambda \mathrm{})=0.0263`$, i.e., it is finite in this limit. Recalling that $`ϵ_g^{reg}=ϵ_g^{reg}(\lambda )`$, then from Eq.(3.8) it follows

$$ϵ_g^{reg}(\lambda )=0.0263q_0^4(\lambda ).$$
(20)
On the other hand,

$$q_0^4(\lambda )=x_0^2(\lambda )\mathrm{\Lambda }_{YM}^4=e^{2\lambda c}\mathrm{\Lambda }_{YM}^4,$$
(21)
by definition and because of (3.11). Thus the vacuum energy becomes

$$ϵ_g^{reg}(\lambda )=0.0263\mathrm{\Lambda }_{YM}^4\times e^{2\lambda c}$$
(22)
and it is badly divergent in the UV limit ($`\lambda \mathrm{}`$) because of the exponential $`\lambda `$-factor. This result should be expected since we have explicitly kept (did not subtract) the UV peace in Eq. (2.1) (for general discussion see Ref. ).

However, despite badly divergent vacuum energy density, our method makes it possible to establish a finite relation between two scale parameters of our model. Indeed, as was mentioned above, the parameter $`b`$ diverges exactly as the UV cutoff $`x_0=e^{\lambda c}`$ (see Eq. (3.11)) itself in order to get final result in the $`\lambda \mathrm{}`$ limit. So let us introduce now the UV finite $`b`$ and $`\mu `$ as

$$\stackrel{~}{b}=b(\lambda )e^{\lambda c}$$
(23)
and

$$\mathrm{\Lambda }_{NP}=\mu (\lambda )e^{\lambda c/2},$$
(24)
respectively. Then from Eq. (3.14) one obtains

$$z_0^{min}=\frac{x_0}{b(\lambda )}=\frac{e^{\lambda c}}{\stackrel{~}{b}e^{\lambda c}}=\stackrel{~}{b}^1=2.2$$
(25)
in the $`\lambda \mathrm{}`$ limit. On the other hand, by definition

$$\mathrm{\Lambda }_{YM}^2=\frac{\mu ^2(\lambda )}{b(\lambda )}=\frac{\mathrm{\Lambda }_{NP}^2e^{\lambda c}}{\stackrel{~}{b}e^{\lambda c}}=\frac{\mathrm{\Lambda }_{NP}^2}{\stackrel{~}{b}}$$
(26)
in the same limit. This finally leads to

$$\mathrm{\Lambda }_{YM}=1.48324\times \mathrm{\Lambda }_{NP}.$$
(27)
A few remarks are in order. Let us underline that the above-displayed numerical result (4.8) was obtained in a gauge invariant way, i.e., this number did not explicitly depend on a gauge choice. Also this result will not be changed in the next-to-leading order approximation to the running coupling constant in the UV limit since they were obtained to-leading dominant order limit. Evidently, it should not, of course, depend on how one introduces the UV cutoff by making use of different scale parameters. For example, it can be defined as $`x_0=q_0^2/\mathrm{\Lambda }_{YM}^2`$ (from the very beginning) while the parameter $`b`$ remains the same. Then the minimization of the effective potential $`\mathrm{\Omega }_g(x_0;b)=ϵ_g^{reg}(x_0;b)\mathrm{\Lambda }_{YM}^4`$ with respect to $`b`$ finally yields the same numerical relation between the nonperturbative and asymptotic scales, (4.8) as it should be in principle (for details see Ref. ). It is worth noting that the minimization with respect to $`x_0`$ itself in this case leads to trivial zero (see Ref. ).

## V Discussion

The main problem now is the numerical value of the YM asymptotic scale parameter, $`\mathrm{\Lambda }_{YM}\mathrm{\Lambda }`$ (for brevity and further aims), since it can not be experimentally determined in the YM theory alone. So its numerical value is to be taken from lattice approach or phenomenology. Unfortunately, lattice estimates of the QCD flavorless asymtotic scale parameter strongly depend on the renormalization scheme chosen for calculation. In Ref.  in the $`\overline{MS}`$ scheme it is expressed in terms of Sommer’s scale $`r_0`$, namely $`\mathrm{\Lambda }_{\overline{MS}}=0.636(54)/r_0`$. Than Eq. (4.8) becomes $`\mathrm{\Lambda }_{NP}=0.4288r_0^4`$. In physical units, using $`r_0=0.5fm`$, one numerically gets, $`\mathrm{\Lambda }_{NP}=171.5MeV`$. However, as it was emphasized in the above mentioned Ref. , the conversion to physical units is rather ambiguous in the pure gauge theory. Approximately the same numerical result for $`\mathrm{\Lambda }_{\overline{MS}}`$ in terms of the string tension was previously obtained in Ref. . In Ref.  the YM scale parameter was calculated first in the $`\stackrel{~}{MOM}`$ scheme and then it was matched to the $`\overline{MS}`$ scheme. The corresponding numerical results (here and everywhere below up to central values, for simplicity) are: $`\mathrm{\Lambda }_{\stackrel{~}{MOM}}=0.88GeV`$ and $`\mathrm{\Lambda }_{\overline{MS}}=0.31GeV`$. Substituting these values into the Eq. (4.8), one arrives at the following estimates, $`\mathrm{\Lambda }_{NP}^{\stackrel{~}{MOM}}=593.3MeV`$ and $`\mathrm{\Lambda }_{NP}^{\overline{MS}}=210MeV`$. The most recent lattice calculation of the flavorless QCD asymptotic scale parameter is being presented in Ref.  are: $`\mathrm{\Lambda }_{MOM}=0.361(6)GeV`$ and $`\mathrm{\Lambda }_{\stackrel{~}{MOM}}=0.349(6)GeV`$ while $`\mathrm{\Lambda }_{\overline{MS}}=0.412GeV`$. Substituting this result into the Eq. (4.8) for example for the last value, one gets, $`\mathrm{\Lambda }_{NP}^{\overline{MS}}=277.7MeV`$.

In Ref.  the vacuum energy density of the covariant model of Euclidean QCD with $`N_f=3`$ quark flavors on a finite $`D`$-dimensional symmetric hypertorus was calculated. Adjusting the recently calculated value of the gluon condensate by Narison  and the numerical value of $`\alpha _s`$ from $`\tau `$ decay, in the above-mentioned paper  it was estimated that $`\mathrm{\Lambda }_{\overline{MS}}^{(3)}500600MeV`$ (compatible with the standard value of the gluon condensate ). From the renormalization group equations, however, it is well-known that the numerical value of the flavorless QCD asymptotic scale parameter should be bigger than its full QCD counterparts, for example, the above estimated $`\mathrm{\Lambda }_{\overline{MS}}^{(3)}`$. Thus roughly estimating $`\mathrm{\Lambda }_{YM}600700MeV`$ in Eq. (4.8), one obtains $`\mathrm{\Lambda }_{NP}404.5472MeV`$. Contrary to the above-discussed cases, in Refs.  the nonperturbative scale parameter was calculated itself in chiral QCD within our approach. Because of the scale-setting scheme in order to relate it to a good physical observeble from full QCD , its numerical value depends on a bounds for the pion decay constant in the chiral limit. The corresponding bounds for the nonperturbative scale parameter, we have obtained, were: $`540.6\overline{\mu }578.4MeV`$. Then from Eq. (4.8) it immediately follows $`801.8\mathrm{\Lambda }_{YM}858MeV`$ by identifying $`\overline{\mu }`$ with $`\mathrm{\Lambda }_{NP}`$ for simplicity. Here we are observing the standard situation when estimates from phenomenology (continuum theory) are usually bigger than lattice estimates of the same scale quantities.

In any case, our solid result is the relation (4.8) which is to be better understood as the ratio between two scale parameters from lattice simulation point of view. In this way, in principle, it should be free from lattice artifacts since it does not require the conversion to physical unuts which still remains rather ambiguous in the pure gauge theory as it was mentioned above.

## VI Conclusions

In summary, we have formulated new quantum model of the QCD vacuum using the effective potential approach for composite operators. It is based on the existence and importance of such kind of the nonperturbative, topologically nontrivial excitations of gluon field configurations (due to self-interaction of massless gluons only, i.e., without explicit involving some extra degrees of freedom) which can be effectively correctly described by the $`q^4`$-type behaviour of the full gluon propagator in the deep IR domain. The UV part of the full gluon propagator was approximated by asymptotic freedom to-leading order log term of the running coupling constant. By explicit retaining the UV tail in the decomposition of the full gluon propagator (2.1), the vacuum energy density remains badly divergent, of course (for general discussion see again Refs.  and Shuryaks’s review in Ref. ). However, despite this, we have formulated a method to establish a finite (in the UV limit) relation (4.8) between two scale parameters of our model. We propose to minimize the properly regularized and subtracted effective potential at a fixed scale (3.8) as a function of the variable $`z_0`$. It has a clear physical meaning when it goes to infinity then the perturbative phase survives only within our approach to nonperturbative QCD. This is possible to do because the above-mentioned effective potential is finite in the UV limit within our model, i.e., this is a specific property of our Ansatz for the full gluon propagator (2.1). Equivalently, one can minimize the auxiliar (again properly regularized and subtracted) effective potential as a function of the parameter which has a clear physical meaning. When it is zero then again the perturbative phase survives only.

For the first time, we have expressed the asymptotic scale parameter as $`purenumber`$ times the nonperturbative scale, namely (4.8). We don’t know any other YM quantum/classical vacuum model which possesses the same property. This can be considered as one more argument in favor of our model. In principle, in QCD one of the final goals is precisely to express all scale parameters (for example, hadron masses) as $`purenumbers`$ times some characteristic scale of the theory. Despite the perturbative vacuum cannot be the QCD true ground state , nevertheless the existence of such kind of relation is a manifestation that ”the problems encountered in perturbation theory are not mere mathematical artifacts but rather signify deep properties of the full theory” . In support of this statement, let us underline one more that the running coupling with scaling violations not only in the UV but in the IR as well is the most important feature of QCD which dratically differs this theory from QED. The former is certainly related to asymptotic freedom while the latter should be apparently related to confinement, in particular quark confinement. Behind both types of scaling violation (in general, with their own scale parameters) is the self-interaction of massless gluons only within our approach to nonperturbative QCD in general and to its ground state in particular. Thus it suggests a unified description of both phenomena, so such kind of relation between two scale parameters, in principle, should be expected. At the same time, it does not mean that they are dependent from each other. It means rather that there exists the procedure (method) how to establish connection between them. Precisely this point we wanted to emphasize in the Introduction.

###### Acknowledgements.

The author would like to thank M. Polikarpov, M. Chernodub, A. Ivanov and especially V.I. Zakharov for useful and critical discussions as well as A.T. Filippov for corespondence. It is a pleasure to thank Gy. Kluge for many remarks and help.This work was supported by HAS-JINR collaboration agreement.


format: text

text:
# Chandra X-ray Observatory Study of the Orion Nebula Cluster and BN/KL Region

## 1 Introduction

X-ray emission among pre-main sequence (PMS) stars with masses $`M2`$ M is known to be elevated $`10^1`$ to $`>10^4`$ times above typical main sequence levels (Feigelson & Montmerle, 1999). High-amplitude X-ray variability and hard spectra, supported by multiwavelength studies, indicate the emission is in most cases due to solar-type magnetic flares where plasma is heated to high temperatures by violent reconnection events in magnetic loops. However, the empirical relationships between PMS X-ray emission and stellar age, mass, radius, and rotation are complex (e.g., Feigelson et al., 1993) and their links to astrophysical theory are not strong. For example, if surface magnetic activity is a manifestation of a simple rotation-driven magnetic dynamo in the stellar interior, and if current theories about the evolution of angular momentum in PMS stars hold (Bouvier et al., 1997), then specific predictions about the evolution of the X-ray luminosity function during the PMS phase can be tested if sufficiently large samples were obtained. X-rays are known to be produced in high-mass star formation regions but little is known about the structure or astrophysical origins of this emission (Churchwell, 1999).

The largest sample of PMS stars observable in a single field is the Orion Nebula Cluster (ONC). The OB members of the ONC illuminate the Orion Nebula (= Messier 42), a blister HII region at the near edge of the Orion A giant molecular cloud. The ONC is a nearby and one of the densest PMS stellar clusters known with $`2000`$ members concentrated in a 1 pc ($`8^{}`$) radius sphere with 80% of the stars younger than 1 My (Hillenbrand, 1997). It was the first star forming region to be detected in the X-ray band (Giacconi et al., 1972) and non-imaging studies soon found the X-ray emission is extended on scales of a parsec or larger den Boggende et al. (1978); Bradt & Kelley (1979). Early explanations for the Orion X-rays included winds from the massive Trapezium stars colliding with each other or the molecular cloud, and hot coronae or magnetic activity in lower mass T Tauri stars. The Einstein (Ku & Chanan, 1979), ROSAT (Gagné et al., 1995; Geier et al., 1995; Alcalá et al., 1996) and ASCA (Yamauchi et al., 1996) imaging X-ray observatories established that both the Trapezium stars and many lower-mass T Tauri stars contribute to the X-ray emission.

These studies could detect only a modest fraction of the ONC stars due to the high stellar concentration and heavy absorption by molecular material along the line of sight. NASA’s Chandra X-ray Observatory launched on 23 July 1999 with its Advanced CCD Imaging Spectrometer (ACIS) detector provides an unprecedented combination of capabilities $``$ a wide spectral bandpass, $`<1`$″ spatial resolution, very low cosmic ray and detector background levels, and moderate spectral resolution at every pixel $``$ that make it particularly well-suited to the study of crowded and absorbed clusters of faint X-ray sources.

After describing the data in §2, we discuss demographics of PMS stellar X-ray emission (§3). The next sections describe two of many possible investigations of this rich dataset: in §4 we examine the evolution of magnetic activity in solar-mass stars, and in §5 we present the X-ray sources from stars in the Becklin-Neugebauer/Kleinman-Low (BN/KL) star forming region in OMC 1. Later papers will provide a more comprehensive analysis of this and additional observations of the Orion Nebula (Feigelson et al., in preparation).

## 2 Observations and data analysis

The ONC was observed with the ACIS-I array detector on board the Chandra X-ray Observatory Weisskopf et al. (1996) for 47.8 ks on 12 October 1999. The array consists of four abutted charge-coupled devices (CCDs) specially designed for X-ray studies (Garmire et al., in preparation111Detailed descriptions of the ACIS instrument and its operation can be found on-line at www.astro.psu.edu/xray/docs/sop and www.astro.psu.edu/xray/docs/cal\_report.). Two chips from the ACIS-S spectroscopic array were also turned on, but their data will not be discussed here. Starting with Level 1 processed event lists provided by the pipeline processing at the Chandra X-ray Center, we removed several types of CCD events likely dominated by particles rather than imaged X-rays: energies above 8 keV and below 0.2 keV; split events except for ASCA grades 0, 2, 3, 4 and 6; cosmic ray events causing stripes along amplifier nodes; hot columns and other low-quality flagged events. The aimpoint on chip I3 was $`5^h35^m15.0^s`$ $`5^{}23^{}20^{\prime \prime }`$.

The resulting image (Fig. 1) shows hundreds of sources concentrated around the bright Trapezium stars. Several instrumental effects are also evident such as the broadening of the point response function off-axis due to the telescope optics, the 11″ gaps between the CCD chips blurred by the satellite aspect dither, and the trailing of events from the brightest sources due to events arriving during the chip readout. The brightest Trapezium sources also suffer photon pileup which distorts their morphology, flux and spectrum.

An additional problem in the current dataset is the degradation of charge transfer efficiency of the ACIS-I CCDs due to high dosage of particle irradiation early in the Chandra mission (Prigozhin et al., 2000). Corrective measures were taken to arrest the degradation as soon as it was discovered and data acquired after 29 January 2000 should have a reduced problem due to changes in on-board operation of the detector to mitigate charge transfer inefficiency. Improved data processing techniques have been developed to treat this problem (Townsley et al., 2000) and are applied in our analysis in §5. The results of other portions of the study should not be seriously affected by the charge transfer degradation.

The excellent point spread function produced by the Chandra mirrors and the satellite aspect solution, combined with the very low background rate in the detector, permits the detection of sources as weak as $`7`$ photons and the resolution of sources as closely spaced as 2″. This high sensitivity and resolution is illustrated in Fig. 2. For a weakly absorbed ONC star with a typical PMS thermal spectrum with energy $`kT=1`$ keV (see Preibisch, 1997), a limit of 7 photons corresponds to a flux of $`9\times 10^{16}`$ erg s-1 cm-2 ($`0.28`$ keV) or, assuming the ONC distance is 450 pc, a luminosity of $`2\times 10^{28}`$ erg s-1. This is $`25`$ times more sensitive than the ROSAT HRI study of the ONC (Gagné et al., 1995) and $`10`$ times more sensitive than the ROSAT PSPC observation which suffered considerable source confusion (Geier et al., 1995). For comparison, the emission of the contemporary Sun in the same X-ray band ranges from $`2\times 10^{25}`$ erg s-1 in its quiet phase to $`1\times 10^{28}`$ erg s-1 during its most violent flares (Peres et al., 2000). Note, however, that the conversion between count rate and luminosity will be much higher for stars which are heavily absorbed. Accurate luminosities for the entire sample require spectral analysis for each star, which is beyond the scope of the present paper.

To locate sources in a methodical fashion, we apply the source detection program wavdetect based on a Mexican hat wavelet decomposition and reconstruction of the image (Freeman et al., 2000). We set the significance criterion at $`1\times 10^5`$, which permits a few false positive detections but gives good sensitivity. Wavelet scales between 1 and 16 pixels are used. The procedure is nearly complete to a limit of 7 photons near the field center increasing to 10 photons towards the edges of the field; many sources below these completeness limits are also located. A careful examination of the field by eye indicates the wavelet detection algorithm is quite reliable. The principal adjustments to the source lists were the removal of some dozens of spurious wavelet sources associated with the trails of the Trapezium sources and the addition of 19 sources just resolved ($`13`$″separation) from brighter sources. Source detection was performed separately in the soft ($`0.22`$ keV) and hard ($`28`$ keV) bands to improve sensitivity to unembedded and embedded sources respectively, after which the source lists were merged.

This procedure located 831 sources above the completeness limit and an additional 142 likely sources below this limit (Table 1). The sources should be designated CXOONC hhmmss.s-ddmmss, where the acronym stands for ‘Chandra X-ray Observatory Orion Nebula Cluster’. Source count rates range from 0.1 to $`>400`$ cts ks-1 (the upper limit here is due to loss of events due to photon pileup in $`\theta ^1`$C Ori), giving a very high dynamic range. This list of 973 ACIS sources was cross-correlated with a large catalog of $`2500`$ objects in the ONC and its vicinity including complete listings of 1578 stars with visible magnitude $`V<20`$ (Hillenbrand, 1997) and 854 stars from K-band surveys with limiting magnitude $`K<18`$ in the central region (Hillenbrand et al., 1998; Hillenbrand & Carpenter, 2000). We find that 860 of the 973 X-ray sources spatially coincide within 1″$``$2″ of $`V`$\- or $`K`$-band cluster members. Systematic and statistical astrometric errors are $`<1`$″ (see §5).

In addition to the optical-infrared stars, most of the known compact radio sources in the ONC region (Felli et al., 1993, and references therein) are detected in the ACIS image. These radio sources are a heterogeneous mixture of ionized stellar winds or nonthermal emission from both unobscured Trapezium O stars and embedded OB stars in the BN/KL region, magnetically active low mass PMS members, and sources without stellar counterparts such as partially ionized globules within the giant Orion HII region (Garay, 1987; Felli et al., 1993). We discuss some of these sources in §5.

## 3 Demographics of X-ray emission in the ONC

Figure 2 illustrates several types of X-ray source identifications. The majority of the 973 X-ray sources are detected in both the soft and hard X-ray bands and are associated with relatively bright stars seen in the $`V`$ and $`K`$ surveys (sources with circle, square, $`V`$ and $`K`$). Faint unobscured stars may be found only in the soft X-ray band (circle, $`V`$ and $`K`$) while embedded stars are often seen only in the hard X-ray and $`K`$ bands (square and $`K`$). Many of the fainter and more heavily absorbed ONC members, and a handful of unabsorbed stars, are not detected in X-rays ($`V`$ and $`K`$ or $`K`$ alone).

One hundred and three ACIS sources found only in the hard band have no optical or infrared counterpart and are distributed approximately evenly across the ACIS field (square). Some of these are likely extragalactic sources, primarily active galactic nuclei, seen through the molecular cloud. Convolving the flux distribution of hard X-ray band source counts at high-Galactic latitudes (Ann Hornschemeier, private communication; see Hornschemeier et al., 2000) with the distribution of hydrogen column densities across the ACIS field (estimated from a velocity-integrated C18O map using standard conversions to total column densities; John Bally, private communication) predicts $`25`$ hard extragalactic sources in the ACIS field. We conclude that most of the unidentified sources are ONC stars, probably faint embedded late-type PMS stars (§5).

Nine hundred and thirty four of the $`V<20`$ stars are sufficiently well-characterized to be placed on the Hertzsprung-Russell (HR) diagram from which, together with theoretical evolutionary tracks (D’Antona & Mazzitelli, 1997), stellar masses and ages can be estimated (Hillenbrand, 1997). Figure 3 shows this HR diagram with symbols coded for detection and non-detection in the ACIS image. The detection rate in the sample is 91% for stars with masses $`M>0.3`$ M and 75% for stars with $`0.1<M<0.3`$ M. The Chandra sensitivity of $`2\times 10^{28}`$ erg s-1 is thus sufficient to detect nearly all optically selected ONC stars with $`V<20`$.

However, the $`V`$-band sample in the HR diagram systematically excludes the most obscured and least massive stars, particularly brown dwarfs, which are more effectively found at near-infrared wavelengths. Figure 4 compares the ACIS sources to the deep $`H`$ and $`K`$ band samples (Hillenbrand et al., 1998; Hillenbrand & Carpenter, 2000). A strong relationship between X-ray detectability and $`K`$ magnitude is present. Over 50% of cluster members with $`M>1`$ M are detected, even if they suffered absorption as high as $`A_V60`$, equivalent to a hydrogen column density $`N_H1\times 10^{23}`$ cm-2. Stars between $`0.1<M<1`$ M are generally detected only if $`A_V10`$, and substellar objects with $`M<0.08`$ M are rarely detected for any $`A_V`$. Astrophysical interpretation of these results is not simple because mass, luminosity, radius, reddening and the effects of circumstellar disks are intertwined across the color-magnitude diagram. To remove the effect of obscuration, we restrict our examination to low-$`A_V`$ stars having infrared colors $`(HK)<0.5`$. Again a strong relationship between ACIS detection fraction and $`K`$ magnitude is found: the detection fraction declines from about 85% for $`K<11`$ to about 50% for $`11K<13`$, about 35% for $`13K<15`$, and $`<8`$% for $`K>15`$.

A few of the 973 ACIS sources coincide with faint sources with inferred $`M0.05`$ M, assuming an age of 1 My (Figures 3 and 4). These candidate PMS brown dwarfs are PC 116, PC 215, H 5096, HC 114 and HC 123 (Prosser et al., 1994; Hillenbrand & Carpenter, 2000). Two of these have spectral type around M6$``$M7. Taken at face value, these brown dwarf candidates roughly doubles the number of PMS candidate brown dwarfs detected in the X-ray band at levels around $`10^{28}`$ erg s-1 with the ROSAT satellite (Neuhäuser & Comerón, 1998; Neuhäuser et al., 1999). However, because these ONC stars reside close to the $`M=0.08`$ M line, infrared excesses from circumstellar disks may be present, photometric variability may occur, and theoretical models are uncertain, they may in fact evolve into hydrogen-burning stars. The strongest result from the ACIS ONC data is that great majority of PMS brown dwarfs found by Hillenbrand & Carpenter (2000) are not X-ray emitters with $`L_x2\times 10^{28}`$ erg s-1, our detection limit, in agreement with the ROSAT studies.

## 4 Evolution of X-ray emission in solar mass stars

Among solar-type main sequence stars, it is well-established that X-ray emission arises from magnetic activity at the stellar surface, where the magnetic fields are generated by a magnetic dynamo in the stellar interior driven by rotation (Rosner et al., 1985). Thus, main sequence stars show a correlation between $`L_x`$ and surface rotational velocity but little dependence between $`L_x`$ and other stellar properties such as bolometric luminosity $`L_{bol}`$ or mass $`M`$. The situation for pre-main sequence stars has been more confusing: despite considerable evidence that their X-ray emission is related to high levels of surface magnetic activity, $`L_x`$ has been found to be strongly correlated with both $`L_{bol}`$ and $`M`$ (e.g., Walter & Kuhi, 1981; Feigelson et al., 1993; Casanova et al., 1995; Neuhäuser et al., 1995; Feigelson & Montmerle, 1999). Furthermore, as PMS stars evolve downward in the Hertzsprung-Russell diagram during their early convective phases, a dependence on $`L_{bol}`$ may mean that X-ray luminosity decreases with age $`t`$ during PMS evolution. $`L_x`$ may also scale with stellar radius (recall $`RL_{bol}^{1/2}`$ along vertical convective tracks) if reconnecting magnetic fields cover the entire surface of the star. A weak relation between $`L_x`$ and stellar rotation is found in some, but not all, PMS samples.

These relationships have never been convincingly explained. If $`L_xL_{bol}`$, $`L_xM`$ or $`L_xt`$ correlations represent physical links, then PMS magnetic activity may not be based on the standard dynamo but rather on more exotic mechanisms. Such mechanisms include surface activity arising from magnetic fields inherited from the star formation processes rather than generated within the PMS stars, or violent reconnection events in magnetic fields connecting the star to the circumstellar disk (Feigelson & Montmerle, 1999).

To untangle a part of this confusion, we consider only a narrow range of stellar masses around 1 M to remove mass as a confounding variable. Figure 5 shows the plot of ACIS count rate $`vs.`$ bolometric luminosity for a complete sample of 17 ONC stars with $`0.8<M<1.2`$ M and low obscuration ($`A_V<3`$). A well-defined sample of Pleiads in the same mass range observed with the ROSAT satellite, selected to have $`0.52(BV)0.74`$, is shown for comparison (Stauffer et al., 1994; Micela et al., 1996). These studies established that rotational velocity is the dominant contributor to the wide dispersion in Pleiad X-ray emission.

Rather than a simple dependence between $`L_x`$ and $`L_{bol}`$, the X-ray emission of ONC solar-mass stars in Fig. 5 might be divided into two epochs. During the first $`2`$ My as the bolometric luminosity fades by a factor of $`20`$, X-ray luminosities appear high and constant around $`L_x2\times 10^{30}`$ erg s-1. This means the X-ray surface flux ($`L_x/4\pi R^2`$) increases as the star descends the Hayashi track, as found in ROSAT studies (Neuhäuser et al., 1995; Kastner et al., 1997). But as the stars age from $`2`$ to $`10`$ My and enter the radiative track, some stars retain their high X-ray levels while in others it drops by a factor of $`50`$ or more. The sample is small, so this pattern is not established with confidence.

This behavior is in accord with current views of a solar-type magnetic dynamo origin for PMS X-ray emission where rotation regulates the level of magnetic activity at the stellar surface, and where rotation in turn is regulated by a magnetic coupling between the star and its circumstellar disk. In this model, the wide dispersion in $`L_x`$ among older PMS stars reflects the wide dispersion in stellar rotation rates which is expected from different durations of star-disk coupling (Bouvier et al., 1997). This same explanation accounts for the wide dispersion of $`L_x`$ linked to stellar rotation found in zero-age main sequence solar-mass Pleiades stars.

This dynamo interpretation of PMS X-ray emission can be directly tested using the sample shown in Fig. 5: older ONC stars with high $`L_x`$ (JW198, JW601 and JW738) should show faster rotational velocities in high-resolution optical spectra than those with low $`L_x`$ (JW278, JW662 and JW868). If this prediction of the standard magnetic dynamo model is not validated in this and similar subsamples of the ONC, then non-solar models for magnetic flaring in T Tauri stars must be considered.

## 5 X-rays from the embedded BN/KL region

A cluster of embedded stars, including the $`L_{bol}10^4`$ L Becklin-Neugebauer (BN) object lie near the density peak of the massive OMC-1 molecular cloud about 1′ northwest of the Trapezium (Genzel & Stutzki, 1989). BN lies in a group of several $`L_{bol}10^3`$ L massive stars (Gezari et al., 1998, and references therein), three of which (BN, IRc2, and infrared Source n) produce compact thermal radio sources (ultracompact HII regions and/or stellar winds) and/or maser outflows. This group of massive young stars illuminates the infrared Kleinman-Low (KL) nebula and powers strong shocked molecular outflows in the region (Stolovy et al., 1998, and references therein). BN/KL lies in the densest and most chemically rich region in the entire Orion giant molecular cloud complex, and is the nearest and best studied example of $`20`$ known molecular hot cores which are thelikely sites of current massive star formation (Kurtz et al., 2000). The region is extremely complex and not fully understood after decades of intensive multiwavelength investigation.

To obtain a reliable detailed understanding of the X-ray emission in this crowded region, the ACIS data were specially processed for this analysis. First, we apply the software corrector for charge transfer inefficiency (§2) described by Townsley et al. (2000). This restores much of the degradation of event energies and grade distributions in a spatially-dependent fashion across the CCD chip. Second, source detection proceeded as described in §2. Third, source positions of several bright X-ray sources were matched to positions of infrared sources (aligned to the Hipparcos frame through the 2MASS survey; Hillenbrand & Carpenter, 2000) and VLA radio sources (Felli et al., 1993; Menten & Reid, 1995). A 0.5″ boresight correction was applied to the ACIS field resulting in a systematic astrometric accuracy better than $`\pm `$0.1″ (estimated 90% confidence level). Individual source positions have estimated statistical uncertainties ranging from $`\pm 0.1`$″ for strong sources to $`\pm 0.5`$″ for weak sources (estimated 90% confidence).

Fourth, spectral analysis on the corrected data was performed with the XSPEC software package (Arnaud, 1996). We present here the equivalent hydrogen column density associated with line-of-sight absorption. Due to uncertainties in the detector spectral response and the complexities of source spectra, we estimate the 90% confidence intervals to be $`\pm 0.5`$ in log$`N_H`$. Recall that log$`N_H=21.2`$ cm-2 corresponds to $`A_V=1`$ for normal elemental abundances and gas-to-dust ratios (Ryter, 1996). Fifth, photon arrival times were tested for variability using a Kolmogorov-Smirnov one-sample test against a model of a constant source. We report source variations significant at the $`P>99.5`$% confidence level.

Figure 6 shows a full-resolution view of the ACIS field around BN/KL in the $`0.28`$ keV band. Unlike Figures 1 and 2, colors here represent the photon energy where red represents hard photons and blue represents soft unabsorbed photons. Twenty-seven sources above the limit of 7 photons (§2) found with wavelet detection and visual examination are indicated; their X-ray properties and counterparts are given in Table 2 and the Table Notes. Two sources may be poorly resolved doubles with separations around 1.5″. Fainter sources are likely present; for example, a group of 5 photons around ($`\mathrm{\Delta }\alpha ,\mathrm{\Delta }\delta `$)=(12″,$``$16″) appears as a source in Table 1 and coincides with an anonymous $`K1819`$ NICMOS star (Stolovy et al., 1998).

In Table 2, running source numbers shown in Figure 6 are provided in column 1. X-ray centroid positions in columns $`23`$ are in epoch J2000 from the ACIS image after boresight correction. Column 4 gives ACIS counts extracted from a 1.5″ circle in the $`0.28`$ keV band with ASCA grades after correction for charge-transfer inefficiency. A background of 2 counts per source has been subtracted. Column densities rounded to the nearest dex(0.5) in column 5 are estimated from XSPEC fits to X-ray pulse height distributions. Column 6 provides offsets in arcseconds between the ACIS source and its counterpart, based on K-band (Hillenbrand & Carpenter, 2000) or radio positions (Felli et al., 1993, adjusted to the astrometry of Menten & Reid 1995).

We divide sources into three classes based on their counterparts: ‘ONC’ ACIS sources are likely members of the Orion Nebula Cluster, generally optically-bright lightly-absorbed stars; ‘IR’ ACIS sources have infrared but not optical counterparts and are likely embedded young stars; and ‘X-ray’ ACIS sources have no photospheric stellar counterparts. Columns $`89`$ of Table 2 indicate the optical counterparts and magnitudes (JW = Jones & Walker 1988, P = Prosser et al. 1994). Columns $`1011`$ give near-infrared counterparts and K-band magnitudes (JW = Jones & Walker 1988, PC = Prosser et al. 1994, HC = Hillenbrand & Carpenter 2000), except for magnitudes of BN and Source n obtained from Lonsdale et al. (1982). The final column points to Table Notes on individual sources where optical counterpart information is from Hillenbrand (1997), infrared information is from Dougados et al. (1993), and radio information is from Felli et al. (1993), Menten & Reid (1995) and Menten (private communication). X-ray variability and spectral properties are also summarized in the Notes.

Orion Nebula Cluster members
The majority of ACIS sources in the vicinity of BN are relatively unobscured ($`A_V3`$) members of the ONC, as expected given the proximity ($`1`$′ or $`0.1`$ pc) from the Trapezium. The X-ray stars are very young with estimated ages below 1 Myr. Most are lower mass stars with $`M0.10.4`$ M, although two X-ray bright 0.8 M stars are present. The poorly studied late-K star JW 470 = CXOONC J053515.35-052215.8 is the most magnetically active star of this class, exhibiting the strongest and most variable X-ray emission. The identification of CXOONC J053515.07-052231.3 with the unobscured optical star PC 086 is uncertain despite the 0.5″ spatial coincidence, as the X-ray source appears heavily absorbed while the optical star is not. An alternative possible counterpart is the 12.4$`\mu `$m source IRc14 lying 2.0″ from the X-ray location (Gezari et al., 1998).

Two ACIS sources in this class, CXOONC J-53514.08-052237.1 = JW 432 and CXOONC J053514.66-052233.9 = JW 448, are radio sources. These and other radio sources discussed here are generally variable and are classified as FOXES, Fluctuating Optical and X-ray Emitting Sources, by Garay (1987). They resemble the several dozen known weak-lined T Tauri stars seen in nearby star forming regions exhibiting nonthermal gyrosynchrotron radio emission and strong X-ray emission (André, 1996; Feigelson & Montmerle, 1999). These radio/X-ray ONC stars have probably lost their interacting disks at a very young age, resulting in rapid rotation and high levels of magnetic flaring.

Embedded infrared stars
Six ACIS sources are associated with embedded 2$`\mu `$m stars without visible counterparts. Four of these are consistent with embedded but otherwise ordinary low-mass ONC members. CXOONC J053514.35-052232.9, however, coincides with the luminous infrared Source n. This $`L_{bol}10^2`$ L member of the BN/KL group lies at the center of the dense molecular hot core and a H2O maser outflow, and appears to be associated with both a nonthermal double radio source and diffuse thermal radio emission (Menten & Reid, 1995; Chandler & Wood, 1997; Gezari et al., 1998). The unabsorbed X-ray luminosity of order $`10^{30}10^{31}`$ erg s-1 is consistent with a magnetically active low-mass protostar or T Tauri star or with an ordinary early-type B star.

Perhaps most remarkable among these sources is CXOONC J053514.06-052222.2 which lies 1.1″ northwest of the BN object itself. This offset appears larger than the expected positional uncertainties, estimated at $`\pm 0.5`$″ at 90% confidence. But only $`18`$ source photons are present and we cannot be completely confident the offset is real. One possibility is that the X-ray source is unrelated to the BN object but arises from an uncatalogued embedded ONC star; the energy distribution of the detected X-ray photons suggests $`A_V50`$. Another possibility is that the emission arises from outflows from the BN object. The photon distribution is consistent with a source structure $`1`$″ ($`500`$ A.U.) in extent. Hydrodynamical calculations of ultra-compact HII regions with mass-loaded winds from massive young stars pushing into a stationary dense medium indicate that, although the emission is dominated by readily-absorbed ultraviolet photons, some high-energy X-rays are produced and escape from the immediate environment (Strickland & Stevens, 1998, see also Arthur et al. 1993). There are also possibilities of non-equilibrium shock processes giving rise to X-ray emission at the interfaces between HII regions and molecular material (Dorland & Montmerle, 1987).

X-ray sources without optical or infrared counterparts
The five ACIS sources in the ‘X-ray’ class with no photospheric counterparts in the deepest available optical, K-band and mid-infrared surveys of the region (Stolovy et al., 1998; Hillenbrand & Carpenter, 2000; Gezari et al., 1998). All five sources are heavily absorbed with log$`N_H23`$ cm-2 or roughly $`A_V50`$, two exhibit short-timescale X-ray variability, and three are associated with compact, likely variable radio sources (Felli et al., 1993). The two ACIS sources not associated with radio sources have unusually hard X-ray spectra similar to those seen in Class I protostars and T Tauri stars during powerful flares (Koyama et al., 1996; Tsuboi et al., 1998). CXC J053514.90-052225.7 = radio D is the most extreme of these sources. Its luminosity would be $`L_x10^{31}`$ erg s-1 if it were unobscured.

Three explanations for these sources might be considered: background Galactic or extragalactic sources unrelated to the Orion star forming region; X-ray luminous low-mass protostars (one of the dozen known X-ray emitting Class I protostars is a nonthermal radio variable Feigelson et al., 1998); or magnetically active low-mass radio-loud (i.e. weak-lined T Tauri) stars. We estimate that only $`01`$ of these six sources may be an extragalactic object seen through the molecular cloud (§3). Protostars are possible, but perhaps their disks or outflows should have been detected in the extensive infrared and molecular studies of the region. We suggest that the third alternative is most plausible. The X-ray and radio properties are consistent with deeply embedded low-mass PMS stars having unusually high $`L_x/L_{bol}`$ ratios such that their photospheric emission falls below current infrared sensitivity limits ($`K19`$; Stolovy et al., 1998).

It is worth remarking on the absence of certain objects in BN/KL region. The failure to detect the widely-discussed multiple infrared sources comprising IRc2 (Dougados et al., 1993) could be due to their high level of absorption ($`A_V58`$; Gezari et al., 1998). No extended emission is evident in this portion of the ACIS image. X-rays are thus not detected from infrared-luminous dust condensations (IRc3, IRc4, IRc7, infrared Shell or Crescent; Stolovy et al., 1998), or the molecular hot core, plateau or H2 outflows (e.g., Wright et al., 1996). A recent supernova explosion in the region (Kundt & Yar, 1997) is clearly excluded by the low background level seen in Figure 6.

Independent of the nature of the host stars, a clear astrophysical implication of the X-ray emission from deeply embedded PMS stars in the BN/KL region and the wider ONC field is that X-ray dissociation regions (Hollenbach & Tielens, 1997) will be present even in the densest molecular cloud cores. These are volumes around X-ray emitting PMS stars where X-ray ionization dominates cosmic ray ionization. X-ray irradiation on molecular cloud interiors is predicted to have a variety of effects on the dynamics, astrochemistry, and dust characteristics of molecular material (Glassgold et al., 2000), although little evidence for such effects has yet been found.

## 6 Conclusions

A 48 ks Chandra ACIS image of the Orion Nebula in the $`0.28`$ keV band with a limiting sensitivity of $`2\times 10^{28}`$ erg s-1, sub-arcsecond spatial resolution and astrometric precision reveals:

1. About one thousand X-ray sources are associated with magnetically active pre-main sequence stars in the Orion Nebula Cluster. Almost 90% of the visible stars with $`V<20`$ are detected across the entire Hertzsprung-Russell diagram. Higher mass stars are detected even when subject to tens of magnitudes of absorption, while the lowest mass stars and substellar objects are typically not detected.
2. The relationships between X-ray emission and other stellar properties are complex and should be elucidated in a more complete analysis of the ONC. Examination of a well-defined subsample of solar mass stars suggests a high constant X-ray luminosity around $`2\times 10^{30}`$ erg s-1, hence increasing surface flux, as stars descend the convective Hayashi track followed by an increased dispersion in X-ray luminosities as they enter the radiative track. The suggested interpretation, that this dispersion in magnetic activity is due to different rotational evolutions of different stars, can be easily tested with measurements of stellar rotational velocities. These results may be relevant to the evolution of magnetic activity in the early Sun.
3. A close examination of X-ray sources in the BN/KL region of OMC 1 shows a diversity of stellar sources: foreground and embedded ONC members; the luminous embedded infrared Source n; several sources not appearing in sensitive optical and K-band surveys; and a faint, heavily absorbed source apparently displaced 1″ from BN itself. The X-ray spectra indicate line-of-sight absorption ranging from $`21\mathrm{log}N_H23`$ cm-2 and one-fourth of the sources are variable within the observation. A considerable fraction of the sources are associated with faint radio sources; this likely arises from nonthermal gyrosynchrotron radiation that accompanies X-ray flaring induced by magnetic reconnection. The ACIS source associated with Source n may be the first unambiguous (i.e. spatially resolved from other young stars) X-ray detection of a massive young stellar object. The X-ray source near BN is enigmatic. It is unclear whether it is a separate star or whether the X-rays are produced in the outflow from this massive young star. The several sources without any photospheric optical or infrared counterpart are likely new embedded, magnetically active stars in OMC 1. The presence of deeply embedded X-ray sources in the BN/KL region support the idea that ionizing X-rays produced by young stars will be present within dense star-forming molecular cloud cores.

We express our appreciation of the many scientists and engineers who brought Chandra to fruition, in particular those at MIT, Penn State and Lockheed-Martin who contributed to the ACIS instrument. Ann Hornschemeier (Penn State), Karl Menten (MPI-Radioastronomie) and John Bally (Colorado) kindly provided results prior to publication, and the referee Ralph Neuhäuser gave many useful comments. This research was funded by NASA contract NAS8-38252 at Penn State. The contributions of S.H.P. were carried out at the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.


text:
# Instabilities of waves in nonlinear disordered media

## Abstract

We develop a self-consistent theory of temporal fluctuations of a speckle pattern resulting from the multiple scattering of a coherent wave in a weakly nonlinear disordered medium. The speckle pattern is shown to become unstable if the nonlinearity exceeds a threshold value. The instability is due to a feedback provided by the multiple scattering and manifests itself in spontaneous fluctuations of the scattered intensity. The development of instability is independent of the sign of nonlinearity.

Interplay between nonlinearity and multiple scattering has received considerable attention during the recent years. Second harmonic generation, phase conjugation, and optical Kerr effect have been studied in the weak scattering limit ($`k\mathrm{}1`$, $`k`$ — the wavenumber in the medium, $`\mathrm{}`$ — the mean free path). On the other hand, the nonlinearity is shown to modify the coherent backscattering peak, and to affect the angular correlations of scattered waves.

In the present Letter, we consider the temporal fluctuations of waves in weakly nonlinear disordered media. We show that, in contrast to linear media, where the fluctuations of scattered waves with time are due to scatterer motion, the wave scattered in a nonlinear medium may exhibit temporal fluctuations even if the scatterers are immobile. These fluctuations are spontaneous and represent a signature of instability of the multiple-scattering speckle pattern. Although unstable regimes are rather common for nonlinear waves in feedback systems, our study is the first to consider the (multiple) scattering as a feedback-providing mechanism in a nonlinear medium. The long-range spatial correlations of intensity fluctuations appear to be crucial for the feedback efficiency.

We consider a scalar plane wave of frequency $`\omega `$ incident upon a nonlinear disordered medium occupying the half-space $`z>0`$. Wave amplitude in the medium $`\psi (𝐫,t)`$ obeys a nonlinear wave equation:

$`\left\{^2+k_0^2\left[\epsilon _0+\delta \epsilon (𝐫)+\epsilon _2\left|\psi (𝐫,t)\right|^2\right]\right\}\psi (𝐫,t)=J(𝐫,t),`$ (1)
where $`k_0`$ is the free-space wavenumber, $`\epsilon _0=n_0^2`$ is the average dielectric function, $`\delta \epsilon (𝐫)`$ is the fluctuating part of the dielectric function, $`\epsilon _2`$ is a nonlinear constant, and $`J(𝐫,t)`$ is a (slowly varying) source term \[$`\left|(/t)J(𝐫,t)\right|\omega \left|J(𝐫,t)\right|`$\]. Eq. (1) describes propagation of waves in centrosymmetric media, where nonlinear terms quadratic in $`\psi `$ are forbidden. Equations of the same type appear for electromagnetic waves in plasma, as well as in the context of the theory of hydrodynamic turbulence. A nonlinear Schrödinger equation similar to Eq. (1) has been recently considered by Spivak and Zyuzin.

We assume the source term $`J(𝐫,t)`$ in Eq. (1) to consist of two parts: $`J_0(𝐫)`$, representing an ideal plane wave incident at $`z=0`$, and $`\delta J(𝐫,t)`$, which is a random noise term. In the present paper, we are only interested in the case of infinitely weak noise: $`\left|\delta J\right|0`$. In a linear medium ($`\epsilon _2=0`$), the solution of Eq. (1) is then unique and time independent: $`\psi (𝐫,t)=\psi (𝐫)`$. Consequently, the autocorrelation function $`g_1(𝐫,\tau )=\psi (𝐫,t)\psi ^{}(𝐫,t+\tau )/I(𝐫,t)1`$, where $`I(𝐫,t)=\left|\psi (𝐫,t)\right|^2`$ is the wave intensity. If $`\epsilon _20`$, however, Eq. (1) may admit multiple solutions, and (even being infinitely weak) the noise $`\delta J`$ can provoke the system to jump from one solution to another. Obviously, such a behavior will result in spontaneous fluctuations of $`\psi (𝐫,t)`$ with time, yielding $`g_1(𝐫,\tau )<1`$. In the general case, the calculation of $`g_1`$ appears to be too involved. We therefore limit ourselves to long times $`\tau T_{jump}`$, where $`T_{jump}`$ is the typical time needed for the system to jump from one solution to another. Calculation of $`T_{jump}`$ is beyond the scope of our analysis, but we expect $`T_{jump}`$ to be inversely proportional to the speed of wave in the medium. For $`\tau T_{jump}`$, $`g_1(𝐫,\tau )`$ should become independent of $`\tau `$, and its asymptotic value can be considered as a measure of the mutual resemblance of solutions of Eq. (1). To simplify the further analysis, we adopt the following four assumptions:

(i). Propagation of waves in a weakly nonlinear disordered medium is diffusive with a mean free path $`\mathrm{}`$ unaffected by nonlinearity. This implies that nonlinear refraction is negligible at distances of order $`\mathrm{}`$, and consequently, that $`\mathrm{\Delta }n^2k\mathrm{}1`$, where $`\mathrm{\Delta }n=n_2I_0`$, $`n_2=\epsilon _2/(2n_0)`$, and $`I_0=I(𝐫,t)`$ in the absence of absorption.

(ii). Statistics of the wave field scattered in a weakly nonlinear disordered medium is close to Gaussian. Consequently, the factorization approximation $`g_2(𝐫,\tau )=I(𝐫,t)I(𝐫,t+\tau )/I(𝐫,t)^21+\left|g_1(𝐫,\tau )\right|^2`$ applies. This assumption is justified under the same conditions as (i), since the Gaussian statistics of the total scattered wave field $`\psi (𝐫,t)`$ is a consequence of the complete randomization of phases of partial waves arriving at $`𝐫`$. The reason for the randomization is that the typical distance $`\mathrm{}`$ between individual scattering events in a multiple-scattering sequence is much larger than the wavelength ($`k\mathrm{}1`$). Obviously, such a mechanism of phase randomization is equally effective in both linear and weakly nonlinear media, as long as (i) holds.

(iii). Intensity of the third harmonic remains always much smaller than the intensity of the fundamental wave. This implies that the medium has a sufficient degree of dispersion for the phase matching condition to be violated: $`\left|\mathrm{\Delta }k\right|\mathrm{}1`$ with $`\mathrm{\Delta }k=k_33k`$ and $`k_3`$ being the wavenumber at frequency $`3\omega `$.

(iv). Functional form of the $`𝐫`$-dependence of $`g_1`$ is the same as in the linear medium: $`g_1(𝐫,\tau )=\mathrm{exp}[\beta (\tau )z/\mathrm{}]`$, where $`\beta 1`$ is some unknown function. To justify this ansatz, we make use of (i) and write $`g_1`$ as a sum of contributions of paths of various lengths $`s`$:

$`g_1(𝐫,\tau )={\displaystyle _0^{\mathrm{}}}dsP(z,s)\mathrm{exp}\left[{\displaystyle \frac{1}{2}}\mathrm{\Delta }\phi ^2(\tau )_s\right]{\displaystyle _0^{s_1(\tau )}}dsP(z,s).`$ (2)
Here $`P(z,s)zs^{3/2}\mathrm{exp}[3z^2/(4\mathrm{}s)]`$ is the relative “weight” of paths of length $`s`$, and $`\mathrm{\Delta }\phi ^2(\tau )_s`$ is the variance of the phase difference between two waves propagating along the same path but separated in time by $`\tau `$, averaged over all possible realizations of disorder and over all paths of the same length $`s`$. From here on, we denote such an averaging by $`\mathrm{}_s`$. As a function of $`s`$, $`P(z,s)`$ is peaked at $`s_0=z^2/(2\mathrm{})`$. The second equality of Eq. (2) is obtained by assuming that $`\mathrm{\Delta }\phi ^2(\tau )_s`$ increases monotonically with $`s`$ and becomes of order unity at $`s=s_1(\tau )s_0`$. The result of integration in Eq. (2) can be written as $`1\beta (\tau )z/\mathrm{}\mathrm{exp}[\beta (\tau )z/\mathrm{}]`$ with $`\beta z/\mathrm{}z(\mathrm{}s_1)^{1/2}=(2s_0/s_1)^{1/2}1`$. In the opposite limit of $`\beta z/\mathrm{}1`$, $`g_1(𝐫,\tau )`$ vanishes and the functional form of its $`z`$-dependence is of no importance. Anyway, it will be seen from the following that the exact functional form of the $`𝐫`$-dependence of $`g_1(𝐫,\tau )`$ is not of crucial importance, since $`g_1`$ is integrated over the whole medium during the calculation of the correlation function of diffusely reflected wave.

Consider now a typical diffusely reflected wave path of length $`s`$. The phase $`\phi (t)`$ accumulated along such a path is

$`\phi (t)=k_0{\displaystyle _0^s}ds^{}\left\{n_0+n_2I[𝐫(s^{}),t]\right\},`$ (3)
where the integration is along the path. The variance of the phase difference $`\mathrm{\Delta }\phi (\tau )=\phi (t+\tau )\phi (t)`$, averaged over all possible paths of the same length $`s`$, is found directly from Eq. (3) as a sum of two contributions:

$`\mathrm{\Delta }\phi ^2(\tau )_s^{(1)}=(\pi /n_0)k_0\mathrm{}n_2^2\mathrm{\Delta }I(𝐫,\tau )^2_ss/\mathrm{},`$ (4)
$`\mathrm{\Delta }\phi ^2(\tau )_s^{(2)}=(k_0\mathrm{})^2n_2^2\mathrm{\Delta }I(𝐫_1,\tau )\mathrm{\Delta }I(𝐫_2,\tau )_s(s/\mathrm{})^2,`$ (5)
where $`\mathrm{\Delta }I(𝐫,\tau )=I(𝐫,t+\tau )I(𝐫,t)`$. $`\mathrm{\Delta }\phi ^2(\tau )_s^{(1)}`$ and $`\mathrm{\Delta }\phi ^2(\tau )_s^{(2)}`$ originate from the short- and the long-range spatial correlations of $`\mathrm{\Delta }I`$, respectively.

If $`\tau T_{jump}`$, the calculation of the r.h.s. averages in Eqs. (4) and (5) is straightforward. To evaluate $`\mathrm{\Delta }I(𝐫,\tau )^2_s`$, we integrate $`\mathrm{\Delta }I(𝐫,\tau )^2=2I(𝐫,t)^2[1g_1^2(𝐫,\tau )]`$ over $`z>0`$ with a weight function $`\rho _s(𝐫)`$ describing the density of paths of length $`s`$ in the medium. To find the average of Eq. (5), we calculate the long-range correlation function $`\mathrm{\Delta }I(𝐫_1,\tau )\mathrm{\Delta }I(𝐫_2,\tau )`$ using the Langevin approach with the time correlation function of Langevin random sources being proportional to $`g_1^2(𝐫,\tau )`$. The result is then averaged over all pairs of points $`𝐫_1`$, $`𝐫_2`$ belonging to the same path with a weight function $`\rho _s(𝐫_1,𝐫_2)`$ which is calculated similarly to $`\rho _s(𝐫)`$. As $`g_1=\mathrm{exp}(\beta z/\mathrm{})`$ enters into the derivations of $`\mathrm{\Delta }I(𝐫,\tau )^2_s`$ and $`\mathrm{\Delta }I(𝐫_1,\tau )\mathrm{\Delta }I(𝐫_2,\tau )_s`$, $`\mathrm{\Delta }\phi ^2(\tau )_s^{(1),(2)}`$ appear to depend explicitly on $`\beta `$. The details of the calculations will be presented elsewhere, the final results are:

$`\mathrm{\Delta }\phi ^2(\tau )_s^{(1)}=2\pi k_0\mathrm{}\mathrm{\Delta }n^2H\left(\beta \sqrt{{\displaystyle \frac{s}{3\mathrm{}}}}\right){\displaystyle \frac{s}{\mathrm{}}},`$ (6)
$`\mathrm{\Delta }\phi ^2(\tau )_s^{(2)}6\mathrm{\Delta }n^2\times \{\begin{array}{cc}\beta \times (s/\mathrm{})^2,\hfill & \beta \sqrt{s/\mathrm{}}1\text{,}\hfill \\ (s/\mathrm{})^{3/2},\hfill & \beta \sqrt{s/\mathrm{}}>1\text{,}\hfill \end{array}`$ (7)
where $`H(x)=\sqrt{\pi }x\mathrm{exp}(x^2)[1\mathrm{Erf}(x)]`$, and $`n_0=1`$ is assumed for simplicity. Applying Eq. (2), we obtain the autocorrelation function of diffusely reflected wave:

$`g_1(\mathrm{},\tau )=\underset{¯}{\mathrm{exp}\left[\beta \left(\tau \right)\right]}={\displaystyle _0^{\mathrm{}}}dsP(\mathrm{},s)\mathrm{exp}\left[{\displaystyle \frac{1}{2}}\mathrm{\Delta }\phi ^2(\tau )_s\right]=\underset{¯}{F\left[\beta \left(\tau \right)\right]},`$ (8)
where $`\mathrm{\Delta }\phi ^2(\tau )_s`$ is a sum of the terms given by Eqs. (6) and (7). Since the integral of Eq. (8) does not contain $`\tau `$ in explicit form, we conclude that $`\beta `$ and $`g_1`$ are independent of $`\tau `$, as expected for $`\tau T_{jump}`$. From here on, $`g_1`$ without arguments denotes the value of the autocorrelation function of diffusely reflected wave at $`\tau T_{jump}`$.

In Eq. (8), we have introduced a new function $`F(\beta )`$ which depends on $`\mathrm{\Delta }n`$ and $`k_0\mathrm{}`$ as well. The underlined part of Eq. (8) can be considered as a nonlinear equation for $`\beta `$, or, equivalently, for $`g_1`$ since $`\beta =\mathrm{ln}g_1`$ by definition. Making a simultaneous plot of $`F(\beta )`$ and $`\mathrm{exp}(\beta )`$ (see the lower solid curve in the inset of Fig. 1), we find that Eq. (8) has two solutions: $`\beta =0`$ and $`\beta >1`$. The first solution is “unstable” with respect to the weak noise $`\delta J`$ in Eq. (1), since $`\delta J0`$ introduces an additional dephasing term into Eq. (8), yielding $`F(0)<1`$. It is then the second solution $`\beta >1`$ which is the physically realizable one. As $`\beta >0`$ implies $`g_1<1`$, our results show that $`\psi (𝐫,t)`$ exhibits spontaneous fluctuations, and that the speckle pattern resulting from the multiple scattering of coherent wave in a weakly nonlinear disordered medium is unstable.

The value of $`g_1`$, characterizing the magnitude of spontaneous fluctuations ($`g_1=1`$ in the absence of fluctuations), can be found by a numerical solution of Eq. (8) (see the dashed line in Fig. 1), as well as analytically. Assuming $`1g_11`$, and noting that the ultimate contribution to $`\mathrm{\Delta }\phi ^2(\tau )_s`$ is given either by Eq. (6) \[for $`s/\mathrm{}(k_0\mathrm{})^2`$\] or by Eq. (7) \[for $`s/\mathrm{}(k_0\mathrm{})^2`$\], we obtain

$`g_11\{\begin{array}{cc}2\left|\mathrm{\Delta }n\right|^{2/3},\hfill & \left|\mathrm{\Delta }n\right|<(k_0\mathrm{})^{3/2}\text{,}\hfill \\ 3\left|\mathrm{\Delta }n\right|\sqrt{k_0\mathrm{}},\hfill & \left|\mathrm{\Delta }n\right|>(k_0\mathrm{})^{3/2}\text{.}\hfill \end{array}`$ (9)
Up to now, we have considered a nonabsorbing medium. Absorption, however, is inevitable in a real medium and should therefore be necessarily taken into account. This can be done by allowing $`\epsilon _0`$ in Eq. (1) to have a small imaginary part. Eq. (6) then becomes

$`\mathrm{\Delta }\phi ^2(\tau )_s^{(1)}=2\pi k_0\mathrm{}\mathrm{\Delta }n^2\left\{H\left[\left(\beta +{\displaystyle \frac{\mathrm{}}{L_a}}\right)\sqrt{{\displaystyle \frac{s}{3\mathrm{}}}}\right]H\left[{\displaystyle \frac{\mathrm{}}{L_a}}\sqrt{{\displaystyle \frac{s}{3\mathrm{}}}}\right]\right\}{\displaystyle \frac{s}{\mathrm{}}},`$ (10)
where $`L_a=[\mathrm{}/(3k_0\mathrm{Im}\epsilon _0)]^{1/2}\mathrm{}`$ is the absorption length, and $`n_0^2=\mathrm{Re}\epsilon _0=1`$. Eq. (7) is still approximately valid for $`(\mathrm{}/L_a)(s/\mathrm{})^{1/2}1`$, while for $`(\mathrm{}/L_a)(s/\mathrm{})^{1/2}>1`$ we get

$`\mathrm{\Delta }\phi ^2(\tau )_s^{(2)}6\mathrm{\Delta }n^2\times \{\begin{array}{cc}\beta \times (L_a/\mathrm{})^3\sqrt{s/\mathrm{}},\hfill & \beta \sqrt{s/\mathrm{}}1\text{,}\hfill \\ (L_a/\mathrm{})^3,\hfill & \beta \sqrt{s/\mathrm{}}>1\text{.}\hfill \end{array}`$ (11)
Graphical solution of Eq. (8) with $`P(\mathrm{},s)s^{3/2}\mathrm{exp}[3\mathrm{}/(4s)s\mathrm{}/(3L_a^2)]`$ is presented in the inset of Fig. 1 for several absorption lengths $`L_a`$ . As follows from the figure, the absorption has to be weak for the second solution of Eq. (8), $`\beta >0`$, to exist. Sufficiently strong absorption prevents spontaneous fluctuations of $`\psi (𝐫,t)`$, stabilizing the speckle pattern. The condition for the onset of the speckle pattern instability is $`(/\beta )F(\beta )=1`$ evaluated at $`\beta =0`$. This yields:

$`p=\mathrm{\Delta }n^2\left({\displaystyle \frac{L_a}{\mathrm{}}}\right)^2\left[k_0\mathrm{}+{\displaystyle \frac{L_a}{\mathrm{}}}\right]1,`$ (12)
where we introduce a control parameter $`p`$. If $`p<1`$, $`g_1=1`$ and the multiple-scattering speckle pattern is stable, while for $`p>1`$, an instability shows up leading to $`g_1<1`$ (see Fig. 1). While $`g_1`$ characterizes the mutual resemblance of solutions of Eq. (1), the number of solutions grows exponentially with $`p^{3/4}`$, as can be estimated using the method of Ref. . The striking feature of Eq. (12) is that $`p`$ can reach unity even for very small $`\left|\mathrm{\Delta }n\right|`$ provided that the extensive parameter $`L_a/\mathrm{}`$ is large enough. The stability condition $`p<1`$ is therefore never satisfied in a nonabsorbing medium. Eq. (12) is consistent with the result of Spivak and Zyuzin, who have shown that the perturbation analysis of the sensitivity of speckle pattern in a nonlinear disordered medium to changes of scattering potential fails for $`\mathrm{\Delta }n^2(L/\mathrm{})^3>1`$, where $`L`$ is the typical size of the medium.

The physical origin of the instability of speckle pattern for $`p>1`$ can be revealed by realizing that the system “coherent wave + nonlinear disordered medium” has a positive distributed feedback. An infinitely weak perturbation of the wave intensity $`I(𝐫,t)`$ produces a change of the local refractive index, altering the phases of waves propagating in the medium, and affecting their mutual interference. Since it is this interference which determines $`I(𝐫,t)`$, the feedback loop is built up and closed. For $`p1`$, the feedback is sufficiently strong to compensate for the (diffusive on average) spreading of the initial intensity perturbation, and the speckle pattern $`I(𝐫,t)`$ becomes unstable.

The above analysis can be generalized to the case of $`\delta \epsilon `$ which varies slowly with time (moving scatterers). An additional contribution to $`\mathrm{\Delta }\phi ^2(\tau )_s`$ arises in this case. If, for instance, the scatterers undergo Brownian motion with a diffusion coefficient $`D`$, this contribution is

$`\mathrm{\Delta }\phi ^2(\tau )_s^{(0)}={\displaystyle \frac{\tau }{\tau _0}}\left\{1+2\mathrm{\Delta }n\left[1H\left({\displaystyle \frac{\mathrm{}}{L_a}}\sqrt{{\displaystyle \frac{s}{12\mathrm{}}}}\right)\right]\right\}{\displaystyle \frac{s}{\mathrm{}}},`$ (13)
where both $`\tau _0=(4k_0^2D)^1`$ and $`\tau `$ are assumed to be much larger than $`T_{jump}`$. The time autocorrelation functions of diffusely reflected waves obtained by solving Eq. (8) with $`\mathrm{\Delta }\phi ^2(\tau )_s`$ given by a sum of contributions of Eqs. (7), (10), (11), and (13) are shown in Fig. 2. The two lower curves correspond to the unstable regime ($`p>1`$), and therefore $`g_1(\tau )<1`$ even if $`\tau /\tau _0=0`$. For the upper curve, the speckle pattern is stable ($`p<1`$) and $`g_1(\tau )=1`$ at $`\tau /\tau _0=0`$.

In conclusion, we have shown that the multiple-scattering speckle pattern may become unstable if the disordered medium exhibits a weak nonlinearity of Kerr type. The instability is due to the distributed feedback provided by the scattering, and its development is independent of the sign of nonlinearity. The instability threshold decreases as absorption is reduced, and the speckle pattern is never stable in a nonabsorbing unbounded medium. Finally, we note that our results may be of particular importance in the context of “random laser”, where absorption may be efficiently compensated by stimulated emission.