text:
# References

Revised version
June 2003

Reviewing the problem of the $`U(1)`$ axial symmetry
and the chiral transition in QCD

M. Marchi, E. Meggiolaro
Dipartimento di Fisica,

Università di Pisa,

Via Buonarroti 2,

I–56127 Pisa, Italy.

Abstract

We discuss the role of the $`U(1)`$ axial symmetry for the phase structure of QCD at finite temperature. We expect that, above a certain critical temperature, also the $`U(1)`$ axial symmetry will be (effectively) restored. We will try to see if this transition has (or has not) anything to do with the usual chiral transition: various possible scenarios are discussed. In particular, supported by recent lattice results, we analyse a scenario in which a $`U(1)`$–breaking condensate survives across the chiral transition. This scenario can be consistently reproduced using an effective Lagrangian model. The effects of the $`U(1)`$ chiral condensate on the slope of the topological susceptibility in the full theory with quarks are studied: we find that this quantity (in the chiral limit of zero quark masses) acts as an order parameter for the $`U(1)`$ axial symmetry above the chiral transition. Further information on the new $`U(1)`$ chiral order parameter is derived from the study (at zero temperature) of the radiative decays of the “light” pseudoscalar mesons in two photons: a comparison of our results with the experimental data is performed.

(PACS codes: 12.38.Aw, 12.39.Fe, 11.15.Pg, 11.30.Rd)

1. Introduction

It is generally believed that a phase transition which occurs in QCD at a finite temperature is the restoration of the spontaneously broken $`SU(L)SU(L)`$ chiral symmetry in association with $`L`$ massless quarks. At zero temperature the chiral symmetry is broken spontaneously by the condensation of $`q\overline{q}`$ pairs and the $`L^21`$ $`J^P=0^{}`$ mesons are just the Goldstone bosons associated with this breaking . At high temperatures the thermal energy breaks up the $`q\overline{q}`$ condensate, leading to the restoration of chiral symmetry. We expect that this property not only holds for massless quarks but also continues for a small mass region. The order parameter for the chiral symmetry breaking is apparently $`\overline{q}q_{i=1}^L\overline{q}_iq_i`$: the chiral symmetry breaking corresponds to the non–vanishing of $`\overline{q}q`$ in the chiral limit $`sup(m_i)0`$. From lattice determinations of the chiral order parameter $`\overline{q}q`$ one knows that the $`SU(L)SU(L)`$ chiral phase transition temperature $`T_{ch}`$, defined as the temperature at which the chiral condensate $`\overline{q}q`$ goes to zero (in the chiral limit $`sup(m_i)0`$), is nearly equal to the deconfining temperature $`T_c`$ (see, e.g., Ref. ). But this is not the whole story: QCD possesses not only an approximate $`SU(L)SU(L)`$ chiral symmetry, for $`L`$ light quark flavours, but also a $`U(1)`$ axial symmetry (at least at the classical level) . The role of the $`U(1)`$ symmetry for the finite temperature phase structure has been so far not well studied and it is still an open question of hadronic physics whether the fate of the $`U(1)`$ chiral symmetry of QCD has or has not something to do with the fate of the $`SU(L)SU(L)`$ chiral symmetry.

In the “Witten–Veneziano mechanism”  for the resolution of the $`U(1)`$ problem, a fundamental role is played by the so–called “topological susceptibility” in a QCD without quarks, i.e., in a pure Yang–Mills (YM) theory, in the large–$`N_c`$ limit ($`N_c`$ being the number of colours):

$$A=\underset{k0}{lim}\underset{N_c\mathrm{}}{lim}\left\{id^4xe^{ikx}TQ(x)Q(0)\right\},$$
(1.1)
where $`Q(x)=\frac{g^2}{64\pi ^2}\epsilon ^{\mu \nu \rho \sigma }F_{\mu \nu }^aF_{\rho \sigma }^a`$ is the so–called “topological charge density”. This quantity enters into the expression for the mass of the $`\eta ^{}`$. Therefore, in order to study the role of the $`U(1)`$ axial symmetry for the full theory at non–zero temperatures, one should consider the YM topological susceptibility $`A(T)`$ at a given temperature $`T`$, formally defined as in Eq. (1.1), where now $`\mathrm{}`$ stands for the expectation value in the full theory at the temperature $`T`$ .

The problem of studying the behaviour of $`A(T)`$ as a function of the temperature $`T`$ was first addressed, in lattice QCD, in Refs. . Recent lattice results  (obtained for the $`SU(3)`$ pure–gauge theory) show that the YM topological susceptibility $`A(T)`$ is approximately constant up to the critical temperature $`T_cT_{ch}`$, it has a sharp decrease above the transition, but it remains different from zero up to $`1.2T_c`$. We recall that, in the Witten–Veneziano mechanism , a (no matter how small!) value different from zero for $`A`$ is related to the breaking of the $`U(1)`$ axial symmetry, since it implies the existence of a would–be Goldstone particle with the same quantum numbers of the $`\eta ^{}`$ (see also Ref. ).

Another way to address the same question is to look at the behaviour at non–zero temperatures of the susceptibilities related to the propagators for the following meson channels  (we consider for simplicity the case of $`L=2`$ light flavours): the isoscalar $`I=0`$ scalar channel $`\sigma `$ (also known as $`f_0`$ in the modern language of hadron spectroscopy), interpolated by the operator $`O_\sigma =\overline{q}q`$; the isovector $`I=1`$ scalar channel $`\delta `$ (also known as $`a_0`$), interpolated by the operator $`\stackrel{}{O}_\delta =\overline{q}\frac{\stackrel{}{\tau }}{2}q`$; the isovector $`I=1`$ pseudoscalar channel $`\pi `$, interpolated by the operator $`\stackrel{}{O}_\pi =i\overline{q}\gamma _5\frac{\stackrel{}{\tau }}{2}q`$; the isoscalar $`I=0`$ pseudoscalar channel $`\eta ^{}`$, interpolated by the operator $`O_\eta ^{}=i\overline{q}\gamma _5q`$. Under $`SU(2)_A`$ transformations, $`\sigma `$ is mixed with $`\pi `$: thus the restoration of this symmetry at $`T_{ch}`$ requires identical correlators for these two channels. Another $`SU(2)`$ chiral multiplet is $`(\delta ,\eta ^{})`$. On the contrary, under the $`U(1)_A`$ transformations, $`\pi `$ is mixed with $`\delta `$: so, an “effective restoration” of the $`U(1)`$ axial symmetry should imply that these two channels become degenerate, with identical correlators. Another $`U(1)`$ chiral multiplet is $`(\sigma ,\eta ^{})`$. (Clearly, if both chiral symmetries are restored, then all $`\pi `$, $`\eta ^{}`$, $`\sigma `$ and $`\delta `$ correlators should become the same.) In practice, one can construct, for each meson channel $`f`$, the corresponding chiral susceptibility

$$\chi _f=d^4xO_f(x)O_f^{}(0),$$
(1.2)
and then define two order parameters: $`\chi _{SU(2)SU(2)}\chi _\sigma \chi _\pi `$, and $`\chi _{U(1)}\chi _\delta \chi _\pi `$. If an order parameter is non–zero in the chiral limit, then the corresponding symmetry is broken. Present lattice data for these quantities seem to indicate that the $`U(1)`$ order parameter survives across $`T_{ch}`$, up to $`1.2T_{ch}`$, where the $`\delta `$$`\pi `$ splitting is small but still different from zero . In terms of the left–handed and right–handed quark fields,***Throughout this paper we use the following notations: $`q_{L,R}\frac{1}{2}(1\pm \gamma _5)q`$, with $`\gamma _5i\gamma ^0\gamma ^1\gamma ^2\gamma ^3`$. one has the following expression for the difference between the correlators for the $`\delta ^+`$ and $`\pi ^+`$ channels:

$`𝒟_{U(1)}(x)O_{\delta ^+}(x)O_{\delta ^+}^{}(0)O_{\pi ^+}(x)O_{\pi ^+}^{}(0)`$ (1.3)
$`=2\left[\overline{u}_Rd_L(x)\overline{d}_Ru_L(0)+\overline{u}_Ld_R(x)\overline{d}_Lu_R(0)\right].`$
(The integral of this quantity, $`d^4x𝒟_{U(1)}(x)`$, is just equal to the $`U(1)`$ chiral susceptibility $`\chi _{U(1)}=\chi _\delta \chi _\pi `$.) What happens below and above $`T_{ch}`$? Below $`T_{ch}`$, in the chiral limit $`sup(m_i)0`$, the left–handed and right–handed components of a given light quark flavour (up or down, in our case with $`L=2`$) can be connected through the $`q\overline{q}`$ chiral condensate, giving rise to a non–zero contribution to the quantity $`𝒟_{U(1)}(x)`$ in Eq. (1.3) (i.e., to the quantity $`\chi _{U(1)}`$). But above $`T_{ch}`$ the $`q\overline{q}`$ chiral condensate is zero: so, how can the quantity $`𝒟_{U(1)}(x)`$ (i.e., the quantity $`\chi _{U(1)}`$) be different from zero also above $`T_{ch}`$, as indicated by present lattice data? The only possibility in order to solve this puzzle seems to be that of requiring the existence of a genuine four–fermion local condensate, which is an order parameter for the $`U(1)`$ axial symmetry and which remains different from zero also above $`T_{ch}`$. This new condensate will be discussed in Section 3. The rest of the paper will be essentially devoted to the analysis of some interesting phenomenological consequences deriving from this hypothesis.

The paper is organized as follows. In Section 2 we discuss the role of the $`U(1)`$ axial symmetry for the phase structure of QCD at finite temperature. One expects that, above a certain critical temperature, also the $`U(1)`$ axial symmetry will be (effectively) restored. We will try to see if this transition has (or has not) anything to do with the usual chiral transition: various possible scenarios are discussed. In particular, supported by the above–mentioned lattice results, in Sections 3 and 4 we analyse a scenario in which a new $`U(1)`$–breaking condensate survives across the chiral transition and it is still present above $`T_{ch}`$. This scenario can be consistently reproduced using an effective Lagrangian model, which also includes the new $`U(1)`$ chiral condensate. This theoretical model was originally proposed in Refs.  and is summarized in Sections 3 and 4 for the convenience of the reader. (See also Refs.  for a recent review on these problems.) In Section 5 (which, together with Section 6, contains the main original results of this paper) we analyse the consequences of our theoretical model on the slope of the topological susceptibility $`\chi ^{}`$, in the full theory with quarks, showing how this quantity is modified by the presence of a new $`U(1)`$ chiral order parameter: we will find that $`\chi ^{}`$ (in the chiral limit $`sup(m_i)0`$) acts as an order parameter for the $`U(1)`$ axial symmetry above $`T_{ch}`$. Further information on the new $`U(1)`$ chiral order parameter is derived in Section 6 from the study (at zero temperature) of the radiative decays of the “light” pseudoscalar mesons in two photons: a comparison of our results with the experimental data is performed. Finally, the conclusions and an outlook are given in Section 7.

2. The phase transitions of QCD

One expects that, above a certain critical temperature, also the $`U(1)`$ axial symmetry will be (effectively) restored.***Of course, the $`U(1)`$ axial symmetry is explicitly broken by the anomaly, which never disappears: the precise meaning of this “effective restoration” will be explained below. We will try to see if this transition has (or has not) anything to do with the usual chiral transition. Let us define the following temperatures:

* $`T_{ch}`$: the temperature at which the chiral condensate $`\overline{q}q`$ goes to zero. The chiral symmetry $`SU(L)SU(L)`$ is spontaneously broken below $`T_{ch}`$ and it is restored above $`T_{ch}`$.
* $`T_\chi `$: the temperature at which the pure–gauge topological susceptibility $`A`$ drops to zero. Present lattice results indicate that $`T_\chi T_{ch}`$ .
* $`T_{U(1)}`$: the temperature at which the $`U(1)`$ axial symmetry is (effectively) restored, meaning that, for $`T>T_{U(1)}`$, there are no $`U(1)`$–breaking condensates. If $`\overline{q}q0`$ also the $`U(1)`$ axial symmetry is broken, i.e., the chiral condensate is an order parameter also for the $`U(1)`$ axial symmetry. Therefore we must have: $`T_{U(1)}T_{ch}`$. Moreover, the Witten–Veneziano mechanism implies that $`T_{U(1)}T_\chi `$, since, after all, the pure–YM topological susceptibility $`A`$ is a $`U(1)`$–breaking condensate.

The following scenario, that we will call “SCENARIO 1”, in which $`T_\chi <T_{ch}`$, is, therefore, immediately ruled out. In this case, in the range of temperatures between $`T_\chi `$ and $`T_{ch}`$ the anomaly effects are absent, but the $`U(1)`$ axial symmetry is still broken by the chiral condensate. In other words, in this range of temperatures the $`U(1)`$ axial symmetry is spontaneously broken (à la Goldstone) and the $`\eta ^{}`$ is the corresponding Goldstone boson, i.e., it is massless in the chiral limit $`sup(m_i)0`$, or, at least, as light as the pion $`\pi `$, when including the quark masses. This scenario was first discussed (and indeed really supported!) in Ref. .

Therefore, we are left essentially with the two following scenarios.
SCENARIO 2: $`T_{ch}T_{U(1)}`$, with $`T_{ch}T_\chi T_{U(1)}`$. If $`T_{ch}=T_\chi =T_{U(1)}`$, then, in the case of $`L=2`$ light flavours, the restored symmetry across the transition is $`U(1)_ASU(2)_LSU(2)_RO(2)O(4)`$, which may yield a first–order phase transition (see, for example, Ref. ).
SCENARIO 3: $`T_{ch}T_{U(1)}`$, that is, the complete $`U(L)_LU(L)_R`$ chiral symmetry is restored only well inside the quark–gluon plasma domain. In the case of $`L=2`$ light flavours, we then have at $`T=T_{ch}`$ the restoration of $`SU(2)_LSU(2)_RO(4)`$. Therefore, we can have a second–order phase transition with the $`O(4)`$ critical exponents. $`L=2`$ QCD at $`TT_{ch}`$ and the $`O(4)`$ spin system should belong to the same universality class. An effective Lagrangian describing the softest modes is essentially the Gell-Mann–Levy linear sigma model, the same as for the $`O(4)`$ spin systems (see Ref. ). If this scenario is true, one should find the $`O(4)`$ critical indices for the $`q\overline{q}`$ chiral condensate and the specific heat: $`\overline{q}q|(TT_{ch})/T_{ch}|^{0.38\pm 0.01}`$, and $`C(T)|(TT_{ch})/T_{ch}|^{0.19\pm 0.06}`$. Present lattice data partially support these results.

3. The $`U(1)`$ chiral order parameter

We make the assumption (discussed in the previous sections) that the breaking/restoration of the $`U(1)`$ chiral symmetry is completely independent of the $`SU(L)SU(L)`$ symmetry. The usual chiral order parameter $`\overline{q}q`$ is an order parameter both for $`SU(L)SU(L)`$ and for $`U(1)_A`$: when it is different from zero, $`SU(L)SU(L)`$ is broken down to $`SU(L)_V`$ and also $`U(1)_A`$ is broken. Thus we need another quantity which could be an order parameter only for the $`U(1)`$ chiral symmetry . The most simple quantity of this kind was found by ’tHooft in Ref. . For a theory with $`L`$ light quark flavours, it is a $`2L`$–fermion interaction that has the chiral transformation properties of:

$$_{eff}\underset{st}{det}(\overline{q}_{sR}q_{tL})+\underset{st}{det}(\overline{q}_{sL}q_{tR}),$$
(3.1)
where $`s,t=1,\mathrm{},L`$ are flavour indices, but the colour indices are arranged in a more general way (see Refs. ). It is easy to verify that $`_{eff}`$ is invariant under $`SU(L)SU(L)U(1)_V`$, while it is not invariant under $`U(1)_A`$. To obtain an order parameter for the $`U(1)`$ chiral symmetry, one can simply take the vacuum expectation value of $`_{eff}`$: $`C_{U(1)}=_{eff}`$. The arbitrarity in the arrangement of the colour indices can be removed if we require that the new $`U(1)`$ chiral condensate is “independent” of the usual chiral condensate $`\overline{q}q`$, as explained in Refs. . In other words, the condensate $`C_{U(1)}`$ is chosen to be a genuine $`2L`$–fermion condensate, with a zero “disconnected part”, the latter being the contribution proportional to $`\overline{q}q^L`$, corresponding to retaining the vacuum intermediate state in all the channels and neglecting the contributions of all the other states. As a remark, we observe that the condensate $`C_{U(1)}`$ so defined turns out to be of order $`𝒪(g^{2L2}N_c^L)=𝒪(N_c)`$ in the large–$`N_c`$ expansion, exactly as the chiral condensate $`\overline{q}q`$.

The existence of a new $`U(1)`$ chiral order parameter has of course interesting physical consequences, which can be revealed by analysing some relevant QCD Ward Identities (WI’s) (see Ref.  and also Ref. ). In the case of the $`SU(L)SU(L)`$ chiral symmetry, one immediately derives the following WI:

$$d^4xT^\mu A_\mu ^a(x)i\overline{q}\gamma _5T^bq(0)=i\delta _{ab}\frac{1}{L}\overline{q}q,$$
(3.2)
where $`A_\mu ^a=\overline{q}\gamma _\mu \gamma _5T^aq`$ are the $`SU(L)`$ axial currents. If $`\overline{q}q0`$ (in the chiral limit $`sup(m_i)0`$), the anomaly–free WI (3.2) implies the existence of $`L^21`$ non–singlet Goldstone bosons, interpolated by the hermitian fields $`O_b=i\overline{q}\gamma _5T^bq`$. Similarly, in the case of the $`U(1)`$ axial symmetry, one finds that:

$$d^4xT^\mu J_{5,\mu }(x)i\overline{q}\gamma _5q(0)=2i\overline{q}q,$$
(3.3)
where $`J_{5,\mu }=\overline{q}\gamma _\mu \gamma _5q`$ is the $`U(1)`$ axial current. But this is not the whole story! One also derives the following WI:

$$d^4xT^\mu J_{5,\mu }(x)O_P(0)=2Li_{eff}(0),$$
(3.4)
where $`_{eff}`$ is the $`2L`$–fermion operator defined by Eq. (3.1), while the hermitian field $`O_P`$ is defined as: $`O_Pi[det(\overline{q}_{sR}q_{tL})det(\overline{q}_{sL}q_{tR})]`$. If the $`U(1)`$–breaking condensate survives across the chiral transition at $`T_{ch}`$, i.e., $`C_{U(1)}=_{eff}(0)0`$ for $`T>T_{ch}`$ (while $`\overline{q}q=0`$ for $`T>T_{ch}`$), then this WI implies the existence of a (would–be) Goldstone boson (in the large–$`N_c`$ limit) coming from this breaking and interpolated by the hermitian field $`O_P`$. Therefore, the $`U(1)_A`$ (would–be) Goldstone boson (i.e., the $`\eta ^{}`$) is an “exotic” $`2L`$–fermion state for $`T>T_{ch}`$.

4. The new chiral effective Lagrangian

The proposed scenario, in which the $`U(1)`$ axial symmetry is (effectively) restored at a temperature $`T_{U(1)}`$ greater than $`T_{ch}`$, can be consistently reproduced using an effective–Lagrangian model. This analysis was originally performed in Refs.  and is summarized here for the convenience of the reader.

It is well known that the low–energy dynamics of the pseudoscalar mesons, including the effects due to the anomaly and the $`q\overline{q}`$ chiral condensate, and expanding to the first order in the light quark masses, can be described, in the large–$`N_c`$ limit, by an effective Lagrangian  written in terms of the mesonic field $`U_{ij}\overline{q}_{jR}q_{iL}`$ (up to a multiplicative constant) and the topological charge density $`Q`$. We make the assumption that there is a $`U(1)`$–breaking condensate which stays different from zero across $`T_{ch}`$, up to $`T_{U(1)}>T_{ch}`$: the form of this condensate has been discussed in the previous section. We must now define a field variable $`X`$, associated with this new condensate, to be inserted in the chiral Lagrangian. The translation from the fundamental quark fields to the effective–Lagrangian meson fields is done as follows. The operators $`i\overline{q}\gamma _5q`$ and $`\overline{q}q`$ entering in the WI (3.3) are essentially equal to (up to a multiplicative constant) $`i(\mathrm{Tr}U\mathrm{Tr}U^{})`$ and $`\mathrm{Tr}U+\mathrm{Tr}U^{}`$ respectively. Similarly, the operators $`_{eff}det(\overline{q}_{sR}q_{tL})+det(\overline{q}_{sL}q_{tR})`$ and $`O_Pi[det(\overline{q}_{sR}q_{tL})det(\overline{q}_{sL}q_{tR})]`$ entering in the WI (3.4) can be put equal to (up to a multiplicative constant) $`X+X^{}`$ and $`i(XX^{})`$ respectively, where $`Xdet\left(\overline{q}_{sR}q_{tL}\right)`$ is the new field variable (up to a multiplicative constant), related to the new $`U(1)`$ chiral condensate, which must be inserted in the chiral effective Lagrangian. It was shown in Refs.  that the most simple effective Lagrangian, constructed with the fields $`U`$, $`X`$ and $`Q`$, is:

$`(U,U^{},X,X^{},Q)={\displaystyle \frac{1}{2}}\mathrm{Tr}(_\mu U^\mu U^{})+{\displaystyle \frac{1}{2}}_\mu X^\mu X^{}`$ (4.1)
$`V(U,U^{},X,X^{})+{\displaystyle \frac{1}{2}}iQ(x)\omega _1\mathrm{Tr}(\mathrm{ln}U\mathrm{ln}U^{})`$
$`+{\displaystyle \frac{1}{2}}iQ(x)(1\omega _1)(\mathrm{ln}X\mathrm{ln}X^{})+{\displaystyle \frac{1}{2A}}Q^2(x),`$
where the potential term $`V(U,U^{},X,X^{})`$ has the form:

$`V(U,U^{},X,X^{})={\displaystyle \frac{1}{4}}\lambda _\pi ^2\mathrm{Tr}[(U^{}U\rho _\pi 𝐈)^2]+{\displaystyle \frac{1}{4}}\lambda _X^2(X^{}X\rho _X)^2`$ (4.2)
$`{\displaystyle \frac{B_m}{2\sqrt{2}}}\mathrm{Tr}(MU+M^{}U^{}){\displaystyle \frac{c_1}{2\sqrt{2}}}[det(U)X^{}+det(U^{})X].`$
I is the identity matrix. $`M`$ represents the quark mass matrix, $`M=\mathrm{diag}(m_1,\mathrm{},m_L)`$, which enters in the QCD Lagrangian as $`\delta _{QCD}^{(mass)}=\overline{q}_RMq_L\overline{q}_LM^{}q_R`$, and $`A`$ is the topological susceptibility in the pure–YM theory. All the parameters appearing in the Lagrangian must be considered as functions of the physical temperature $`T`$. In particular, the parameters $`\rho _\pi `$ and $`\rho _X`$ are responsible for the behaviour of the theory respectively across the $`SU(L)SU(L)`$ and the $`U(1)`$ chiral phase transitions, as follows:

$`\rho _\pi (T<T_{ch})`$ $``$ $`{\displaystyle \frac{1}{2}}F_\pi ^2>0,\rho _\pi (T>T_{ch})<0;`$
$`\rho _X(T<T_{U(1)})`$ $``$ $`{\displaystyle \frac{1}{2}}F_X^2>0,\rho _X(T>T_{U(1)})<0.`$ (4.3)
The parameter $`F_\pi `$ is the well–known pion decay constant, while the parameter $`F_X`$ is related to the new $`U(1)`$ axial condensate and will be discussed at length in the rest of the paper. For $`T<T_{ch}`$, $`\rho _\pi >0`$ and therefore, by virtue of the form (4.2) of the potential,***As it has been stressed in Ref.  (see also Ref. ), the linear $`\sigma `$–type model (4.1)–(4.2) contains redundant scalar fields (we are only interested in the pseudoscalar Goldstone, or would–be Goldstone, bosons), which can be eliminated by taking the limits $`\lambda _\pi ^2,\lambda _X^2+\mathrm{}`$ in Eq. (4.2): so an expansion is performed not only in powers of the light quark masses $`m_i`$, but also in powers of $`1/\lambda _\pi ^2`$ and $`1/\lambda _X^2`$. one finds that $`U0`$, or, in other words $`\overline{q}q0`$ (being $`U_{ij}\overline{q}_{jR}q_{iL}`$, up to a multiplicative constant): i.e., the $`SU(L)SU(L)`$ chiral symmetry is broken. Instead, for $`T>T_{ch}`$, $`\rho _\pi <0`$ and then, always from Eq. (4.2), one has that $`U=0`$, i.e., $`\overline{q}q=0`$. The $`U(1)`$ chiral condensate stays different from zero also in the region of temperatures $`T_{ch}<T<T_{U(1)}`$, where, on the contrary, the $`SU(L)SU(L)`$ chiral symmetry is restored. In fact, for $`T<T_{U(1)}`$, $`\rho _X>0`$ and therefore, from Eq. (4.2), one finds that $`X0`$, or, in other words, $`C_{U(1)}=_{eff}0`$ (being $`Xdet\left(\overline{q}_{sR}q_{tL}\right)`$, up to a multiplicative constant). The $`U(1)`$ chiral symmetry is (effectively) restored above $`T_{U(1)}`$, where $`\rho _X<0`$ and then, from Eq. (4.2), $`X=0`$, i.e., $`C_{U(1)}=_{eff}=0`$. According to what we have said in the Introduction and in Section 2, we also assume that the topological susceptibility $`A(T)`$ of the pure–YM theory drops to zero at a temperature $`T_\chi `$ greater than $`T_{ch}`$ (but smaller than, or equal to, $`T_{U(1)}`$).

One can study the mass spectrum of the theory for $`T<T_{ch}`$ and $`T_{ch}<T<T_{U(1)}`$. First of all, let us see what happens for $`T<T_{ch}`$, where both the $`q\overline{q}`$ chiral condensate and the $`U(1)`$ chiral condensate are present. Integrating out the field variable $`Q`$ and taking only the quadratic part of the Lagrangian, one finds that, in the chiral limit $`sup(m_i)0`$, there are $`L^21`$ zero–mass states, which represent the $`L^21`$ Goldstone bosons coming from the breaking of the $`SU(L)SU(L)`$ chiral symmetry down to $`SU(L)_V`$. Then there are two singlet eigenstates with non–zero masses:

$`\eta ^{}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{F_\pi ^2+LF_X^2}}}(\sqrt{L}F_XS_X+F_\pi S_\pi ),`$
$`\eta _X`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{F_\pi ^2+LF_X^2}}}(F_\pi S_X+\sqrt{L}F_XS_\pi ),`$ (4.4)
where $`S_\pi `$ is the usual $`SU(L)`$–singlet meson field associated with $`U`$, while $`S_X`$ is the meson field associated with $`X`$ (see Refs.  and Eqs. (5.4) below). The field $`\eta ^{}`$ has a “light” mass, in the sense of the $`N_c\mathrm{}`$ limit, being

$$m_\eta ^{}^2=\frac{2LA}{F_\pi ^2+LF_X^2}=𝒪(\frac{1}{N_c}).$$
(4.5)
This mass is intimately related to the anomaly and they both vanish in the $`N_c\mathrm{}`$ limit. On the contrary, the field $`\eta _X`$ has a sort of “heavy hadronic” mass of order $`𝒪(N_c^0)`$ in the large–$`N_c`$ limit. We immediately see that, if we put $`F_X=0`$ in the above–written formulae (i.e., if we neglect the new $`U(1)`$ chiral condensate), then $`\eta ^{}=S_\pi `$ and $`m_\eta ^{}^2`$ reduces to $`\frac{2LA}{F_\pi ^2}`$, which is the “usual” $`\eta ^{}`$ mass in the chiral limit . Yet, in the general case $`F_X0`$, the two states which diagonalize the squared mass matrix are linear combinations of the “quark–antiquark” singlet field $`S_\pi `$ and of the “exotic” field $`S_X`$. Both the $`\eta ^{}`$ and the $`\eta _X`$ have the same quantum numbers (spin, parity and so on), but they have a different quark content: one is mostly $`i(\overline{q}_Lq_R\overline{q}_Rq_L)`$, while the other is mostly $`i[det(\overline{q}_{sL}q_{tR})det(\overline{q}_{sR}q_{tL})]`$. What happens when approaching the chiral transition temperature $`T_{ch}`$? We know that $`F_\pi (T)0`$ when $`TT_{ch}`$. From Eq. (4.5) we see that $`m_\eta ^{}^2(T_{ch})=\frac{2A}{F_X^2}`$ and, from the first Eq. (4.4), $`\eta ^{}(T_{ch})=S_X`$. We have continuity in the mass spectrum of the theory through the chiral phase transition at $`T=T_{ch}`$. In fact, if we study the mass spectrum of the theory in the region of temperatures $`T_{ch}<T<T_{U(1)}`$ (where the $`SU(L)SU(L)`$ chiral symmetry is restored, while the $`U(1)`$ chiral condensate is still present), one finds that there is a singlet meson field $`S_X`$ (associated with the field $`X`$ in the chiral Lagrangian) with a squared mass given by (in the chiral limit): $`m_{S_X}^2=\frac{2A}{F_X^2}`$. This is nothing but the would–be Goldstone particle coming from the breaking of the $`U(1)`$ chiral symmetry, i.e., the $`\eta ^{}`$, which, for $`T>T_{ch}`$, is a sort of “exotic” matter field of the form $`i[det(\overline{q}_{sL}q_{tR})det(\overline{q}_{sR}q_{tL})]`$. Its existence could be proved perhaps in the near future by heavy–ion experiments.

5. A relation between $`\chi ^{}`$ and the new $`U(1)`$ chiral condensate

In this section and in the following one we want to describe some methods which provide us with some information about the parameter $`F_X`$. This quantity is a $`U(1)`$–breaking parameter: indeed, from Eq. (4.3), $`\rho _X=\frac{1}{2}F_X^2>0`$ for $`T<T_{U(1)}`$, and therefore, from Eq. (4.2), $`X=F_X/\sqrt{2}0`$. Remembering that $`Xdet\left(\overline{q}_{sR}q_{tL}\right)`$, up to a multiplicative constant, we find that $`F_X`$ is proportional to the new $`2L`$–fermion condensate $`C_{U(1)}=_{eff}`$ introduced above.
In the same way, the pion decay constant $`F_\pi `$, which controls the breaking of the $`SU(L)SU(L)`$ symmetry, is related to the $`q\overline{q}`$ chiral condensate by a simple and well–known proportionality relation (see Ref.  and references therein): $`\overline{q}_iq_i_{T<T_{ch}}\frac{1}{2}B_mF_\pi `$. Considering, for simplicity, the case of $`L`$ light quarks with the same mass $`m`$, one immediately derives from this equation the so–called *Gell-Mann–Oakes–Renner relation* :

$$m_{NS}^2F_\pi ^2\frac{2m}{L}\overline{q}q_{T<T_{ch}},$$
(5.1)
where, as usual, $`\overline{q}q_{i=1}^L\overline{q}_iq_i`$, and, moreover, $`m_{NS}^2=mB_m/F_\pi `$, $`m_{NS}`$ being the mass of the non–singlet pseudoscalar mesons. Eq. (5.1) relates, on the left–hand side, the pion decay constant $`F_\pi `$ and the mass $`m_{NS}`$ of the non–singlet mesons with, on the right–hand side, the chiral condensate $`\overline{q}q`$ and the quark mass $`m`$.
It is not possible to find, in a simple way, the analogous relation between $`F_X`$ and the new condensate $`C_{U(1)}=_{eff}`$, since the QCD Lagrangian does not contain any term proportional to the $`2L`$–fermion operator $`_{eff}`$.

Alternatively, the quantity $`F_X`$ can be written in terms of a certain two–point Green function of the topological charge–density operator $`Q(x)`$ in the full theory with quarks.
If we want to derive the two–point function of $`Q(x)`$, we need to consider the effective Lagrangian in the form (4.1), where the field variable $`Q(x)`$ has not yet been integrated. Therefore:

$$\chi (k)id^4xe^{ikx}TQ(x)Q(0)=(𝒦^1(k))_{11},$$
(5.2)
where $`𝒦^1(k)`$ is the inverse of the matrix $`𝒦(k)`$ associated with the quadratic part of the Lagrangian (4.1) in the momentum space, for the ensemble of fields $`(Q(x),\mathrm{})`$.
In particular, for $`T<T_{ch}`$, one has to consider the following quadratic Lagrangian, in the chiral limit $`sup(m_i)0`$:

$`_2={\displaystyle \frac{1}{2}}{\displaystyle \underset{a=1}{\overset{L^21}{}}}_\mu \pi _a^\mu \pi _a+{\displaystyle \frac{1}{2}}_\mu S_\pi ^\mu S_\pi +{\displaystyle \frac{1}{2}}_\mu S_X^\mu S_X`$ (5.3)
$`{\displaystyle \frac{1}{2}}c\left({\displaystyle \frac{\sqrt{2L}}{F_\pi }}S_\pi {\displaystyle \frac{\sqrt{2}}{F_X}}S_X\right)^2+{\displaystyle \frac{1}{2A}}Q^2`$
$`\omega _1{\displaystyle \frac{\sqrt{2L}}{F_\pi }}S_\pi Q(1\omega _1){\displaystyle \frac{\sqrt{2}}{F_X}}S_XQ,`$
where $`c\frac{c_1}{\sqrt{2}}\left(\frac{F_X}{\sqrt{2}}\right)\left(\frac{F_\pi }{\sqrt{2}}\right)^L`$. We have used for $`U`$ and $`X`$ the following exponential form, valid for $`T<T_{ch}`$ :

$`U`$ $`=`$ $`{\displaystyle \frac{F_\pi }{\sqrt{2}}}\mathrm{exp}\left({\displaystyle \frac{i\sqrt{2}}{F_\pi }}\mathrm{\Phi }\right),\mathrm{with}:\mathrm{\Phi }={\displaystyle \underset{a=1}{\overset{L^21}{}}}\pi _a\tau _a+{\displaystyle \frac{S_\pi }{\sqrt{L}}}\text{I};`$
$`X`$ $`=`$ $`{\displaystyle \frac{F_X}{\sqrt{2}}}\mathrm{exp}\left({\displaystyle \frac{i\sqrt{2}}{F_X}}S_X\right),`$ (5.4)
where the matrices $`\tau _a`$ ($`a=1,\mathrm{},L^21`$) are the generators of the algebra of $`SU(L)`$ in the fundamental representation, with normalization: $`\mathrm{Tr}(\tau _a)=0,\mathrm{Tr}(\tau _a\tau _b)=\delta _{ab}.`$
As already pointed out in Section 4, $`S_\pi `$ is the usual “quark–antiquark” $`SU(L)`$–singlet meson field associated with $`U`$, while $`S_X`$ is the “exotic” $`2L`$–fermion meson field associated with $`X`$. The $`\pi _a`$ ($`a=1,\mathrm{},L^21`$) are the self–hermitian fields describing the $`L^21`$ pions: they are massless in the chiral limit $`sup(m_i)0`$.
From the quadratic part of the Lagrangian (5.3) in the momentum space, we derive the following matrix $`𝒦(k)`$ for the ensemble of fields $`(Q,S_X,S_\pi )`$ \[the contribution of the pion fields $`\pi _a`$ is simply diagonal, diag($`k^2,\mathrm{},k^2`$), and therefore can be trivially factorized out\]:

$$𝒦(k)=\left(\begin{array}{ccc}\frac{1}{A}& \frac{\sqrt{2}(1\omega _1)}{F_X}& \frac{\omega _1\sqrt{2L}}{F_\pi }\\ & & \\ \frac{\sqrt{2}(1\omega _1)}{F_X}& k^2\frac{2c}{F_X^2}& \frac{2c\sqrt{L}}{F_\pi F_X}\\ & & \\ \frac{\omega _1\sqrt{2L}}{F_\pi }& \frac{2c\sqrt{L}}{F_\pi F_X}& k^2\frac{2Lc}{F_\pi ^2}\end{array}\right).$$
(5.5)
The calculation of the right–hand side of Eq. (5.2) can then be performed explicitly, using Eq. (5.5), obtaining:

$`\chi (k)={\displaystyle \frac{A\left[k^4\frac{2c(F_\pi ^2+LF_X^2)}{F_\pi ^2F_X^2}k^2\right]}{\left[k^4\frac{2c(F_\pi ^2+LF_X^2)}{F_\pi ^2F_X^2}k^2\right]2A\left[\frac{(1\omega _1)^2F_\pi ^2+L\omega _1^2F_X^2}{F_\pi ^2F_X^2}\right]k^2+\frac{4LAc}{F_\pi ^2F_X^2}}}.`$ (5.6)
One immediately sees that $`\chi \chi (0)=0`$, i.e., the topological susceptibility in the full theory with quarks vanishes in the chiral limit $`sup(m_i)0`$, as expected.
But the most interesting result is found when considering the so–called “slope” of the topological susceptibility, defined as:

$`\chi ^{}{\displaystyle \frac{1}{8}}{\displaystyle \frac{}{k_\mu }}{\displaystyle \frac{}{k^\mu }}\chi (k)|_{k=0}={\displaystyle \frac{i}{8}}{\displaystyle d^4xx^2TQ(x)Q(0)}.`$ (5.7)
Moreover, whenever $`\chi (k)`$ is a function of $`k^2`$ (such as in the theory at $`T=0`$, by virtue of the Lorentz invariance), one can also write:

$`\chi ^{}={\displaystyle \frac{d}{dk^2}}\chi (k)|_{k=0}.`$ (5.8)
In the theory at finite temperature $`T0`$ the Lorentz invariance is broken down to the $`O(3)`$ invariance under spatial rotations only, so that $`\chi (k)`$ is, in general, a function of $`\stackrel{}{k}^2`$ and $`k^0`$. However, the propagator $`\chi (k)`$ obtained in Eq. (5.6) is a function of $`k^2=(k^0)^2\stackrel{}{k}^2`$ also in the case of non–zero temperature. Therefore, we can explicitly calculate the quantity $`\chi ^{}`$ in our effective model in the chiral limit $`sup(m_i)0`$ (and we shall call “$`\chi _{ch}^{}`$” its value), using Eqs. (5.6) and (5.8), obtaining:

$`\chi _{ch}^{}={\displaystyle \frac{1}{2L}}(F_\pi ^2+LF_X^2){\displaystyle \frac{1}{2L}}F_\eta ^{}^2,`$ (5.9)
where $`F_\eta ^{}\sqrt{F_\pi ^2+LF_X^2}`$ is the decay constant of the $`\eta ^{}`$ (at the leading order in the $`1/N_c`$ expansion), modified by the presence of the new $`U(1)`$ chiral order parameter .

For the benefit of the reader, we here briefly repeat the arguments developed in Ref. , leading to the result $`F_\eta ^{}\sqrt{F_\pi ^2+LF_X^2}`$: this will also provide us with an alternative derivation of Eq. (5.9). We shall consider the $`T=0`$ case for simplicity.
It turns out that $`\eta ^{}`$ is just the meson state, with a squared mass of order $`1/N_c`$, whose contribution to the full topological susceptibility $`\chi `$ exactly cancels out (in the chiral limit of massless quarks) the pure–gauge part $`A`$ of $`\chi `$, so making $`\chi =0`$: this is the so–called Witten’s mechanism. To see how this picture comes out in our theory, one first determines the $`U(1)`$ axial current, starting from our effective Lagrangian. This is easily done remembering how the fields $`U`$ and $`X`$ transform under a $`U(1)`$ chiral transformation and one ends up with the following expression :

$$J_{5,\mu }=i[\mathrm{Tr}(U^{}_\mu UU_\mu U^{})+L(X^{}_\mu XX_\mu X^{})].$$
(5.10)
After having inserted here the expressions (5.4) in place of $`U`$ and $`X`$, the current $`J_{5,\mu }`$ takes the following form:

$$J_{5,\mu }=\sqrt{2L}F_\eta ^{}_\mu \eta ^{},$$
(5.11)
where the field $`\eta ^{}`$ is defined by the first Eq. (4.4) and the relative coupling between $`J_{5,\mu }`$ and $`\eta ^{}`$, i.e., the $`SU(L)`$–singlet ($`\eta ^{}`$) decay constant defined as $`0|J_{5,\mu }(0)|\eta ^{}(p)=i\sqrt{2L}p_\mu F_\eta ^{}`$, is given by:

$$F_\eta ^{}=\sqrt{F_\pi ^2+LF_X^2}.$$
(5.12)
Let us now recall the Witten’s argument and write the two–point function (at four–momentum $`k`$) of the topological charge density $`Q(x)`$ as a sum over one–hadron poles, i.e., one–hadron intermediate states:

$$\chi (k)=id^4xe^{ikx}TQ(x)Q(0)=A_0(k)+\underset{mesons}{}\frac{|0|Q|n|^2}{k^2m_n^2},$$
(5.13)
where $`A_0(k)`$ is the pure Yang–Mills contribution from the glueball intermediate states and it is the leading–order term in $`1/N_c`$ \[being of order $`𝒪(N_c^0)`$\]. In the chiral limit in which we have $`L`$ massless quarks, the full topological susceptibility $`\chi \chi (k=0)`$ must vanish: so there must be a meson state, with squared mass $`m_n^2=𝒪(1/N_c)`$ \[since $`A_0(0)=𝒪(N_c^0)`$, while $`|0|Q|n|^2=𝒪(1/N_c)`$\], which exactly cancels out $`A_0(0)`$. This is the meson that we usually call $`\eta ^{}`$: the other meson states in (5.13) have squared masses of order $`𝒪(N_c^0)`$, so that their contributions to the summation in (5.13) are suppressed by a factor of $`1/N_c`$. Therefore we obtain that:

$$\frac{|0|Q|\eta ^{}|^2}{m_\eta ^{}^2}=A,$$
(5.14)
where $`AA_0(0)`$ is the pure Yang–Mills topological susceptibility in the large–$`N_c`$ limit. In the chiral limit of $`L`$ massless quarks, the topological charge density $`Q(x)`$ is directly related to the four–divergence of the axial current $`J_{5,\mu }`$, via the anomaly equation, $`^\mu J_{5,\mu }(x)=2LQ(x)`$, so that $`0|Q|\eta ^{}=\frac{1}{\sqrt{2L}}m_\eta ^{}^2F_\eta ^{}`$, which can be substituted into Eq. (5.14) to give:

$$A=\frac{m_\eta ^{}^2F_\eta ^{}^2}{2L}.$$
(5.15)
This equation relates the mass $`m_\eta ^{}`$ of the $`\eta ^{}`$ state, its decay constant $`F_\eta ^{}`$ and the pure–gauge topological susceptibility $`A`$. And in fact Eq. (5.15) is verified when putting for $`m_\eta ^{}`$ and $`F_\eta ^{}`$ their values determined above in Eqs. (4.5) and (5.12).
The dominance of the $`\eta ^{}`$ state in the sum over the meson states at the right–hand side of Eq. (5.13) can also be used to evaluate the slope of the topological susceptibility:

$`\chi _{ch}^{}={\displaystyle \frac{d}{dk^2}}\chi (k)|_{k=0}{\displaystyle \frac{|0|Q(0)|\eta ^{}|^2}{m_\eta ^{}^4}}={\displaystyle \frac{1}{2L}}F_\eta ^{}^2={\displaystyle \frac{1}{2L}}(F_\pi ^2+LF_X^2),`$ (5.16)
thus recovering the result (5.9), derived above from our effective Lagrangian. This is perfectly natural: using the effective Lagrangian at tree–level (i.e., using its “free” propagators) one gets the results in the one–hadron pole (i.e., one–hadron intermediate state) approximation.

Summarizing, we have found that the value of $`\chi ^{}`$ in the chiral limit $`sup(m_i)0`$ is shifted from the “original” value $`\frac{1}{2L}F_\pi ^2`$ (derived in the absence of an extra $`U(1)`$ chiral condensate: see Refs. , where $`\chi _{ch}^{}`$ is shown to be a relevant quantity in the discussion of the so–called “proton–spin crisis” problem, and also Ref. ) to the value $`\frac{1}{2L}F_\eta ^{}^2=\frac{1}{2L}(F_\pi ^2+LF_X^2)`$, which also depends on the quantity $`F_X`$, proportional to the extra $`U(1)`$ chiral condensate.
Therefore, a measure of this quantity $`\chi _{ch}^{}`$, e.g., in lattice gauge theory, could provide an estimate for the $`\eta ^{}`$ decay constant $`F_\eta ^{}`$, and, as a consequence, for $`F_X`$. At present, lattice determinations of $`\chi ^{}`$ only exist for the pure–gauge theory at $`T=0`$, with gauge group $`SU(2)`$  and $`SU(3)`$ .***It has been recently pointed out in Ref.  that there can be ambiguities in the definition of $`\chi _{ch}^{}`$ in a lattice regularized theory. However, these ambiguities do not apply to the leading term \[of order $`𝒪(N_c)`$\] given by Eq. (5.9), but they can affect the next–to–leading term \[of order $`𝒪(N_c^0)`$\] in the $`1/N_c`$ expansion.

All the above refers to the theory at $`T=0`$ (or, more generally, for $`T<T_{ch}`$).
When approaching the chiral transition at $`T=T_{ch}`$, one expects that $`F_\pi `$ vanishes, while $`F_X`$ remains different from zero and the quantity $`\chi _{ch}^{}`$ tends to the value:

$$\chi _{ch}^{}\underset{TT_{ch}}{}\frac{1}{2}F_X^2.$$
(5.17)
Indeed, the quantity $`\chi (k)=id^4xe^{ikx}TQ(x)Q(0)`$ can also be evaluated in the region of temperatures $`T_{ch}<T<T_{U(1)}`$, proceeding as for the case $`T<T_{ch}`$, obtaining the result (already derived in Ref. ):

$$\chi (k)=A\frac{k^2}{k^2\frac{2A}{F_X^2}},$$
(5.18)
in the chiral limit $`sup(m_i)0`$.
Therefore, in the region of temperatures $`T_{ch}<T<T_{U(1)}`$, $`\chi _{ch}^{}`$ is given by:

$$\chi _{ch}^{}=\frac{d}{dk^2}\chi (k)|_{k=0}=\frac{1}{2}F_X^2,$$
(5.19)
consistently with the results (5.9) and (5.17) found above: i.e., $`\chi _{ch}^{}`$ varies with continuity across $`T_{ch}`$. This means that $`\chi _{ch}^{}`$ acts as a sort of order parameter for the $`U(1)`$ axial symmetry above $`T_{ch}`$: if $`\chi _{ch}^{}`$ is different from zero above $`T_{ch}`$, this means that the $`U(1)`$–breaking parameter $`F_X`$ is different from zero.

6. Radiative decays of the pseudoscalar mesons

Further information on the quantity $`F_X`$ (i.e., on the new $`U(1)`$ chiral condensate, to which it is related) can be derived from the study of the radiative decays of the “light” pseudoscalar mesons in two photons. These decays were also studied in Ref. , using an effective Lagrangian model, in which only the $`q\overline{q}`$ chiral condensate was considered. In this section we want to find the decay rates of the processes $`\pi ^0,\eta ,\eta ^{},\eta _X\gamma \gamma `$ (where $`\eta _X`$ represents the pseudoscalar meson state, introduced in Section 4, having the same quantum numbers of the $`\eta ^{}`$, but a larger mass and a different quark content) and see which are the effects due to the new $`U(1)`$ axial condensate, in the realistic case of $`L=3`$ light quarks and in the simple case of zero temperature ($`T=0`$).
To this purpose, we have to introduce the electromagnetic interaction in our effective model. First of all, in order to make the Lagrangian invariant under local $`U(1)`$ electromagnetic transformations \[$`qq^{}=e^{i\theta e𝐐}q`$, in terms of the quark fields; the matrix $`𝐐`$ is defined in Eq. (6.2) below\], we have to replace the derivative of the field $`U`$ with the corresponding covariant derivative $`D_\mu U`$, which, by virtue of the transformation property of the field $`U`$ \[$`UU^{}=e^{i\theta e𝐐}Ue^{i\theta e𝐐}`$\], has the following form:

$$D_\mu U=_\mu U+ieA_\mu [𝐐,U],$$
(6.1)
where $`A_\mu `$ is the electromagnetic field (which transforms as: $`A_\mu A_\mu ^{}=A_\mu _\mu \theta `$) and $`𝐐`$ is the quark charge matrix (in units of $`e`$, the absolute value of the electron charge):

$$𝐐=\left(\begin{array}{ccc}\frac{2}{3}& & \\ & \frac{1}{3}& \\ & & \frac{1}{3}\end{array}\right).$$
(6.2)
Instead,the field $`X`$ is invariant under a $`U(1)`$ electromagnetic gauge transformation and therefore its covariant derivative just coincides with the ordinary four–derivative: $`D_\mu X=_\mu X`$.
In addition, we have to reproduce the effects of the electromagnetic anomaly, whose contribution to the four–divergence of the $`U(1)`$ axial current ($`J_{5,\mu }=\overline{q}\gamma _\mu \gamma _5q`$) and of the $`SU(3)`$ axial currents ($`A_\mu ^a=\overline{q}\gamma _\mu \gamma _5\frac{\tau _a}{\sqrt{2}}q`$, where the matrices $`\tau _a`$, with $`a=1,\mathrm{},8`$, are the generators of the algebra of $`SU(3)`$ in the fundamental representation, already introduced in the previous section) is given by:

$$(^\mu J_{5,\mu })_{anomaly}^{e.m.}=2\mathrm{T}\mathrm{r}(𝐐^2)G,(^\mu A_\mu ^a)_{anomaly}^{e.m.}=2\mathrm{T}\mathrm{r}\left(𝐐^2\frac{\tau _a}{\sqrt{2}}\right)G,$$
(6.3)
where $`G\frac{e^2N_c}{32\pi ^2}\epsilon ^{\mu \nu \rho \sigma }F_{\mu \nu }F_{\rho \sigma }`$ ($`F_{\mu \nu }`$ being the electromagnetic field–strength tensor), thus breaking the corresponding chiral symmetries. We observe that $`\mathrm{Tr}(𝐐^2\tau _a)0`$ only for $`a=3`$ or $`a=8`$.
We must look for an interaction term $`_I`$ (constructed with the chiral Lagrangian fields and the electromagnetic operator $`G`$) which, under a $`U(1)`$ axial transformation $`qq^{}=e^{i\alpha \gamma _5}q`$, transforms as:

$$U(1)_A:_I_I+2\alpha \mathrm{Tr}(𝐐^2)G,$$
(6.4)
while, under $`SU(3)`$ axial transformations of the type $`qq^{}=e^{i\beta \gamma _5\tau _a/\sqrt{2}}q`$ (with $`a=3,8`$), transforms as:

$$SU(3)_A:_I_I+2\beta \mathrm{Tr}\left(𝐐^2\frac{\tau _a}{\sqrt{2}}\right)G.$$
(6.5)
The electromagnetic anomaly term for the field $`U`$, originally proposed in Ref. , has the following form: $`\frac{1}{2}iG\mathrm{Tr}[𝐐^2(\mathrm{ln}U\mathrm{ln}U^{})]`$. By virtue of the transformation property of the field $`U`$ under a $`U(L)U(L)`$ chiral transformation ($`q_LV_Lq_L`$, $`q_RV_Rq_R`$ $``$ $`UV_LUV_R^{}`$, where $`V_L`$ and $`V_R`$ are arbitrary unitary matrices ), it is immediate to see that this term satisfies both the transformation properties (6.4) and (6.5).
If we also consider an analogous term for the field $`X`$, of the form $`\frac{1}{2L}iG\mathrm{Tr}(𝐐^2)(\mathrm{ln}X\mathrm{ln}X^{})`$ (where, in our case, $`L=3`$), then, by virtue of the transformation property of the field $`X`$ under a $`U(L)U(L)`$ chiral transformation ($`Xdet(V_L)det(V_R)^{}X`$ ), one can see that this term satisfies the transformation property (6.4) under a $`U(1)`$ axial transformation, while it is invariant under $`SU(3)`$ axial transformations. Therefore, if we try to consider a linear combination of the two terms (inspired by what was done for the $`U(1)`$ axial anomaly terms in Eq. (4.1) ), i.e.,

$$_I=\frac{1}{2}i\omega _2G\mathrm{Tr}[𝐐^2(\mathrm{ln}U\mathrm{ln}U^{})]+\frac{1}{2L}i(1\omega _2)G\mathrm{Tr}(𝐐^2)(\mathrm{ln}X\mathrm{ln}X^{}),$$
(6.6)
we immediately see that the property (6.4) is satisfied for every value of the parameter $`\omega _2`$, while the property (6.5) is satisfied only for $`\omega _2=1`$.
In conclusion, the electromagnetic anomaly interaction term is simply given by:

$$_I=\frac{1}{2}iG\mathrm{Tr}[𝐐^2(\mathrm{ln}U\mathrm{ln}U^{})].$$
(6.7)
Therefore, we have to consider the following effective chiral Lagrangian, which includes the electromagnetic interaction terms described above:

$`(U,U^{},X,X^{},Q,A^\mu )={\displaystyle \frac{1}{2}}\mathrm{Tr}(D_\mu UD^\mu U^{})+{\displaystyle \frac{1}{2}}_\mu X^\mu X^{}`$ (6.8)
$`V(U,U^{},X,X^{})+{\displaystyle \frac{1}{2}}iQ\omega _1\mathrm{Tr}(\mathrm{ln}U\mathrm{ln}U^{})`$
$`+{\displaystyle \frac{1}{2}}iQ(1\omega _1)(\mathrm{ln}X\mathrm{ln}X^{})+{\displaystyle \frac{1}{2A}}Q^2`$
$`+{\displaystyle \frac{1}{2}}iG\mathrm{Tr}[𝐐^2(\mathrm{ln}U\mathrm{ln}U^{})]{\displaystyle \frac{1}{4}}F_{\mu \nu }F^{\mu \nu },`$
where the potential term $`V(U,U^{},X,X^{})`$ is the one written in Eq. (4.2).

We can now describe the predictions of our model about the decay amplitudes and rates of the processes $`\pi ^0,\eta ,\eta ^{},\eta _X\gamma \gamma `$. It is well known that the two final photons in the pseudoscalar–meson decays are associated with the pseudoscalar operator $`\epsilon ^{\mu \nu \rho \sigma }F_{\mu \nu }F_{\rho \sigma }`$: so, e.g., the $`\pi ^0`$ decay is reproduced by an interaction of the form $`\pi _3\epsilon ^{\mu \nu \rho \sigma }F_{\mu \nu }F_{\rho \sigma }`$. The covariant derivative does not contain this term (since the diagonal matrix $`𝐐`$ and the generator $`\tau _3`$ commute); instead, it is produced by the electromagnetic anomaly interaction written in Eq. (6.7) (indeed the chiral symmetry is broken by the anomaly). The same arguments can be applied to the other decays. Therefore, the decay amplitude of the generic process “$`meson\gamma \gamma `$” is given by:

$$A(meson\gamma \gamma )=\gamma \gamma |_I(0)|meson.$$
(6.9)
Substituting the exponential form (5.4) (with $`L=3`$) of the field $`U`$ (we remind the reader that we are considering the case of zero temperature) into Eq. (6.7), we have the following interaction Lagrangian:

$$_I=G\frac{1}{3F_\pi }\left(\pi _3+\frac{1}{\sqrt{3}}\pi _8+\frac{2\sqrt{2}}{\sqrt{3}}S_\pi \right).$$
(6.10)
We immediately consider the realistic case of $`L=3`$ quarks with masses different from zero. In this case one has to consider the following quadratic Lagrangian, which can be obtained substituting the exponential expressions (5.4) of the fields $`U`$ and $`X`$ into Eq. (4.1):

$`_2={\displaystyle \frac{1}{2}}\mathrm{Tr}(_\mu \mathrm{\Phi }^\mu \mathrm{\Phi })+{\displaystyle \frac{1}{2}}_\mu S_X^\mu S_X{\displaystyle \frac{B_m}{2F_\pi }}\mathrm{Tr}(M\mathrm{\Phi }^2)`$ (6.11)
$`c\left({\displaystyle \frac{1}{F_\pi }}\mathrm{Tr}\mathrm{\Phi }{\displaystyle \frac{1}{F_X}}S_X\right)^2A\left({\displaystyle \frac{\omega _1}{F_\pi }}\mathrm{Tr}\mathrm{\Phi }+{\displaystyle \frac{1\omega _1}{F_x}}S_X\right)^2,`$
where $`c\frac{c_1}{\sqrt{2}}\left(\frac{F_X}{\sqrt{2}}\right)\left(\frac{F_\pi }{\sqrt{2}}\right)^3`$.
Substituting Eq. (5.4) (with $`L=3`$) into Eq. (6.11), one immediately sees that the fields $`\pi _1,\pi _2,\pi _4,\pi _5,\pi _6,\pi _7`$ are diagonal, while the fields $`\pi _3,\pi _8,S_\pi ,S_X`$ mix together. However, neglecting the experimentally small mass difference between the quarks *up* and *down* (i.e., neglecting the experimentally small violations of the $`SU(2)`$ isotopic spin), also $`\pi _3`$ becomes diagonal and can be identified with the physical state $`\pi ^0`$. The other physical states can be found by the diagonalization of the following squared mass matrix, written for the ensemble of fields $`(\pi _8,S_\pi ,S_X)`$ (originally derived in Ref. ):

$$𝒦=\left(\begin{array}{ccc}\frac{2B}{3}(\stackrel{~}{m}+2m_s)& \frac{2B\sqrt{2}}{3}(\stackrel{~}{m}m_s)& 0\\ & & \\ \frac{2B\sqrt{2}}{3}(\stackrel{~}{m}m_s)& \frac{6(A\omega _1^2+c)}{F_\pi ^2}+\stackrel{~}{m}_0^2& \frac{2\sqrt{3}[A\omega _1(1\omega _1)c]}{F_\pi F_X}\\ & & \\ 0& \frac{2\sqrt{3}[A\omega _1(1\omega _1)c]}{F_\pi F_X}& \frac{2[A(1\omega _1)^2+c]}{F_X^2}\end{array}\right),$$
(6.12)
where:

$$B\frac{B_m}{2F_\pi },\stackrel{~}{m}\frac{m_u+m_d}{2},\stackrel{~}{m}_0^2\frac{2}{3}B(2\stackrel{~}{m}+m_s).$$
(6.13)
The fields $`(\pi _8,S_\pi ,S_X)`$ can be written in terms of the eigenstates $`\eta `$, $`\eta ^{}`$, $`\eta _X`$ as follows:

$$\left(\begin{array}{c}\pi _8\\ S_\pi \\ S_X\end{array}\right)=𝐂\left(\begin{array}{c}\eta \\ \eta ^{}\\ \eta _X\end{array}\right),$$
(6.14)
where $`𝐂`$ is the following $`3\times 3`$ orthogonal matrix:

$$𝐂=\left(\begin{array}{ccc}\alpha _1& \alpha _2& \alpha _3\\ \beta _1& \beta _2& \beta _3\\ \gamma _1& \gamma _2& \gamma _3\end{array}\right)=\left(\begin{array}{ccc}\mathrm{cos}\stackrel{~}{\phi }& \mathrm{sin}\stackrel{~}{\phi }& 0\\ & & \\ \mathrm{sin}\stackrel{~}{\phi }\frac{F_\pi }{\sqrt{F_\pi ^2+3F_X^2}}& \mathrm{cos}\stackrel{~}{\phi }\frac{F_\pi }{\sqrt{F_\pi ^2+3F_X^2}}& \frac{\sqrt{3}F_X}{\sqrt{F_\pi ^2+3F_X^2}}\\ & & \\ \mathrm{sin}\stackrel{~}{\phi }\frac{\sqrt{3}F_X}{\sqrt{F_\pi ^2+3F_X^2}}& \mathrm{cos}\stackrel{~}{\phi }\frac{\sqrt{3}F_X}{\sqrt{F_\pi ^2+3F_X^2}}& \frac{F_\pi }{\sqrt{F_\pi ^2+3F_X^2}}\end{array}\right).$$
(6.15)
Here $`\stackrel{~}{\phi }`$ is a mixing angle, which can be related to the masses of the quarks *up, down, strange* (and therefore to the masses of the octet mesons) by the following relation:

$`\mathrm{tan}\stackrel{~}{\phi }={\displaystyle \frac{\sqrt{2}}{9A}}BF_\pi \sqrt{F_\pi ^2+3F_X^2}(m_s\stackrel{~}{m})={\displaystyle \frac{F_\pi \sqrt{F_\pi ^2+3F_X^2}}{6\sqrt{2}A}}(m_\eta ^2m_\pi ^2),`$ (6.16)
where $`m_\pi ^2=2B\stackrel{~}{m}`$ and $`m_\eta ^2=\frac{2}{3}B(\stackrel{~}{m}+2m_s)`$. The masses $`m_\eta `$, $`m_\eta ^{}`$, $`m_{\eta _X}`$ are given by the squared root of the eigenvalues of the squared mass matrix written in Eq. (6.12) .

The interaction Lagrangian (6.10), written in terms of the physical fields $`\pi ^0,\eta ,\eta ^{}`$ and $`\eta _X`$, reads as follows:

$$_IG\frac{1}{3F_\pi }\left(\pi ^0+a_1\eta +a_2\eta ^{}+a_3\eta _X\right),$$
(6.17)
where $`a_i=\frac{1}{\sqrt{3}}(\alpha _i+2\sqrt{2}\beta _i)`$ $`(\mathrm{for}i=1,2,3)`$, so that:

$`a_1`$ $`=`$ $`\sqrt{{\displaystyle \frac{1}{3}}}\left(\mathrm{cos}\stackrel{~}{\phi }+2\sqrt{2}\mathrm{sin}\stackrel{~}{\phi }{\displaystyle \frac{F_\pi }{\sqrt{F_\pi ^2+3F_X^2}}}\right),`$ (6.18)
$`a_2`$ $`=`$ $`\sqrt{{\displaystyle \frac{1}{3}}}\left(2\sqrt{2}\mathrm{cos}\stackrel{~}{\phi }{\displaystyle \frac{F_\pi }{\sqrt{F_\pi ^2+3F_X^2}}}\mathrm{sin}\stackrel{~}{\phi }\right),`$ (6.19)
$`a_3`$ $`=`$ $`2\sqrt{2}\left({\displaystyle \frac{F_X}{\sqrt{F_\pi ^2+3F_X^2}}}\right).`$ (6.20)
With simple calculations we find the following expressions for the decay amplitudes:

$`A(\pi ^0\gamma \gamma )`$ $`=`$ $`{\displaystyle \frac{e^2N_c}{12\pi ^2F_\pi }}I,`$ (6.21)
$`A(\eta \gamma \gamma )`$ $`=`$ $`{\displaystyle \frac{e^2N_c}{12\pi ^2F_\pi }}\sqrt{{\displaystyle \frac{1}{3}}}\left(\mathrm{cos}\stackrel{~}{\phi }+2\sqrt{2}\mathrm{sin}\stackrel{~}{\phi }{\displaystyle \frac{F_\pi }{\sqrt{F_\pi ^2+3F_X^2}}}\right)I,`$ (6.22)
$`A(\eta ^{}\gamma \gamma )`$ $`=`$ $`{\displaystyle \frac{e^2N_c}{12\pi ^2F_\pi }}\sqrt{{\displaystyle \frac{1}{3}}}\left(2\sqrt{2}\mathrm{cos}\stackrel{~}{\phi }{\displaystyle \frac{F_\pi }{\sqrt{F_\pi ^2+3F_X^2}}}\mathrm{sin}\stackrel{~}{\phi }\right)I,`$ (6.23)
$`A(\eta _X\gamma \gamma )`$ $`=`$ $`{\displaystyle \frac{e^2N_c}{12\pi ^2F_\pi }}2\sqrt{2}\left({\displaystyle \frac{F_X}{\sqrt{F_\pi ^2+3F_X^2}}}\right)I,`$ (6.24)
where $`I\epsilon _{\mu \nu \rho \sigma }k_1^\mu ϵ_1^\nu k_2^\rho ϵ_2^\sigma `$ ($`k_1`$, $`k_2`$ being the four–momenta of the two final photons and $`ϵ_1`$, $`ϵ_2`$ their polarizations). Consequently we derive the following decay rates (in the real case $`N_c=3`$):

$`\mathrm{\Gamma }(\pi ^0\gamma \gamma )`$ $`=`$ $`{\displaystyle \frac{\alpha ^2m_\pi ^3}{64\pi ^3F_\pi ^2}},`$ (6.25)
$`\mathrm{\Gamma }(\eta \gamma \gamma )`$ $`=`$ $`{\displaystyle \frac{\alpha ^2m_\eta ^3}{192\pi ^3F_\pi ^2}}\left(\mathrm{cos}\stackrel{~}{\phi }+2\sqrt{2}\mathrm{sin}\stackrel{~}{\phi }{\displaystyle \frac{F_\pi }{\sqrt{F_\pi ^2+3F_X^2}}}\right)^2,`$ (6.26)
$`\mathrm{\Gamma }(\eta ^{}\gamma \gamma )`$ $`=`$ $`{\displaystyle \frac{\alpha ^2m_\eta ^{}^3}{192\pi ^3F_\pi ^2}}\left(2\sqrt{2}\mathrm{cos}\stackrel{~}{\phi }{\displaystyle \frac{F_\pi }{\sqrt{F_\pi ^2+3F_X^2}}}\mathrm{sin}\stackrel{~}{\phi }\right)^2,`$ (6.27)
$`\mathrm{\Gamma }(\eta _X\gamma \gamma )`$ $`=`$ $`{\displaystyle \frac{\alpha ^2m_{\eta _X}^3}{8\pi ^3F_\pi ^2}}\left({\displaystyle \frac{F_X}{\sqrt{F_\pi ^2+3F_X^2}}}\right)^2,`$ (6.28)
where $`\alpha =e^2/4\pi 1/137`$ is the fine–structure constant.

If we put $`F_X=0`$ (i.e., if we neglect the new $`U(1)`$ chiral condensate), the expressions (6.22)$`÷`$(6.23) written above reduce to the corresponding amplitudes derived in Ref.  using an effective Lagrangian which includes only the usual $`q\overline{q}`$ chiral condensate (so there is no field $`\eta _X`$!):

$`A(\eta \gamma \gamma )`$ $`=`$ $`{\displaystyle \frac{e^2N_c}{12\pi ^2F_\pi }}\sqrt{{\displaystyle \frac{1}{3}}}(\mathrm{cos}\phi +2\sqrt{2}\mathrm{sin}\phi )I,`$ (6.29)
$`A(\eta ^{}\gamma \gamma )`$ $`=`$ $`{\displaystyle \frac{e^2N_c}{12\pi ^2F_\pi }}\sqrt{{\displaystyle \frac{1}{3}}}(2\sqrt{2}\mathrm{cos}\phi \mathrm{sin}\phi )I.`$ (6.30)
The mixing angle $`\phi `$ is now defined as follows \[see Eqs. (6.14)–(6.15) with $`F_X=0`$\]:

$`\pi _8`$ $`=`$ $`\eta \mathrm{cos}\phi \eta ^{}\mathrm{sin}\phi ,`$
$`S_\pi `$ $`=`$ $`\eta \mathrm{sin}\phi +\eta ^{}\mathrm{cos}\phi ,`$ (6.31)
and is related to the quark masses (i.e., to the masses of the octet mesons) by the following relation \[see Eq. (6.16) with $`F_X=0`$\]:

$$\mathrm{tan}\phi \frac{\sqrt{2}}{9A}BF_\pi ^2(m_s\stackrel{~}{m})=\frac{F_\pi ^2}{6\sqrt{2}A}(m_\eta ^2m_\pi ^2).$$
(6.32)
Therefore, the introduction of the new condensate (while leaving the $`\pi ^0\gamma \gamma `$ decay rate unaffected, as it must!) modifies the decay rates of $`\eta `$ and $`\eta ^{}`$ (and, moreover, we also have to consider the particle $`\eta _X`$). In particular, it modifies the $`\eta ^{}`$ decay constant, already in the chiral limit $`sup(m_i)0`$. Indeed, in this limit, Eq. (6.30) \[with $`\phi =0`$: see Eq. (6.32)\] differs from the corresponding Eq. (6.23) \[with $`\stackrel{~}{\phi }=0`$: see Eq. (6.16)\] by the substitution of the pion decay constant $`F_\pi `$ with the quantity $`F_\eta ^{}\sqrt{F_\pi ^2+3F_X^2}`$, which (as it has been explained in the previous section for the general case of $`L`$ light flavours) can be identified with the $`\eta ^{}`$ decay constant (in the large–$`N_c`$ limit).

In conclusion, a study of the radiative decays $`\eta \gamma \gamma `$, $`\eta ^{}\gamma \gamma `$ and a comparison with the experimental data can provide us with further information about the parameter $`F_X`$ and the new exotic condensate. For example, from Eqs. (6.26) and (6.27), and using the experimental values for the various quantities which there appear, i.e.,

$`F_\pi =92.4(4)\mathrm{MeV},`$
$`m_\eta =547.30(12)\mathrm{MeV},`$
$`m_\eta ^{}=957.78(14)\mathrm{MeV},`$
$`\mathrm{\Gamma }(\eta \gamma \gamma )=0.46(4)\mathrm{KeV},`$
$`\mathrm{\Gamma }(\eta ^{}\gamma \gamma )=4.26(19)\mathrm{KeV},`$ (6.33)
we can extract the following values for the quantity $`F_X`$ and for the mixing angle $`\stackrel{~}{\phi }`$:

$$F_X=27(9)\mathrm{MeV},\stackrel{~}{\phi }=16(3)^0.$$
(6.34)
The value of $`F_X`$ is not far from the upper limit $`|F_X|20\mathrm{MeV}`$ obtained from the generalized Witten–Veneziano formula derived in Ref. . Moreover, the values of $`F_X`$ and $`\stackrel{~}{\phi }`$ so found are perfectly consistent with the relation (6.16) for the mixing angle.***We use for the pure–YM topological susceptibility the value $`A=(180\pm 5\mathrm{MeV})^4`$, obtained from lattice simulations . We thus see that $`F_X`$ is small when compared with $`F_\pi `$. In particular, the corrections due to the presence of the new $`U(1)`$ condensate are proportional to the ratio $`F_X^2/F_\pi ^2`$ \[see Eqs. (6.16), and (6.26)$`÷`$(6.28)\] and they are of the order of some $`\%`$.

7. Conclusions

In this paper we have discussed the role of the $`U(1)`$ axial symmetry in QCD both at zero and at finite temperature. One expects that, above a certain critical temperature $`T_{U(1)}`$, also the $`U(1)`$ axial symmetry will be (effectively) restored. We have tried to see if this transition has (or has not) anything to do with the usual $`SU(L)SU(L)`$ chiral transition: various possible scenarios have been discussed in Section 2. In particular, supported by recent lattice results on the pure–YM topological susceptibility and the so–called “chiral susceptibilities” (which have been discussed at length in the Introduction), we have analysed a scenario in which a new $`U(1)`$–breaking condensate survives across the chiral transition at $`T_{ch}`$, staying different from zero up to $`T_{U(1)}>T_{ch}`$. This scenario can be consistently reproduced using an effective Lagrangian model, which also includes the new $`U(1)`$ chiral condensate: this theoretical model was originally proposed in Refs.  and it has been briefly summarized in Sections 3 and 4 for the convenience of the reader. This scenario could perhaps be verified in the near future by heavy–ion experiments, by analysing the pseudoscalar–meson spectrum in the singlet sector.

In Section 5 (which, together with Section 6, contains the main original results of this paper) we have analysed the consequences of our theoretical model on the slope of the topological susceptibility $`\chi ^{}`$, in the full theory with quarks, showing how this quantity is modified by the presence of the new $`U(1)`$ chiral order parameter: we have found that $`\chi ^{}`$ (in the chiral limit $`sup(m_i)0`$) acts as an order parameter for the $`U(1)`$ axial symmetry above $`T_{ch}`$. This prediction of our model could be tested in the near–future Monte Carlo simulations on the lattice (but see also Ref. ).

Finally, in Section 6, we have found that the existence of the new $`U(1)`$ chiral condensate can be directly investigated by studying (at $`T=0`$) the radiative decays of the pseudoscalar mesons $`\eta `$ and $`\eta ^{}`$ in two photons. A first comparison of our results with the experimental data has been performed at the end of Section 6: the results are encouraging, pointing towards a certain evidence of a non–zero $`U(1)`$ axial condensate (i.e., $`F_X0`$). However, one should keep in mind that our results have been derived from a very simplified model, obtained doing a first–order expansion in $`1/N_c`$ and in the quark masses. We expect that such a model can furnish only qualitative or, at most, “semi–quantitative” predictions. Higher–order terms in $`1/N_c`$ could give rise to corrections (in the “real world” with $`N_c=3`$) of the same order of magnitude of (or even larger than!) those induced by the new $`U(1)`$ axial condensate (having the form $`F_X^2/F_\pi ^2`$). Moreover, when going beyond the leading order in $`1/N_c`$, it becomes necessary to take into account questions of renormalization–group behaviour of the various quantities and operators involved in our theoretical analysis. This issue has been widely discussed in the literature, both in relation to the analysis of $`\chi _{ch}^{}`$, in the context of the proton–spin crisis problem , and also in relation to the study of the $`\eta ^{}(\eta )`$ radiative decays . Further studies are therefore necessary in order to continue this analysis from a more quantitative point of view.
Last, but not least, it would be also very interesting (for a comparison with future heavy–ion experiments) to extend our present analysis of the radiative decays to the non–zero–temperature case. We expect that some progress will be done along this line in the near future.


format: text

text:
# A Dynamical Model for the Globular Cluster G1

## 1 Introduction

We report results from a series of $`N`$-body simulations for the globular cluster G1 in M31. G1 is one of the brightest and most massive globular clusters in the local group. Its total luminosity ($`M_V=10.94`$ mag) and central velocity dispersion ($`\sigma _0=25.1\pm 1.7`$ km/sec) are larger than those of any galactic globular cluster (Meylan et al., 2001; Djorgovski et al., 1997).

Meylan et al. (2001) used virial mass estimates and mass estimates from King-Michie models for G1. They obtained total masses in the range $`7.3\times 10^6M_{}`$ to $`15\times 10^6M_{}`$ and (for their model 4) a core radius, half-mass radius, and tidal radius of 0.53, 13.2, and 187 pc, respectively. The estimated half-mass relaxation time was 50 Gyr, much longer than the Hubble time.

Gebhardt et al. (2002) have reported evidence for an intermediate-mass black hole of $`2.0_{0.8}^{+1.4}\times 10^4M_{}`$ in the center of G1. Based on velocity profiles obtained with the Hubble space telescope’s STIS instrument, they constructed orbit-based axisymmetric models. Varying $`M/L`$ and the mass of the central black hole, they found a best fit for $`M/L=2.5`$ and $`M_{BH}=2\times 10^4M_{}`$. A model without a central black hole was rejected at a $`2\sigma `$ level.

The presence of such a black hole would be very interesting, for at least two reasons. First, it would lie neatly on the $`M_{BH}\sigma `$ relation for galaxies (Gebhardt et al., 2000; Ferrarese and Merritt, 2000). Second, G1 would then be a good example of the type of cluster postulated by Ebisuzaki et al. (2001), some of which may find their way into the center of a galaxy by dynamical friction, where their intermediate-mass black holes may then merge to provide the seeds for supermassive black holes.

However, before embracing such an exciting conclusion it is all the more important to ensure that more conventional explanations of the observational data are ruled out. To this end, we have tried to construct the best possible evolutionary model for G1 as a large globular cluster that is still in the early stages of core collapse, without harboring an intermediate-mass black hole. We have run a set of models with varying initial density profiles, half mass radii, total masses, and global $`M/L`$ until we found a model that gave the best fit to the light and velocity profiles of G1.

In §2 we describe our numerical method. In §3 we present the results of simulations starting with a single non-rotating cluster, and in §4 we show what happens when we consider G1 to be the rotating product of a merger of two smaller globular clusters. We briefly summarize in §5.

## 2 Modeling Method

In order to model the evolution of G1 using $`N`$-body simulations, we face a scaling and a fitting problem: we can only handle $`10^5`$ particles while G1 contains $`10^7`$ stars; and we do not know which values to assign to the initial cluster model parameters such as the total mass and the half-mass radius. We solve the scaling problem by scaling the dynamical parameters in such as way as to reproduce in our model simulations the correct two-body relaxation time scales inferred for G1 from observations. We solve the fitting problem by carrying out a large enough number of runs to allow us to isolate simulations that closely reproduce the observational data. Without the use of the GRAPE-6 computers (Makino et al., 2003), it would have been unpractical to run the several dozen runs needed to determine our best fits.

We used Aarseth’s $`N`$-body code NBODY4 (Aarseth, 1999). All simulated clusters contained $`N=65,536`$ stars initially, with a range of masses following Kroupa (2001)’s mass function with lower and upper mass limits of 0.1 and 30 $`\text{M}_{}`$, respectively. Our simulations did not contain primordial binaries, which is a reasonable simplification for a cluster that is still quite far from core collapse. We did not include M31’s galactic tidal field, which would have a negligible influence at the position of G1, at least 40 kpc from the center of M31. Since tidal effects are unimportant, we are left with two evolution mechanisms: stellar evolution and two-body relaxation.

Stellar evolution was modeled according to the fitting formulae of Hurley et al. (2000), using a metallicity of $`\text{[Fe/H]}=0.95`$, similar to the mean metallicity of G1 as determined by Meylan et al. (2001). We assumed a retention fraction of neutron stars of 15%.

All simulations were carried out for 13 Gyr and the final density and velocity profiles were obtained from 10 snapshots spanning a 500 Myr period centered at $`T=12`$ Gyr. For comparison of our models with the observations of G1, we assume a distance of 770 kpc to M31, so one arcsecond corresponds to 3.7 pc. Typically, about 1% of the stars escaped from the cluster during a simulation. Only bound stars were used for the comparison with observations.

We have to scale the parameters of our simulations in order to match the most important stellar evolution and stellar dynamical parameters of the actual G1 cluster. In order to match the relaxation time of G1, we have to increase the radius of our cluster by

$$r_{hS}=r_{hG1}\left(\frac{N_{G1}}{N_S}\right)^{1/3}\left(\frac{\mathrm{ln}(\gamma N_S)}{\mathrm{ln}(\gamma N_{G1})}\right)^{2/3},$$
(1)
where subscripts $`G1`$ and $`S`$ denote, respectively, the actual values for G1 and those used in our simulations. Effects which strongly depend on the number of particles in the cluster are unimportant before the cluster goes into core collapse, so our models should give a valid description of the dynamical evolution of G1 up to the present time.

In the first set of simulations, we started from isotropic King model conditions with no initial mass-segregation, and dimensionless central concentrations in the range $`4.0W_011.0`$. For each choice of initial density profile, we ran full simulations for a number of choices for the initial physical half-mass radius $`r_h(t=0)`$ and mass $`M(t=0)`$ of G1 until we could fit the surface density profile of Meylan et al. (2001) over a maximum range in radius while simultaneously obtaining an optimal fit to the observed velocity profile. For the surface density, we scaled our predicted profile by a multiplicative factor to obtain the best fit (in practice changing the M/L predicted by our assumed IMF by a factor of $`1.52`$). For the velocity profile, we used the symmetrized profile shown by Gebhardt et al. (2002) in their Fig. 1 and the ground-based value of Djorgovski et al. (1997), who measured a velocity of $`25.1\pm 1.7`$ km/sec inside an aperture of $`1\stackrel{}{\mathrm{.}}15\times 7\stackrel{}{\mathrm{.}}0`$. For each run, a best fit was determined by a $`\chi ^2`$ test against the combined data. With improved estimates for $`r_h(t=0)`$ and $`M(t=0)`$, a new initial half-mass radius could be calculated and a new simulation was performed. Simulations were performed until the half-mass radius changed by less than 5% between successive iterations. A more detailed description of our simulations and their results will be presented in a forthcoming paper (Baumgardt et al., 2003b).

## 3 Single Cluster Simulations

Figure 1 shows the data for the best fit from among the runs where we started with a single cluster in the form of a King model; here the initial central potential depth was $`W_0=7.5`$. The top panel shows the inferred projected luminosity density. We can see that the fit is very good for $`r<15`$ pc. The reason why the model density drops off sharply at large radius is because the initial King model had a tidal radius of 32 pc if we scale it to G1. Two-body relaxation begins to produce an extended halo with a surface density slope $`4`$, but a Hubble time is too short to let this effect propagate very far into the observed halo. Starting from deeper King models ($`W_0=8`$ or higher) does not solve this problem: such models predict too high a surface density around $`r=10`$ pc, while still falling short at larger radii. The implication is that G1 must have started with a density distribution more extended than any King model that can be fit to the bulk of the stars.

The bottom panel shows the velocity dispersion inferred from our $`W_0=7.5`$ model, as compared to the dispersion observed by Gebhardt et al. (2002). For larger $`W_0`$ values, our models produce velocities that are too high at the largest observed radii. Models with slightly lower concentration give a somewhat better fit, but when we require the model to reproduce the density as well, the combined requirements clearly point to $`W_0=7.5`$ as producing the best agreement, and one which falls within the observational errors everywhere except near the tidal radius artificially imposed by the initial conditions; we address this limitation in the next section.

Note that our model cluster has a mass smaller than those found from multi-mass King model fits by Meylan et al. (2001). Their extreme values stem from the implicit King-model requirement that a cluster has complete mass-segregation, which is unphysical in a massive cluster like G1 where the relaxation time is much longer than a Hubble time.

To sum up, an evolutionary model starting from a King model without initial mass segregation reproduces both the luminosity profile and the velocity dispersion profile of G1 rather well. The fits are not perfect, though, on two counts. First, the best-fit model still produces too steep a luminosity profile at larger radii. Second, since we start from a spherically symmetric non-rotating model, in principle we cannot fit the observed rotation profile or ellipticity. The question is whether we can introduce rotation while simultaneously at least preserving, and hopefully improving, the reasonable fits obtained so far. In the next section we answer this question affirmatively.

## 4 Merger Simulations

Currently favored scenarios for the formation of star clusters are the collapse of giant molecular clouds or the collision of smaller clouds (Fall and Rees, 1985; Fujimoto and Kumai, 1997). A collision scenario could easily explain the apparent rotation of G1. It might also account for the run of surface density in the halo, since simulations of the merging of two stellar systems usually give surface density profiles $`\mathrm{\Sigma }(R)R^{3.0}`$ (Sugimoto and Makino, 1989; Makino, Akiyama, & Sugimoto, 1990; Okumura et al., 1991).

Based on these theoretical hints, we have carried out a series of simulations starting with an early merger of two star clusters. For the sake of simplicity, we have restricted ourselves to the merger of two identical King model clusters on parabolic orbits with initial separation $`r_i=20`$ and pericenter separation in $`N`$-body units of $`p=1`$ (another simulation with $`p=2`$ gave similar results). The chosen pericenter distance corresponds to approximately 1.3 half-mass radii for the initial clusters. We used $`N=80,000`$ equal-mass stars in our merger simulation without including stellar evolution, a reasonable approximation given that our merger was postulated to occur during formation of the clusters. After the merger product had undergone its violent relaxation, we randomly selected 65,536 stars from among all the stars still bound to the final merger product, assigned masses drawn from a Kroupa (2001) IMF to them, and then started our dynamical evolution simulations for a duration of $`T=13`$ Gyr.

Fig. 2 shows the final density and velocity dispersion profiles for our best-fit simulations which started with a collision between two $`W_0=6.5`$ initial King models. Note that our simulations now reproduce the observed extended halo very well. The agreement between the observed and model velocity dispersions is also very good.

Figure 3 compares the rotation and ellipticity profiles of our merger model and the observations. We measured the rotation profile from two directions perpendicular to each other and the minor axis and took the mean of the two directions. The profiles were determined from the radial velocities of all bright stars located in an area between angles of $`10\mathrm{°}`$ and $`40\mathrm{°}`$ with respect to the major axis, in order to make an optimal comparison with Gebhardt et al. (2002), who performed their HST/STIS spectroscopy at an angle of $`25\mathrm{°}`$ against the major axis. The agreement between simulated and observed rotation profiles is very good.

Similarly, the ellipticity profile of our merger run is in good agreement with the ellipticity profile of G1 as determined by Meylan et al. (2001), (as can be seen in the lower panel of Fig. 3). The $`N`$-body run starts with a near constant ellipticity of about $`ϵ=0.25`$. After $`T=12`$ Gyr, the cluster core has become almost spherical due to relaxation effects, while the halo ellipticity has remained unchanged. The observations show a similar drop of $`ϵ`$ toward the core.

Our success in modeling G1 as a merger product does not necessarily imply that a merger history is the only way to explain its current state. For example, it is also possible that G1 is a heavily stripped remnant of a dwarf spheroidal. What is important is that the observed rotation can be well modeled under at least one set of reasonable assumptions, as we have shown here. The presence of rotation does not invalidate attempts at modeling under the simpler assumption of spherical symmetry; rather it invites a further fine-tuning of the already good agreement of spherical models.

Table 1 summarizes our results for the best-fitting models. It shows $`W_0`$ for each initial model, as well as the half-mass radius $`r_h`$ at $`T=12`$ Gyr, and the inferred total mass $`M`$ and $`M/L`$ required to give the best-fit velocity dispersions. The half-mass radius shown is the half-mass radius of our models after they were scaled back to G1 by eq. 1. The errors given for the $`M`$ and $`M/L`$ values are the statistical errors from the $`\chi ^2`$-fit. The global mass-to-light ratios we obtain in our best fits lie around $`M/L4`$, relatively large but still within the range of mass-to-light ratios observed for galactic globular clusters (Pryor and Meylan, 1993). Column 6 gives the probability $`P_V`$ that our velocity distribution agrees with the observations of Gebhardt et al. (2002) and Djorgovski et al. (1997), determined from a $`\chi ^2`$-test against the combined data.

The last three columns of Table 1 give the $`M/L`$ values inside the core (defined as the region containing the innermost 1% of bright stars) and the half-mass and core relaxation times calculated from the cluster parameters at $`T=12`$ Gyrs and eq. 2-62 of Spitzer (1987). Since the half-mass relaxation time is much longer than a Hubble time, G1 has not yet reached core collapse and the core of G1 is still dominated by low-mass main-sequence stars. Nevertheless, core $`M/L`$ values are smaller than global ones since mass-segregation has caused bright stars, which are more massive than average, to sink into the center. Shortly before core collapse, bright stars will be depleted from the center by the heavier neutron stars and massive white dwarfs, as in the simulations of M15 of Baumgardt et al. (2003a).

## 5 Conclusions

We have constructed evolutionary models for the massive globular cluster G1. Starting from a $`W_0=7.5`$ King model we can reproduce both the observed luminosity profile for $`r<15`$ pc and the observed velocity dispersion profile. A model starting from the merger of two $`W_0=6.5`$ King models fares even better: it can reproduce the luminosity, velocity dispersion, rotation profiles and ellipticity for the entire range of observations.

Our simulations were motivated by the recent claim of evidence for a massive central black hole. Given that our dynamical model without central black hole can fit all available observational data very well, there seems to be no need to invoke the presence of an intermediate-mass black hole. Note that we obtained an excellent fit by varying only the following basic parameters: the central potential $`W_0`$, the initial total mass $`M(0)`$ and total mass-to-light ratio $`M/L`$, and the initial half-mass radius $`r_h(0)`$. Our conclusions are therefore robust, and independent of any fine tuning in initial conditions.

The 2-sigma evidence presented by Gebhardt et al. (2002) for a massive black hole is not supported by direct observation of luminosity profile, velocity dispersion and rotation. It must have come from the data not presented in their paper (e.g. the higher order moments of the velocity profiles, together with multiparameter fits to orbit families). Without independent checks or further observational support, we consider the evidence for such a black hole to be inconclusive.

This work is the first example of the successful detailed dynamical modeling of the evolution of a globular cluster with rotation. We have shown how $`N`$-body simulations have matured as the most powerful tool to interpret detailed observational data, obviating the need for simplifying assumptions such as spherical symmetry or the use of static (e.g. multi-mass King) models.

## Acknowledgments

The authors thank Toshi Fukushige and Yoko Funato for stimulating discussions and Karl Gebhardt for sharing his velocity data on G1 with us. We are especially grateful to Sverre Aarseth for making the NBODY4 code available to us and his constant, untiring help with the code. This work is supported in part by Grant-in-Aid for Scientific Research B (13440058) of the Ministry of Education, Culture, Science and Technology, Japan, by NASA ATP grant NAG5-10775 and by the Royal Dutch Academy of Science (KNAW) and Dutch organization for Scientific Research (NWO).


text:
# Instability of speckle patterns in random media with noninstantaneous Kerr nonlinearity

## Abstract

Onset of the instability of a multiple-scattering speckle pattern in a random medium with Kerr nonlinearity is significantly affected by the noninstantaneous character of the nonlinear medium response. The fundamental time scale of the spontaneous speckle dynamics beyond the instability threshold is set by the largest of times $`T_\mathrm{D}`$ and $`\tau _{\mathrm{NL}}`$, where $`T_\mathrm{D}`$ is the time required for the multiple-scattered waves to propagate through the random sample and $`\tau _{\mathrm{NL}}`$ is the relaxation time of the nonlinearity. Inertial nature of the nonlinearity should complicate the experimental observation of the instability phenomenon.

Accepted for publication in Optics Letters © 2003 Optical Society of America, Inc.

Instabilities, self-oscillations and chaos are widely encountered in nonlinear optical systems gibbs85 ; agrawal95 ; arrechi99 . In particular, the so-called modulation instability (MI), consisting in the growth of small perturbations of the wave amplitude and phase, resulting from the combined effect of nonlinearity and diffraction, can be observed for both coherent agrawal95 ; shih02  and incoherent soljacic00 ; kip00 ; klinger01 ; chen02 ; hall02  light. In the latter case, propagation of a powerful, spatially incoherent quasimonochromatic light beam through an optically homogeneous nonlinear crystal results in ordered or disordered spatial patterns soljacic00 ; kip00 ; klinger01 ; chen02 ; hall02 . Recently, instabilities have also been predicted skip00 ; skip03  for diffuse light in random media with Kerr nonlinearity (sample size $`L`$ mean free path $`\mathrm{}`$ wavelength $`\lambda `$). The strength of the nonlinearity has to exceed a threshold for the instability to develop, similarly to the case of incoherent MI (coherent MI has no threshold), but in contrast to the latter, (a) initially coherent wave loses its spatial coherence inside the medium, due to the multiple scattering on heterogeneities of the refractive index, (b) temporal coherence of the incident light is assumed to be perfect (in the case of incoherent MI, the coherence time of the incident beam is much shorter than the response time of the medium), (c) incident light beam is completely destroyed by the scattering after a distance $`\mathrm{}`$ and light propagation is diffusive in the bulk of the sample, and (d) scattering provides a distributed feedback mechanism, absent in the case of MI.

Since the noninstantaneous nature of nonlinearity proved to be important for the development of optical instabilities gibbs85 ; agrawal95 ; arrechi99 ; shih02 ; soljacic00 ; kip00 ; klinger01 ; chen02 ; hall02 , it is natural to consider its effect on the instability of diffuse waves in random nonlinear media. Up to now, the latter phenomenon has only been studied assuming the instantaneous medium response skip00 ; skip03 . In the present letter, we show that while the absolute instability threshold does not depend on the nonlinearity response time $`\tau _{\mathrm{NL}}`$, the dynamics of the speckle pattern above the threshold is extremely sensitive to the relation between $`\tau _{\mathrm{NL}}`$ and the time $`T_\mathrm{D}=L^2/D`$ that takes the diffuse wave to propagate through the random sample ($`D=c\mathrm{}/3`$ is the diffusion constant and $`c`$ is the speed of light in the average medium). We derive explicit expressions for the characteristic time scale $`\tau `$ of spontaneous intensity fluctuations beyond the instability threshold and the maximal Lyapunov exponent $`\mathrm{\Lambda }_{\mathrm{max}}`$, characterizing the predictability of the system behavior. Finally, we discuss the experimental implications of our results.

We consider a plane monochromatic wave of frequency $`\omega _0`$ incident on a random sample of typical size $`L\mathrm{}`$. The real dielectric constant of the sample can be written as $`\epsilon (𝐫,t)=\epsilon _0+\delta \epsilon (𝐫)+\mathrm{\Delta }\epsilon _{\mathrm{NL}}(𝐫,t)`$, where $`\epsilon _0`$ and $`\delta \epsilon (𝐫)`$ are the average value and the fluctuating part of the linear dielectric constant, respectively, and $`\mathrm{\Delta }\epsilon _{\mathrm{NL}}(𝐫,t)`$ denotes the nonlinear part of $`\epsilon `$. Without loss of generality, we assume $`\epsilon _0=1`$ and $`\delta \epsilon (𝐫)\delta \epsilon (𝐫_1)=4\pi /(k_0^4\mathrm{})\delta (𝐫𝐫_1)`$. Under conditions of weak scattering $`k_0\mathrm{}1`$ (where $`k_0=\omega _0/c`$), the average intensity $`I(𝐫)`$ of the scattered wave inside the sample can be described by a diffusion equation, the short-range correlation function of intensity fluctuations $`\delta I(𝐫,t)=I(𝐫,t)I(𝐫)`$ decreases to zero at distances of order $`\lambda `$, and the long-range correlation of $`\delta I(𝐫,t)`$ can be found from the Langevin equation zyuzin87  (see Ref. rossum99,  for a review of wave diffusion in random media):

$`{\displaystyle \frac{}{t}}\delta I(𝐫,t)D^2\delta I(𝐫,t)=𝐣_{\mathrm{ext}}(𝐫,t).`$ (1)
In a nonlinear medium, the random external Langevin currents $`𝐣_{\mathrm{ext}}(𝐫,t)`$ obey a dynamic equation skip03

$`{\displaystyle \frac{}{t}}𝐣_{\mathrm{ext}}(𝐫,t)`$ $`=`$ $`{\displaystyle _V}\mathrm{d}^3𝐫^{}{\displaystyle _0^{\mathrm{}}}d\mathrm{\Delta }t𝐪(𝐫,𝐫^{},\mathrm{\Delta }t)`$ (2)
$`\times `$ $`{\displaystyle \frac{}{t}}\mathrm{\Delta }\epsilon _{\mathrm{NL}}(𝐫^{},t\mathrm{\Delta }t),`$
where the integration is over the volume $`V`$ of the random medium and $`𝐪(𝐫,𝐫^{},\mathrm{\Delta }t)`$ is a random, sample-specific response function with zero average and a correlation function given by the diagrams of Fig. 1.

Supplemented by a Debye-relaxation relation for $`\mathrm{\Delta }\epsilon _{\mathrm{NL}}(𝐫,t)`$:

$`\tau _{\mathrm{NL}}{\displaystyle \frac{}{t}}\mathrm{\Delta }\epsilon _{\mathrm{NL}}(𝐫,t)=\mathrm{\Delta }\epsilon _{\mathrm{NL}}(𝐫,t)+2n_2I(𝐫,t),`$ (3)
where $`n_2`$ is the nonlinear coefficient, Eqs. (1) and (2) can now be used to perform a linear stability analysis of the multiple-scattering speckle pattern in a random medium with noninstantaneous nonlinearity. Assuming $`\delta I(𝐫,t)=\delta I(𝐫,\alpha )\mathrm{exp}(\alpha t)`$ with $`\alpha =i\mathrm{\Omega }+\mathrm{\Lambda }`$ ($`\mathrm{\Omega }>0`$), and similarly for $`𝐣_{\mathrm{ext}}(𝐫,t)`$ and $`\mathrm{\Delta }\epsilon _{\mathrm{NL}}(𝐫,t)`$, we use Eqs. (2) and (3) to express the correlation function of Langevin currents $`j_{\mathrm{ext}}^{(i)}(𝐫,\alpha )j_{\mathrm{ext}}^{(j)}(𝐫_1,\alpha )`$ in terms of the intensity correlation function $`\delta I(𝐫,\alpha )\delta I^{}(𝐫_1,\alpha )`$ and Eq. (1) to express the latter correlation function in terms of the former. For large enough sample size $`L/\mathrm{}k_0\mathrm{}`$ and moderate frequency $`\mathrm{\Omega }T_\mathrm{D}[L/(k_0\mathrm{}^2)]^2`$ we can neglect the short-range correlation of intensity fluctuations and write the condition of consistency of the two obtained equations \[and hence the condition of consistency of Eqs. (13)\] as

$`pF(\mathrm{\Omega }T_\mathrm{D},\mathrm{\Lambda }T_\mathrm{D})H(\mathrm{\Omega }\tau _{\mathrm{NL}},\mathrm{\Lambda }\tau _{\mathrm{NL}}),`$ (4)
where a numerical factor of order unity has been omitted, $`p=\mathrm{\Delta }n_{\mathrm{NL}}^2(L/\mathrm{})^3`$ is the effective nonlinearity parameter, $`\mathrm{\Delta }n_{\mathrm{NL}}=n_2I(𝐫)`$ is the average value of the nonlinear part of refractive index, and $`I(𝐫)`$ is assumed to be approximately constant inside the sample. The function $`F`$ in Eq. (4) is the same as in the case of instantaneous nonlinearity skip03 , while $`H(x,y)=x^2+(1+y)^2`$ originates from the noninstantaneous nature of the nonlinear response. In the limit of $`\tau _{\mathrm{NL}}0`$ (instantaneous nonlinearity), $`H1`$ and Eq. (4) reduces to Eq. (7) of Ref. skip03, .

For given $`\mathrm{\Omega }`$ and $`p`$, Eq. (4) determines the corresponding value of the Lyapunov exponent $`\mathrm{\Lambda }`$. $`\mathrm{\Lambda }>0`$ signifies the instability of the speckle pattern with respect to the excitations at frequency $`\mathrm{\Omega }`$. If $`p<1`$, then $`\mathrm{\Lambda }<0`$ for all $`\mathrm{\Omega }`$. As $`p`$ exceeds $`1`$, $`\mathrm{\Lambda }`$ becomes positive inside some frequency band $`(0,\mathrm{\Omega }_{\mathrm{max}})`$ (see the inset of Fig. 2). We show the frequency-dependent threshold following from Eq. (4) in Fig. 2 for different values of $`\tau _{\mathrm{NL}}/T_\mathrm{D}`$. The threshold value of $`p`$ at $`\mathrm{\Omega }\tau _{\mathrm{NL}}^1`$ deviates from the result obtained for the instantaneous nonlinearity (dashed line in Fig. 2), while the absolute threshold is $`p=1`$ at $`\mathrm{\Omega }=0`$, independent of $`\tau _{\mathrm{NL}}/T_\mathrm{D}`$.

To study the onset of the instability, we consider the limit of $`0<p11`$ in more detail. In this limit, a series expansion of the function $`F`$ in Eq. (4) can be found: $`F(x,y)1+C_1x^2+y(C_2C_3x^2)`$, where we keep only the terms quadratic in $`x`$ and linear in $`y`$, and the numerical constants are $`C_10.031`$, $`C_20.59`$, $`C_30.045`$. At a given frequency $`\mathrm{\Omega }`$, the threshold value of $`p`$ then becomes $`p1+\left[C_1+(\tau _{\mathrm{NL}}/T_\mathrm{D})^2\right](\mathrm{\Omega }T_\mathrm{D})^2`$, while at a given $`p`$ the maximum excited frequency is $`\mathrm{\Omega }_{\mathrm{max}}=T_\mathrm{D}^1\left[C_1+(\tau _{\mathrm{NL}}/T_\mathrm{D})^2\right]^{1/2}(p1)^{1/2}`$ and the shortest typical time of spontaneous intensity fluctuations is $`\tau =1/\mathrm{\Omega }_{\mathrm{max}}`$.

As follows from the above analysis, a continuous low-frequency spectrum of frequencies $`(0,\mathrm{\Omega }_{\mathrm{max}})`$ is excited at $`p>1`$, and the Lyapunov exponent $`\mathrm{\Lambda }`$ decreases monotonically with $`\mathrm{\Omega }`$ (see the inset of Fig. 2). This allows us to hypothesize that at $`p=1`$ the speckle pattern undergoes a transition from a stationary to chaotic state. Such a behavior should be contrasted from the “route to chaos” through a sequence of bifurcations, characteristic of many nonlinear systems gibbs85 ; agrawal95 ; arrechi99 . Given initial conditions at $`t=0`$ and deterministic dynamic equations, the evolution of a chaotic system can only be predicted for $`t<1/\mathrm{\Lambda }_{\mathrm{max}}`$. We find the maximal Lyapunov exponent $`\mathrm{\Lambda }_{\mathrm{max}}=T_\mathrm{D}^1(p1)/(C_2+2\tau _{\mathrm{NL}}/T_\mathrm{D})`$. It is now easy to see that the fundamental time scale of the speckle pattern dynamics at $`p>1`$ is set by the largest of the times $`T_\mathrm{D}`$ and $`\tau _{\mathrm{NL}}`$. Indeed, in the case of fast nonlinearity ($`\tau _{\mathrm{NL}}T_\mathrm{D}`$) we have $`\mathrm{\Omega }_{\mathrm{max}}T_\mathrm{D}^1(p1)^{1/2}`$, $`\tau T_\mathrm{D}(p1)^{1/2}`$, and $`\mathrm{\Lambda }_{\mathrm{max}}T_\mathrm{D}^1(p1)`$. In the opposite limit of slow nonlinear response ($`\tau _{\mathrm{NL}}T_\mathrm{D}`$) we find $`\mathrm{\Omega }_{\mathrm{max}}\tau _{\mathrm{NL}}^1(p1)^{1/2}`$, $`\tau \tau _{\mathrm{NL}}(p1)^{1/2}`$, and $`\mathrm{\Lambda }_{\mathrm{max}}\tau _{\mathrm{NL}}^1(p1)`$, respectively.

The origin of the instability of the speckle pattern $`I(𝐫,t)`$ in a nonlinear random medium and its sensitivity to $`\tau _{\mathrm{NL}}`$ can be qualitatively understood by considering $`I`$ as a result of interference of many partial waves traveling along various diffusion paths. An infinitely small perturbation of $`I`$ at some point $`𝐫_0`$ at time $`t_0`$ diffuses throughout the medium \[according to Eq. (1)\] and, after some time $`\mathrm{\Delta }t`$, modifies the speckle pattern inside a sphere of radius $`R(D\mathrm{\Delta }t)^{1/2}`$ around $`𝐫_0`$. Due to the nonlinearity of the medium, this leads to a change of the nonlinear part of the dielectric constant $`\mathrm{\Delta }\epsilon _{\mathrm{NL}}`$ inside the same sphere after a time $`\tau _{\mathrm{NL}}`$ \[according to Eq. (3)\] and, consequently, to a modification of the phases of the partial diffuse waves that interfere to produce a new speckle pattern $`I(𝐫,t)`$. If the strength of the nonlinearity exceeds the threshold given by Eq. (4), $`I(𝐫_0,t)I(𝐫_0,t_0)`$ can be larger than the initial perturbation. The latter is hence amplified and the speckle pattern develops spontaneous dynamics. Noninstantaneous nature of the nonlinearity leads to the damping of the high-frequency components in the spectrum of the initial intensity perturbation, raising the instability threshold for $`\mathrm{\Omega }\tau _{\mathrm{NL}}^1`$, as we show in Fig. 2. This phenomenon is similar in its origin to the increase of the speed of coherent MI patterns with decreasing $`\tau _{\mathrm{NL}}`$ in a homogeneous nonlinear medium shih02 . Similarly to the case of coherent MI shih02 , the instability is arrested in the limit of $`\tau _{\mathrm{NL}}\mathrm{}`$, in contrast to the case of incoherent MI, where $`\tau _{\mathrm{NL}}\mathrm{}`$ (more precisely, $`\tau _{\mathrm{NL}}`$ than the coherence time of the incident beam) is an explicit assumption of the theoretical models soljacic00 ; kip00 ; hall02 .

To conclude, we discuss the experimental implications of our results. $`\mathrm{\Delta }n_{\mathrm{NL}}`$ up to $`10^3`$ can be realized in inorganic photorefractive crystals kip00 ; klinger01 ; chen02 . If scattering centers (optical defects) are introduced in the crystal, $`L/\mathrm{}10^2`$ will suffice to reach the instability threshold $`p1`$. Alternatively, in nematic liquid crystals both nonlinearity ($`\mathrm{\Delta }n_{\mathrm{NL}}`$ up to 0.1, see, e.g., Ref. muenster97, ) and scattering ($`L/\mathrm{}>10`$ has been reported in Ref. kao96, ) can be sufficiently strong to obtain $`p1`$. The nonlinearity is rather slow in the above cases, $`\tau _{\mathrm{NL}}`$ exceeds $`T_\mathrm{D}`$ ($`T_\mathrm{D}10÷100`$ ns for $`\mathrm{}0.1÷1`$ mm and $`L1`$ cm), and the results obtained in the present letter are, therefore, of particular relevance. Finally, although the noninstantaneous nature of nonlinearity does not affect the absolute instability threshold, it should complicate the experimental observation of the instability phenomenon. Indeed, in a real experiment only sufficiently fast fluctuations can be detected. Inasmuch as $`\mathrm{\Omega }_{\mathrm{max}}`$ decreases with $`\tau _{\mathrm{NL}}`$, larger $`p`$ (and thus stronger nonlinearities) will be required for the instability to be detectable in a medium with noninstantaneous nonlinearity.