text:
# Suppression of inhomogeneous broadening in rf spectroscopy of optically trapped atoms

## Abstract

We present a novel method for reducing the inhomogeneous frequency broadening in the hyperfine splitting of the ground state of optically trapped atoms. This reduction is achieved by the addition of a weak light field, spatially mode-matched with the trapping field and whose frequency is tuned in-between the two hyperfine levels. We experimentally demonstrate the new scheme with 85Rb atoms, and report a 50-fold narrowing of the rf spectrum.

The time available for measuring an atomic transition in a perturbation-free environment has been substantially increased by the achievement of ultracold atomic samples. Using cold atoms in atomic fountains, measurements of the 9.2 GHz ”clock” transition in Cesium were performed with measurement times as long as 450 msKasevich89 . A possible way to increase the measurement time beyond the practical limit imposed by the height of a fountain is to use optically trapped atomsChu86 ; Grimm00 . Far-off-resonance optical traps (FORTs)Miller93 , which are based on the conservative dipole force created by cycles of absorption and stimulated emission, are possible candidates for these measurements. The great disadvantage of a trap, is that the trapping potential acts also as a perturbation for the atomic level, and in particular an optical trap introduces a relative ac Stark shift of the hyperfine levelsDavidson95 , which results in a systematic shift in the clock frequency measurement. When ensemble-averaged, the spatial dependence of the potential also results in a inhomogeneous broadening of the transition, and hence in a loss of atomic coherence at a much faster rate than the spontaneous photon scattering rateDavidson95 . These effects can be reduced by increasing the trap detuning Davidson95 ; Adams95 ; Miller93 ; Takekoshi96  and by using blue-detuned optical traps, in which the atoms are confined mainly in the darkDavidson95 ; Friedman02 . However, the residual frequency shifts are still the main limiting factor for precision spectroscopy in optical traps.

ac Stark shifts from the trapping beams are also detrimental to achieving high phase space densities in optical traps, since they shift the atom’s resonance frequency away from the cooling beams frequency Gordon80 ; Kuppens00 . Recently, it was shown that in some cases a FORT laser frequency can be chosen to couple the states involved in cooling to some other states, in order to suppress the frequency shift of the cooling transition and allow simultaneous trapping and Doppler coolingKatori99 . In this way, a phase-space density exceeding 0.1 was achieved.

In this letter, we demonstrate a method for reducing the inhomogeneous broadening in the spectroscopic measurement of the hyperfine splitting of the ground state of optically trapped atoms. This reduction is achieved by the addition of a very weak light field (the so-called compensating beam), spatially mode-matched with the trapping field and whose frequency is tuned between the two hyperfine levels. An experimental realization of the new scheme is presented with 85Rb atoms, and a reduction by a factor 50 is reported. This method can be applied to red- or blue-detuned FORTs, and hence can be used as an additional way to further increase the long atomic coherence times of the latter.

A ground state $`|g_i`$ of an atom exposed to a light field with intensity $`I`$ and frequency $`\omega `$, is ac Stark shifted by an amount given by:

$$\mathrm{\Delta }E_i=\frac{3\pi c^2\gamma }{2\omega _0^3}I\times \frac{c_{ij}^2}{\delta _{ij}}$$
(1)
where $`\gamma `$ is the natural width of the transitionGrimm00 . The summation takes into account the contributions of the different coupled excited levels $`|e_j`$, each with its respective transition coefficient $`c_{ij}`$, and detuning $`\delta _{ij}=\omega \omega _{ij}`$.

Specifically for the D2 line in 85Rb, the fine structure of the excited state ($`\mathrm{\Delta }_F^{}7`$ THz), and the hyperfine structure of the ground state ($`\mathrm{\Delta }_{HF}3`$ GHz) and excited state ($`\mathrm{\Delta }_{HF}^{}`$, tens of MHz) obey $`\mathrm{\Delta }_F^{}\mathrm{\Delta }_{HF}\mathrm{\Delta }_{HF}^{}`$. For linearly polarized light in the vicinity of this line, and as long as the detuning of the light is large as compared to the excited state hyperfine splitting $`\mathrm{\Delta }_{HF}^{}`$, a general result can be derived from Eq. 1 for the shift of a ground state with total angular momentum $`F`$, exposed to a light field $`I\left(𝐫\right)`$:

$$\mathrm{\Delta }E\left(𝐫\right)=\frac{\pi c^2\gamma }{\omega _0^3}\frac{I\left(𝐫\right)}{\delta _F}$$
(2)
where $`\delta _F=\omega \omega _F`$ is the detuning of the laser from the $`|5S_{1/2},F|5P_{3/2}`$ transition. Note, that Eq. 2 is a reasonable approximation even for a detuning $`\delta _F`$ comparable with the optical frequency $`\omega _0`$ Adams95 . We are interested in rf spectroscopy, where the energy difference between the two ground state hyperfine levels, $`|F=2`$ and $`|F=3`$, is measured. In the presence of the light this energy difference is modified by:

$$\stackrel{~}{\mathrm{\Delta }}_{HF}\left(𝐫\right)\mathrm{\Delta }_{HF}=\frac{\pi c^2\gamma }{\omega _0^3}\frac{\mathrm{\Delta }_{HF}}{\delta ^2}\left[\frac{1}{1\left(\frac{\mathrm{\Delta }_{HF}}{2\delta }\right)^2}\right]I\left(𝐫\right)$$
(3)
where $`\stackrel{~}{\mathrm{\Delta }}_{HF}\left(𝐫\right)`$ is the spatially dependent hyperfine splitting in the presence of the light and $`\delta (\delta _{F=2}+\delta _{F=3})/2`$ measures the laser detuning from the center of the ground state hyperfine splitting. Equation 3 indicates that for $`\left|\delta \right|>\frac{\mathrm{\Delta }_{HF}}{2}`$, i.e. for $`\delta _{F=2}`$ and $`\delta _{F=3}`$ both positive (or both negative), the ground state energy splitting is always reduced by the presence of a light field (See Fig. 1). When the detuning is ”between” the hyperfine levels, $`\left|\delta \right|<\frac{\mathrm{\Delta }_{HF}}{2}`$, the energy splitting is enlarged.

For a FORT with detuning $`\delta \frac{\mathrm{\Delta }_{HF}}{2}`$ Eqs. 2,3 yield:

$$\stackrel{~}{\mathrm{\Delta }}_{HF}\left(𝐫\right)\mathrm{\Delta }_{HF}\left(\frac{\mathrm{\Delta }_{HF}}{\delta }\right)U\left(𝐫\right)$$
(4)
where $`U\left(𝐫\right)`$ is the spatially dependent dipole potential that forms the trap, in frequency units. The fact that the relative ac Stark shift is $`\mathrm{\Delta }_{HF}\mathrm{}\delta `$ smaller than the dipole potential, is the main motivation for using FORTs for precision spectroscopy. For example in Ref. Davidson95  the relative ac Stark shifts were only $`\mathrm{\Delta }_{HF}\mathrm{}\delta 210^4`$ times the dipole potential. Note that, however small, this relative ac Stark shift is still much larger than the spontaneous photon scattering rate Davidson95 .

For a trapped atomic ensemble, $`I\left(𝐫\right)`$ will cause a shift in the ensemble averaged ground state hyperfine energy splitting, $`\stackrel{~}{\mathrm{\Delta }}_{HF}\mathrm{\Delta }_{HF}`$. In addition, the spatial dependence of $`I\left(𝐫\right)`$ will result in an inhomogeneous broadening, $`\sqrt{\stackrel{~}{\mathrm{\Delta }}_{HF}^2\stackrel{~}{\mathrm{\Delta }}_{HF}^2}`$, which depends also on the atoms’ temperatureNodynamics . For example, for a thermal ensemble of atoms with $`\frac{3}{2}k_bT`$ kinetic energy in an harmonic trap, we have $`\stackrel{~}{\mathrm{\Delta }}_{HF}\mathrm{\Delta }_{HF}=\left(\frac{\mathrm{\Delta }_{HF}}{\delta }\right)\left[U_0+\frac{3}{2}k_bT\right]`$ and $`\sqrt{\stackrel{~}{\mathrm{\Delta }}_{HF}^2\stackrel{~}{\mathrm{\Delta }}_{HF}^2}=\left(\frac{\mathrm{\Delta }_{HF}}{\delta }\right)\sqrt{\frac{3}{2}}k_bT`$ , where $`U_0`$ is the dipole potential at the trap’s bottom Kaplan02 .

In order to cancel these shifts, we introduce an additional laser beam, with intensity $`I^{}\left(𝐫\right)`$ and frequency between the resonant frequencies of the two ground state hyperfine levels, say in the middle, i.e. $`\mathrm{\Delta }_{HF}/2`$ and $`+\mathrm{\Delta }_{HF}/2`$ from the lower and higher hyperfine level respectivelyNomiddle . For a FORT, the total shift is obtained by adding the shifts from the trap and the compensating beam,

$$\stackrel{~}{\mathrm{\Delta }}_{HF}\left(𝐫\right)\mathrm{\Delta }_{HF}=\frac{\pi c^2\gamma \mathrm{\Delta }_{HF}}{\omega _0^3}\left[\frac{I\left(𝐫\right)}{\delta ^2\left(\frac{\mathrm{\Delta }_{HF}}{2}\right)^2}\frac{I^{}\left(𝐫\right)}{\left(\frac{\mathrm{\Delta }_{HF}}{2}\right)^2}\right]$$
(5)
If the compensating beam is spatially mode-matched with the trap beam, i.e. $`I^{}\left(𝐫\right)=\eta \times I\left(𝐫\right)`$, then a complete cancellation of the inhomogeneous broadening will occur for

$$\eta =\frac{(\frac{\mathrm{\Delta }_{HF}}{2})^2}{\delta ^2\left(\frac{\mathrm{\Delta }_{HF}}{2}\right)^2}(\frac{\mathrm{\Delta }_{HF}}{2\delta })^2$$
(6)
Fig. 1 shows the shift of the hyperfine levels (a) and the hyperfine energy difference (b) caused by the trapping beam (dashed line) and compensating beam (dotted line). In the presence of both beams, the levels are shifted by the same amount (full line). As a specific example, with a FORT detuned by 5 nm we have $`\eta 3.6\times 10^7`$. Hence, with a typical FORT power of 50 mW, the required compensating beam power is 20 nWnoted1 .

Note, that the dipole potential created by the compensating beam is $`U^{}\left(𝐫\right)=\pm \frac{1}{2}\frac{\mathrm{\Delta }_{HF}}{\delta }U\left(𝐫\right)`$ for atoms in the upper an lower hyperfine level, respectively, and hence is negligible when compared with the potential of the FORT. Moreover, the photon scattering rate from the compensating beam $`\gamma _s^{}`$ is given by $`\mathrm{}\gamma _s^{}2\gamma U^{}/\mathrm{\Delta }_{HF}`$, which can be also written as $`\mathrm{}\gamma _s^{}\frac{\gamma }{\delta }U\left(𝐫\right)\mathrm{}\gamma _s`$. Hence, the scattering rate from the nearly resonant compensating beam, is of similar magnitude to that of the far-of-resonance trapping beam, $`\gamma _s`$.

We implement the proposed scheme with a FORT (a red-detuned gaussian beam Chu86 ), created by focusing a 50 mW laser, detuned 5 nm below resonance, to a waist of $`W_0=50\mu m`$, resulting in a potential depth of $`U_0200E_{rec}`$ ($`E_{rec}`$ is the recoil energy) and oscillation frequencies of $`\nu _r=2.3\times 10^3`$ Hz and $`\nu _z=13`$ Hz in the radial and axial directions, respectively. An additional laser, with frequency locked close to the middle of the ground state hyperfine splitting, is combined with the FORT laser. To achieve optimal spatial mode-match, both lasers are coupled into a polarization-conserving single-mode optical fiber, and the fiber’s output is passed through a polarizer and focused into the vacuum chamber. Two servo loops are used to control and stabilize the power of the lasers: The first one ensures a $`1\%`$ stability of the trap laser. More importantly, for complete compensation of the relative ac Sark shifts, a second servo loop ensures a $`0.1\%`$ stability of the power ratio $`\eta `$ throughout the entire duration of the experimentGratings .

The loading procedure is similar to that described in Ozeri99 . Briefly, the FORT is loaded by overlapping it with a compressed 85Rb magneto-optical trap (MOT). The MOT beams are shut off after 650 ms of loading, 50 ms of compression, and 5 ms of polarization gradient cooling, leaving $`10^5`$ confined atoms with a temperature of $`10\mu K`$.

We perform a Rabi spectroscopy measurement on the trapped atoms by driving the ground state $`|F=2,m_F=0`$ $`|F=3,m_F=0`$ transition, which is insensitive to magnetic fields, to the first order. A bias magnetic field of $``$ 80mG is applied parallel to the FORT’s polarization axis and to the rf magnetic field direction, in order to Zeeman shift the magnetic sensitive $`m_F0`$ levels out of resonance with the rf pulse. A typical sequence is as follows: First, the atoms are prepared in the $`|F=2`$ ground state by turning on the MOT beams, without a repump beamOzeri99 , for 1 ms. Then, an rf pulse in applied at a variable frequency, using an Anritsu 69317B Signal Generator locked to a high stability oscillator. The intensity and duration of the pulse are adjusted to maximize the $`|F=3`$ population when on-resonance ($`\pi `$ pulse condition). Following the rf pulse, $`N_3`$ (the population in the $`|F=3`$ level) is measured by detecting the fluorescence during a short pulse of a laser beam resonant with the cycling transition $`|5S_{1/2},F=3|5P_{3/2},F=4`$. The population of the $`|F=2`$ level ($`N_2`$) is then measured by turning on the repumping beam (which is resonant with $`|5S_{1/2},F=2|5P_{3/2},F=3`$) and applying an additional detection pulse. The normalized signal $`N_3/(N_2+N_3)`$ is insensitive to shot-to-shot fluctuations in atom number as well as slow drifting fluctuations of the detection laser frequency and intensityKhaykovich99 .

Figure 2 shows results for the Rabi spectrum with a 3 ms long $`\pi `$ pulse. A constant background resulting from spontaneous Raman scattering Cline94  is substracted. The spectrum of free-falling atoms shows no inhomogeneous broadening and a rms width, $`\sigma `$, which is Fourier limited to $`110Hz`$. A shift in the peak frequency ($`756Hz`$), and a broadening of the line (to $`\sigma =320Hz`$) are seen in the spectrum of trapped atoms, in fair agreement with the calculated trap depth and atomic temperature. This inhomogeneous broadening is not significantly affected by the duration of the pulse Nodynamics . The addition of the weak compensating beam, nearly cancels the broadening of the spectrum, as well as its shift from the free-atom line center.

Figure 3 shows the measured rms width and shift of the trapped atoms as a function of compensating beam power, again for a 3 ms rf pulse. The spectrum width is minimized to the Fourier broadening limit at a compensating beam power which corresponds also to a minimal shift from the free-atoms line center.

Figure 4 shows the measured rf spectrum for a 25 ms long pulse. A measurement of free atoms with this pulse length is not possible in our setup since the atoms fall due to gravity, and leave the interaction region. A Fourier limited $`\sigma =13Hz`$ is measured, representing a 25-fold reduction in the line broadening, as compared to the line broadening of trapped atoms. We performed a similar measurement with a 50 ms pulse, and observed a nearly Fourier limited width (a 50-fold narrowing), at the expense of a much larger spontaneous photon scattering and hence a smaller signalBackground5 . For even larger measurement times spontaneous photon scattering prevents further narrowing of the line.

We measure the spin relaxation rateCline94 ; Ozeri99  to be $`3\times 10^3s^1`$ for atoms trapped in the FORT. The addition of the compensating beam induces an increase of only $`20\%`$.

In summary, we perform an rf spectroscopy measurement of the hyperfine splitting of the ground state of optically trapped atoms. We demonstrate a novel scheme for eliminating the trap-induced inhomogeneous broadening of the transition, by adding a weak ”compensating” laser, spatially mode-matched with the trapping laser and with a proper detuning and intensity. Despite being tuned close to resonance, this laser induces a negligible change in the dipole potential, and does not considerably increase the spontaneous scattering rate. With the suppression of inhomogeneous broadening, the atomic coherence time is now limited by the much smaller spontaneous scattering time.

Whereas in a conventional optical trap the ac Stark shift of the line center strongly depends on the temperature of the atoms, which may drift considerably, in the compensated trap the suppression of the line shift is equally effective for all temperatures. Hence, it provides a mean of achieving a higher stability of the line center than that achieved by simply stabilizing the trapping laser detuning and intensity. For relative spectroscopic measurements, such as the proposed measurement of the electron’s permanent electric dipole moment (EDM)Chin01 , only stability (and not absolute accuracy) of the line center is of importance. For example, for a 10 $`\mu `$K deep YAG-laser trap, and a compensating beam with a 15 KHz (time averaged) frequency stability, locking the relative intensity between both beams to a $`1:10^5`$ stability, will result in $`10^{14}`$ stability of the rf line center.

Finally, a weak compensating beam, spatially mode-matched with the trapping beam and properly tuned near resonance between the upper level of a laser cooling transition and another excited level, can suppress spatially dependent frequency shifts of the cooling transition, and allow simultaneous trapping and cooling with more flexibility than the single frequency method of Ref. Katori99 .

This work was supported in part by the Israel Science Foundation, the Minerva Foundation, and the United States-Israel Binational Science Foundation. MFA acknowledges help from the Nachemsohn Dansk-Israelsk Studienfond.


format: text

text:
# 𝑈_𝐴(1) Symmetry Breaking and the Scalar Sector of QCD

## I Introduction

The light quark sector of QCD acquires a $`U(N_f)_LU(N_f)_R`$ symmetry and scale invariance in the limit of $`N_f`$ massless quarks. This limit is a reasonable starting point for constructing effective theories. The $`u`$ and $`d`$ quark masses are small enough compared to the naive confinement scale ($`\mathrm{\Lambda }_{QCD}`$) to make the procedure reliable. Incorporating the $`s`$ quark  makes sense once one realizes that the relevant scale is $`\mathrm{\Lambda }_\chi 4\pi f_\pi `$ . The conventional approach to low energy QCD is to assume that the octet of lowest lying pseudoscalar mesons $`\{\pi ,K,\eta \}`$ are approximate Goldstone Bosons (GB) associated with the spontaneous breakdown of the $`SU(3)_LSU(3)_R`$ symmetry. This belief is supported by the smallness of $`\pi `$, $`K`$ and $`\eta `$ meson masses as compared to the typical hadronic scale of $`1GeV`$. The approach leads to a so far successful description of pseudoscalar meson physics where the relevant QCD Green functions are systematically expanded in powers of $`/\mathrm{\Lambda }_\chi `$, $`m_q/\mathrm{\Lambda }_\chi `$, where $`\mathrm{\Lambda }_\chi `$ stands for the chiral symmetry breaking scale, in the so-called Chiral Perturbation Theory ($`\chi PT`$) expansion . The remaining pseudoscalar, the $`\eta ^{}`$ meson, has a large mass which has been related to the breakdown of the $`U_A(1)`$ symmetry .

The symmetries of massless QCD are reduced to smaller symmetries in many ways. Both scale invariance and $`U_A(1)`$ symmetry are broken at the quantum level . The residual $`SU(3)_LSU(3)_RU(1)_V`$ symmetry is spontaneously reduced down to $`SU(3)_V`$ (the vector $`U(1)`$ symmetry being trivial for mesons). Finally, if one introduces mass terms, $`SU(3)_V`$ reduces down to isospin, or completely breaks down, depending on the choice for the quark masses.

The most important contribution to the breakdown of the $`U(1)_A`$ symmetry comes from Euclidean classical field configurations with non-trivial topology (instantons). These effects are modulated by the $`e^{8\pi ^2/g_s^2}`$ factor, where $`g_s`$ denotes the strong coupling constant, thus being small at high energies but becoming important at low energies .

Physics beyond the low energy region (where $`\chi PT`$ is strictly valid) and up to the $`\varphi (1020)`$ mass requires to introduce all the relevant degrees of freedom, namely, the $`\eta ^{}`$ meson, the well established $`a_0(980)`$ and $`f_0(980)`$ scalar mesons and the lowest lying vector meson nonet. The formulation of effective theories including the $`\eta ^{}`$ and its relation to the $`U(1)_A`$ anomaly was firstly done in . In some of these pioneering works  the starting point is a Lagrangian which exhibits chiral symmetry realized in a linear way. This amounts to consider scalar fields as the chiral partners of pseudoscalars. In particular in  the relation to the instanton induced quark interaction discovered by ’t Hooft is established by expanding one of the terms in the Lagrangian. Scalar fields are then integrated out and authors focus on the properties of the vacuum and on the pseudoscalar phenomenology.

Over the past few years the problem of the scalar mesons has been under intense debate. The identification of $`\overline{q}q`$ scalar mesons is still awaiting for a final answer. This problem is related to the unique QCD prediction for the existence of bound states of gluons (glueballs) and the possibility for the existence of hybrid states. The well established low mass scalar mesons are: the $`f_0(980)`$, $`f_0(1370)`$ and the recently resurrected $`f_0(4001200)`$ (or $`\sigma `$) in the $`I=0`$ sector; the $`a_0(980)`$ and $`a_0(1450)`$ on the isovector side, and the isospinor $`K_0^{}(1430)`$ . In addition to these states there are claims for the existence of signals of other scalar mesons below 1 GeV . In particular, the approach summarized in  which is based on a non-linear chiral Lagrangian including scalar degrees of freedom concludes the existence of a light isoscalar scalar $`\sigma `$ and an isovector scalar $`\kappa `$ with a mass around 850 MeV. This is consistent with results in  though these approaches differ in the mixing of scalars and its interpretation. The differences come from the fact that $`U_A(1)`$ anomaly effects are not taken into account in , while they are considered as fundamental in . As we shall see below the picture arising in the later case is completely consistent and most of the controversial properties of scalar mesons find a natural and simple explanation in this framework.

Many other schemes have been put forth trying to understand the lowest lying scalar meson properties. A $`\overline{q}q`$ structure for the $`a_0(980)`$ and $`f_0(980)`$ was put on doubt since the corresponding quark models were not able to explain the tiny coupling of these mesons to two photons  and the quark content suggested by the mass spectrum was incompatible with the main decay channels . More complicated models such as four quark states  and molecules  were advocated to explain these points. This is a very active field and no definitive conclusion has been reached as to which states are to be considered as $`\overline{q}q`$, multi-quark, molecule, gluonia or hybrid states .

In this work we intend to shed some light on this problem. To this end let us first exhibit some of the properties of scalar (and pseudoscalar) mesons which seem to be in conflict with a $`\overline{q}q`$ interpretation.

On the base of a naive constituent quark model we expect the $`{}_{}{}^{1}S_{0}^{}`$ $`\overline{q}q`$ states ($`q=u,d`$) to have a mass $`2m_q630`$MeV where $`m_q`$ stands for $`u`$ or $`d`$ constituent quark mass, which we consider as equal and evaluate to one third of the proton mass. In a similar way the $`\overline{q}s`$ ($`\overline{s}s`$) states should have a mass $`870MeV`$ ( $`1010`$ MeV) if we evaluate the constituent strange quark mass to $`m_sm_\mathrm{\Omega }/3555`$MeV. In this picture, a $`{}_{}{}^{3}S_{1}^{}`$ $`\overline{q}q`$ state would be slightly heavier than the $`{}_{}{}^{1}S_{0}^{}`$ state and a $`{}_{}{}^{3}P_{0}^{}`$ $`\overline{q}q`$ state should have an even larger mass.

The lowest lying vector mesons fulfill these naive expectations, it is in this sense a “well-behaved” sector. The pseudoscalar and scalar sectors, however, are not. The pions and kaons are lighter than expected from naive quark model considerations. The smallness of the pion and kaon masses is qualitatively understood on the base of spontaneous breaking of $`\chi `$-ral symmetry. However, even with this explanation in mind, there exist at least one problem in this sector which is usually overlooked and which has to do with the quark content of pseudoscalars. Consider the isoscalar-pseudoscalar states with well defined quark content (flavor states): $`\eta _{ns}(\overline{u}u+\overline{d}d)/\sqrt{2}`$ and $`\eta _s\overline{s}s`$. Since they are not the physical states, they must be mixed somehow to yield the physical $`\eta `$ and $`\eta ^{}`$ mesons. Independently of the mechanism which mixes the flavor isoscalar-pseudoscalar states, a straightforward analysis yields the relations

$$m_{\eta _{ns}}^2=m_\eta ^2\mathrm{cos}^2\varphi _P+m_\eta ^{}^2\mathrm{sin}^2\varphi _P,m_{\eta _s}^2=m_\eta ^2\mathrm{sin}^2\varphi _P+m_\eta ^{}^2\mathrm{cos}^2\varphi _P,$$
(1)
where the information on the mixing mechanism is hidden in the mixing angle in the flavor basis $`\varphi _P`$. This angle has been estimated as $`\varphi _P39.5^{}`$ , corresponding to an angle $`\theta _P15.2^{}`$ in the naive single-angle description of pseudoscalar mesons mixing in the usual singlet-octet basis. Introducing this information in Eq.(1) we obtain $`m_{\eta _{ns}}741MeV`$ and $`m_{\eta _s}816MeV`$. The obvious question here is: Why is the $`\eta _{ns}`$ ( a $`\overline{q}q`$ state differing from the pion only in isospin quantum numbers) much heavier than the pion, even more massive than the $`\overline{q}s`$ ($`K`$) and as heavy as the (purely strange) $`\eta _s`$ state ?

Concerning the scalar sector things are much more involved since even an undoubtful identification of the corresponding nonet is still missing. Assume that the $`a_0(980)`$ and $`f_0(980)`$ mesons are the $`I=1`$ and $`I=0`$ members of the $`\overline{q}q`$ scalar nonet. Their nearby degeneracy suggest that they are the scalar mesons analogous to the $`\rho (770)`$ and $`\omega (780)`$ vector mesons which are also almost degenerate in mass. On this identification grounds the $`f_0(980)`$ should be predominantly non-strange but then: How can one explain its strong coupling to $`\overline{K}K`$?. This could be qualitatively understood if the $`f_0(980)`$ had a strong strange component but then one does not understand why this meson is almost degenerate in mass with the non-strange $`a_0(980)`$. On the other hand the $`a_0(1450)`$ and the $`f_0(1350)`$ seem to be heavy enough to be out of the scope of the naive quark model considerations. The same can be said for the vector $`K_0^{}(1430)`$. Another possibility is the existence of light scalar mesons: $`\sigma (400600)`$ and $`\kappa (900)`$. It seems that we have now a general consensus on the existence of a broad scalar structure in the low energy region  though it is not clear yet what its nature is and which mechanism makes this meson so light. The $`\kappa (900)`$ is a more controversial object  and its existence has not been firmly established.

As mentioned above, the vector mesons spectrum is in qualitative agreement with naive quark model expectations. Thus, whatever the explanation for the scalar and pseudoscalar mesons spectrum may be, it should have no effect on the vector sector.

In this letter we suggest that all the points raised above can be understood by considering the instanton induced quark interaction which is also responsible for the $`U_A(1)`$ symmetry breaking. The appealing property of this interaction is that it does not affect the vector sector since vector mesons consist of quark-anti-quark pairs that are either both left-handed or both right-handed and do not match the quantum numbers of the instanton induced interaction . This, and the fact that it is able to explain why the $`\eta ^{}`$ is not a Goldstone Boson make this interaction a natural candidate for explaining the unusual properties of pseudoscalar mesons and their chiral partners, the scalar mesons. Thus, we require a model considering both nonets in a chiral symmetric way and which takes into account t’Hooft interaction.

## II The Model

A model with the characteristics mentioned at the end of the previous section has been reconsidered recently . The model, incorporates pseudoscalar and scalar degrees of freedom, exhibits chiral symmetry and incorporates $`U_A(1)`$ breaking in a phenomenological way. A $`U(2)U(2)`$ version of this model was formulated in connection with the explanation of the $`U_A(1)`$ anomaly . Amazingly, similar models were put forward even before the birth of QCD . The Lagrangian is given by

$$=_{sym}+_{SB}$$
(2)
where $`_{sym}`$ denotes the $`U(3)\times U(3)`$ symmetric Lagrangian:

$$_{sym}=<\frac{1}{2}\left(_\mu M\right)\left(^\mu M^{}\right)>\frac{\mu ^2}{2}X(\sigma ,P)\frac{\lambda }{4}Y(\sigma ,P)\frac{\lambda ^{}}{4}X^2(\sigma ,P)$$
(3)
with $`M=\sigma +iP`$, and $`X,Y`$ stand for the $`U(3)\times U(3)`$ chirally symmetric terms:

$$X(\sigma ,P)=<MM^{}>,Y(\sigma ,P)=<(MM^{})^2>.$$
(4)
We closely follow the conventions in Ref.. In particular, we use $`F\frac{1}{\sqrt{2}}\lambda _if_i`$ with $`F=\sigma ,P`$; $`f_i=\sigma _i,p_i`$; $`i=\mathrm{0..8}`$ and $`\lambda _i`$ denote Gell-Mann matrices. For further details the interested reader is referred to Ref.. Chiral and $`U_A(1)`$ symmetries are explicitly broken by

$$_{SB}=<c\sigma >\beta Z(\sigma ,P)$$
(5)
where $`c\frac{1}{\sqrt{2}}\lambda _ic_i`$, with $`c_i`$ constant and

$$Z(\sigma ,P)=\text{det}\left(M\right)+\text{det}\left(M^{}\right)$$
(6)
The most general form of $`c`$ which preserves isospin and gives PCAC is such that the only non-vanishing coefficients are $`c_0`$ and $`c_8`$. The former gives, by hand, the pseudo-scalar nonet a common mass, while the later breaks the $`SU(3)`$ symmetry down to isospin. These parameters can be related to quark masses in QCD.

The $`Z`$ term in Eq.(5) corresponds to the instanton induced quark interaction in the case when the instanton angle $`\theta _{inst}`$ is aligned to zero. This interaction has the form of a determinant in flavor space and breaks $`U(3)_LU(3)_R`$ into $`SU(3)_LSU(3)_RU(1)_V`$ .

The linear $`\sigma `$ term in Eq. (5) induces $`\sigma `$-vacuum transitions which give to $`\sigma `$ fields a non-zero vacuum expectation value (hereafter denoted by $`\{\}`$). Linear terms can be eliminated from the theory by performing a shift to a new scalar field $`S=\sigma V`$ such that $`\{S\}=0`$, where $`V\{\sigma \}`$. This shift generates new three-meson interactions and mass terms.

For the sake of simplicity let us write $`V`$=Diag$`(a,a,b)`$ where $`a,b`$ are related to $`\{\sigma \}`$ through

$$a=\frac{1}{\sqrt{3}}\{\sigma _0\}+\frac{1}{\sqrt{6}}\{\sigma _8\},b=\frac{1}{\sqrt{3}}\{\sigma _0\}\frac{2}{\sqrt{6}}\{\sigma _8\}$$
(7)
Meson masses generated by this procedure are:

i) Non-mixed sectors.

$`m_\pi ^2=\xi +2\beta b+\lambda a^2,`$ $`m_K^2=\xi +2\beta a+\lambda (a^2ab+b^2),`$ (8)
$`m_a^2=\xi 2\beta b+3\lambda a^2,`$ $`m_\kappa ^2=\xi 2\beta a+\lambda (a^2+ab+b^2),`$ (9)
where $`a`$ and $`\kappa `$ denote the scalar mesons analogous to $`\pi `$ and $`K`$ respectively and we used the shorthand notation $`\xi =\mu ^2+\lambda ^{}(2a^2+b^2)`$.

ii) Mixed sectors.

$`_2^P`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(m_{0P}^2P_0^2+m_{8P}^2P_8^2+2m_{08P}^2P_0P_8\right),`$ (10)
$`_2^S`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(m_{0S}^2S_0^2+m_{8S}^2S_8^2+2m_{08S}^2S_0S_8\right)`$ (11)
where

$`m_{8P}^2`$ $`=`$ $`\xi +{\displaystyle \frac{1}{3}}[\lambda \left(a^2+2b^2\right)+2\beta \left(4ab\right)],`$ (12)
$`m_{0P}^2`$ $`=`$ $`\xi +{\displaystyle \frac{1}{3}}[\lambda \left(2a^2+b^2\right)4\beta \left(2a+b\right)],`$ (13)
$`m_{08P}^2`$ $`=`$ $`{\displaystyle \frac{\sqrt{2}}{3}}\left(ab\right)[\lambda \left(a+b\right)+2\beta ],`$ (14)
$`m_{8S}^2`$ $`=`$ $`\xi +{\displaystyle \frac{1}{3}}[2\beta \left(4ab\right)+3\lambda \left(a^2+2b^2\right)+4\lambda ^{}\left(ab\right)^2],`$ (15)
$`m_{0S}^2`$ $`=`$ $`\xi +{\displaystyle \frac{1}{3}}[4\beta \left(2a+b\right)+3\lambda \left(2a^2+b^2\right)+2\lambda ^{}\left(2a+b\right)^2],`$ (16)
$`m_{08S}^2`$ $`=`$ $`{\displaystyle \frac{\sqrt{2}}{3}}\left(ab\right)[2\beta +3\lambda \left(a+b\right)+2\lambda ^{}\left(2a+b\right)].`$ (17)
As we are interested in the quark content of fields it is convenient to analyze the mixed sector in the flavor ($`\{|s>,|ns>\}`$) basis  defined by:

$$\eta _{ns}=\sqrt{\frac{1}{3}}P_8+\sqrt{\frac{2}{3}}P_0,\eta _s=\sqrt{\frac{2}{3}}P_8+\sqrt{\frac{1}{3}}P_0$$
(19)
with analogous relations for the scalar ($`\sigma _{ns},\sigma _s`$) mixed fields. In this representation, the mass terms in the Lagrangian read

$`_2^P={\displaystyle \frac{1}{2}}\left(m_{\eta _{ns}}^2\eta _{ns}^2+m_{\eta _s}^2\eta _s^2+2m_{\eta _{sns}}^2\eta _{ns}\eta _s\right)`$ (20)
$`_2^S={\displaystyle \frac{1}{2}}\left(m_{\sigma _{ns}}^2\sigma _{ns}^2+m_{\sigma _s}^2\sigma _s^2+2m_{\sigma _{sns}}^2\sigma _{ns}\sigma _s\right)`$ (21)
where

$`\begin{array}{ccccccc}m_{\eta _{ns}}^2\hfill & =& \xi 2\beta b+\lambda a^2,\hfill & & m_{\sigma _{ns}}^2\hfill & =& \xi +2\beta b+3\lambda a^2+4\lambda ^{}a^2,\hfill \\ m_{\eta _s}^2\hfill & =& \xi +\lambda b^2,\hfill & & m_{\sigma _s}^2\hfill & =& \xi +3\lambda b^2+2\lambda ^{}b^2,\hfill \\ m_{\eta _{sns}}^2\hfill & =& 2\sqrt{2}\beta a,\hfill & & m_{\sigma _{sns}}^2\hfill & =& 2\sqrt{2}\left(\beta +\lambda ^{}b\right)a.\hfill \end{array}`$ (25)
Many of our initial concerns can be qualitatively understood from Eqs. (8,14,25). The first point to be noticed is that, in this model, the $`U_A(1)`$ anomaly gets coupled to the v.e.v.’s of scalar fields by the spontaneous breaking of chiral symmetry and contributes, via this effect, to the masses of all fields entering the theory except the strange fields . In the remaining of the paper we will call this the “anomaly-vacuum” effect. The mass terms arising as a consequence of this effect, have exactly the same origin as the mass terms coming from the $`\varphi ^4`$ terms in the $`SU(2)\times SU(2)`$ Linear Sigma Model of Gell-Mann and Levy. The novelty here is the existence of the three-meson determinantal instanton induced interaction which couples either two pseudoscalars to one scalar or three scalars (other possibilities being forbidden by parity). In the case when one of the scalars in the vertex has the same quantum numbers as the vacuum (the isospin singlets $`\frac{1}{\sqrt{2}}(\overline{u}u+\overline{d}d)`$ or $`\overline{s}s`$ ) a mass term is generated by replacing this leg by the corresponding vacuum expectation value. In particular the vertex with a strange pseudoscalar a non-strange pseudoscalar and an isosinglet scalar field (which the determinantal structure dictates to be the non-strange isosinglet) gives rise to a mass term proportional to the strength of the anomaly, which mixes the strange and non-strange pseudoscalar fields as shown in Eq.(refmixed2. Similar phenomena occur in the scalar sector. The strength of the mixing due to the anomaly is exactly of the same size but opposite sign in the pseudoscalar and scalar sectors. In the latter case, there exist an additional contribution coming from the $`\varphi ^4`$ interaction in the Lagrangian Eq. (3) whose strength is measured by $`\lambda ^{}`$.

The extraction of the parameters of the model has been done in  using well known information on the pseudoscalar sector. In this concern, it is worth mentioning that the outcome of the model strongly depends on the input used. In Ref.  pseudoscalar meson masses ($`m_\pi ,m_\eta ,m_\eta ^{},m_K`$) and the pion decay constant ($`f_\pi `$) were used in order to fix the parameters entering the model ($`\xi ,\lambda ,\beta ,a,x=(ba)/2a`$). Ref.  uses $`m_\pi ,m_K,m_\eta ^2+m_\eta ^{}^2,f_K`$ and $`f_\pi `$ as input while Ref.  uses the pseudoscalar mixing angle in the singlet-octet basis $`\theta _P`$ and $`m_\pi ,m_K,m_\eta ,m_\eta ^{}`$. The outcome is different in all these three cases. In particular, the last approach yields a heavy scalar nonet. It must be pointed out, however, that according to recent work , the proper description of pseudoscalar mixing in the singlet-octet basis requires two mixing angles. In this concern, the use of the strange-non-strange basis is more appropriate since, in this basis, pseudoscalar mixing can be described using a single angle . As we shall see below, the pseudoscalar spectrum is consistent, within the model, with a small mixing angle in the singlet-octet basis. An explanation for the physics behind this quantity requires the improvement of the model. The mixing angle of pseudoscalars in the flavor basis($`\varphi _P`$) and its scalar analogous ($`\varphi _S`$) can be extracted from (20) by diagonalizing the Lagrangian. A straightforward calculation yields

$$\mathrm{sin}2\varphi _P=\frac{2m_{\eta _{sns}}^2}{m_\eta ^{}^2m_\eta ^2},\mathrm{sin}2\varphi _S=\frac{2m_{\sigma _{sns}}^2}{m_\sigma ^{}^2m_\sigma ^2}$$
(26)
where $`\sigma ^{}`$ and $`\sigma `$ denote the physical isoscalar-scalar mesons. Alternatively,

$$\mathrm{cos}2\varphi _P=\frac{m_{\eta _s}^2m_{\eta _{ns}}^2}{m_\eta ^{}^2m_\eta ^2},\mathrm{cos}2\varphi _S=\frac{m_{\sigma _s}^2m_{\sigma _{ns}}^2}{m_\sigma ^{}^2m_\sigma ^2}.$$
(27)
The important point is that the parameter which measures the strength of the $`U_A(1)`$ breaking, $`\beta `$, turns out to be negative with the conventions in Eq.(5): $`\beta 1.5GeV`$. The value of $`a`$ can be fixed from $`a=f_\pi /\sqrt{2}=65.7MeV`$, whereas the value of $`b`$, which we rewrite in terms of $`x=(ba)/2a`$ can be extracted by the procedure used in  which yields $`x=0.22`$ or from the input in  which gives $`x=0.39`$. The values of $`\lambda `$ and $`\lambda ^{}`$ are positive. Their actual values depend on the input used but in general $`\lambda ^{}`$ turns out to be small . Using the approach of Ref. , Eqs. (26,27) yield

$$\mathrm{sin}2\varphi _P=0.9202\mathrm{cos}2\varphi _P=0.3911,$$
(28)
which imply a value $`\varphi _P56.7^{}`$ for the pseudoscalar mixing angle in the $`sns`$ basis. This corresponds to a small angle in the singlet-octet basis ($`\theta _P=\varphi _P54.7^{}+2^{}`$) consistent with results in  ($`\theta 5^{}`$) <sup>*</sup><sup>*</sup>*In Ref.  the pseudoscalar mixing angle in the $`sns`$ basis was extracted from the sine relation Eq.(26) and estimated as $`\varphi _P33.3^{}`$. A careful analysis of the diagonalization process shows that the correct value arises from the cosine relation and its actual value is $`\varphi _P56.7^{}`$. The problem when using the sine relation is that it does not distinguish between $`\varphi _P`$ and $`\frac{\pi }{2}\varphi _P`$ which is also a solution. We appreciate illuminating correspondence with Prof. G. ’t Hooft which helped to clarify the input scheme dependence and lead us to reconsider the extraction of the pseudoscalar mixing angle.. The difference in these approaches comes from the use of $`f_K`$ as input in , instead of the combination $`m_\eta ^{}^2m_\eta ^2`$ used in .

The extraction of the scalar mixing angle requires to fix the coupling $`\lambda ^{}`$ which enters the pseudoscalar and the unmixed scalar sectors only in the combination $`\xi =\mu ^2+\lambda ^{}(2a^2+b^2)`$ as can be seen from Eqs.(8,25). Thus, fixing this parameter necessarily requires to use information on the mixed scalar sector. In  the masses of physical scalar mesons where studied as a function of this parameter. This analysis leads to the identification of the isoscalar-scalar mesons with the $`\sigma (400600)`$ and $`f_0(980)`$. This identification yield $`\lambda ^{}4`$ which predicts a scalar mixing angle $`\varphi _S14^{}`$. This result is consistent with the analysis of $`f_0\gamma \gamma `$  and the recently measured $`\varphi \pi ^0\pi ^0\gamma `$ . The small scalar mixing angle is consistent with a mostly $`\overline{s}s`$ structure for the $`f_0(980)`$ meson and a nearby $`(\overline{u}u+\overline{d}d)/\sqrt{2}`$ for the sigma meson. The isovector and isospinor scalar fields are identified with the $`a_0(980)`$ and $`\kappa (900)`$ in the model ( the procedure followed in  identifies $`\kappa `$ with $`K_0^{}(1430)`$, this is a consequence of using different input and to the high dependence of the outcome on the input scheme mentioned above).

## III Pseudoscalar spectrum

The purpose of this letter is to remark the role played by the $`U_A(1)`$ anomaly in structuring the mass spectrum of scalar (and pseudoscalar) mesons. To this end, let us qualitatively analyze the splittings arising from from Eqs. (8,25) in order to disentangle the different contributions to the meson masses. In this concern, it is important to recall that the extraction of the parameters yields a negative sign for $`\beta `$ in Eq.(5) which quantify the strength of the anomaly. From Eqs.(8,25) we obtain the relations

$`m_{\eta _{ns}}^2m_\pi ^2`$ $`=`$ $`4\beta b,`$ (29)
$`m_{\eta _{ns}}^2m_K^2`$ $`=`$ $`2\beta \left(a+b\right)\lambda b\left(ba\right),`$ (30)
which reveal that, in the absence of the anomaly, the pion and non-strange eta are degenerate and the kaon is heavier than both of them, which is consistent with the expectations of naive quark models. Thus, the $`\eta _{ns}\pi `$ splitting is due to the “anomaly-vacuum” effect which pushes the pion and kaon masses down and the non-strange eta mass up as can be seen from Eqs. (8,25). As to the $`\eta _{ns}K`$ splitting there is an additional contribution coming from the invariant in (3) whose strength is measured by the $`\lambda `$ coupling constant. From now on we will call such kind of contributions as the “normal” effects. They are “normal” in the sense that, in the absence of the anomaly, the splittings between pseudoscalar mesons are proportional to this coupling constant times the SU(3) symmetry breaking factor $`(ba)`$ which can be related to the difference of strange and non-strange quark masses. In the case at hand, the “normal” effect tends to make $`K`$ heavier than $`\eta _{ns}`$, but the “anomaly-vacuum” effect goes in the opposite direction. The latter effect is stronger than the former thus rendering the non-strange eta heavier than the kaon.

Similar results are obtained for the $`\eta _s\pi `$ and $`\eta _sK`$ splittings as can be seen from the following relations

$`m_{\eta _s}^2m_\pi ^2`$ $`=`$ $`2\beta b+\lambda \left(b+a\right)\left(ba\right)`$ (31)
$`m_{\eta _s}^2m_K^2`$ $`=`$ $`2\beta a+\lambda a\left(ba\right).`$ (32)
In this case, both effects interfere constructively, reinforcing the corresponding splittings. The individual effects of the anomaly can be read from Eqs.(8,25). The anomaly leaves the strange eta untouched while, as noticed above pushes the pion and kaon masses down.

The next interesting effect in the pseudoscalar spectrum has to do with the $`K\pi `$ splitting. From the relation

$$m_K^2m_\pi ^2=\left(2\beta +\lambda b\right)\left(ba\right),$$
(33)
we see that $`K\pi `$ mass splitting is proportional to the $`SU(3)`$ symmetry breaking. It must be noticed , however, that the “normal” effect is enlarged by the “anomaly-vacuum” effect thus making pions much lighter than kaons. Finally the $`\eta _s\eta _{ns}`$ splitting is

$$m_{\eta _s}^2m_{\eta _{ns}}^2=2\beta b+\lambda \left(b+a\right)\left(ba\right).$$
(34)
In the absence of the anomaly the normal pattern appears, i.e. the strange field is heavier than the non-strange field. Turning on the anomaly modifies this pattern. The anomaly does not affect the strange pseudoscalar mass but it does push up the non-strange pseudoscalar. If we consider the mixing angle predicted by the model, the later turns out to be heavier than the former.

## IV Scalar spectrum

Let us now turn to the scalar sector. The first point to be emphasized is that, the “anomaly-vacuum” contribution systematically has the opposite sign in the pseudoscalar and scalar sectors ( we will call this the “sign”effect ). The reader can be easily convinced of this by looking at the mass relations in Eqs.(8,25). There are two important consequences of the “sign” effect. The first one is that the $`a_0`$ and $`\kappa `$ masses are pushed up by the anomaly (in contrast to the pseudoscalar analogous which are pushed down) and the non-strange sigma is pushed down (unlike the non-strange pseudoscalar which is pushed up) . The second effect has to do with the mixing of strange and non-strange isoscalar scalar fields. The “anomaly-vacuum” effect in the mixing of scalars is exactly of the same size as in pseudoscalars but with the opposite sign . In this case, there exists an additional contribution coming from one of the chiral invariants in the Lagrangian, whose strength is measured by the $`\lambda ^{}`$ coupling. As discussed in  this term corresponds to disconnected quark diagrams, thus being suppressed by the Okubo-Zweig-Izuka rule. This term, however, becomes relevant to the mixing of scalars since it has the opposite sign to the anomaly and interferes destructively with it, thus rendering the scalars less strongly mixed than pseudoscalars.

The $`\sigma _{ns}a`$ and $`\sigma _{ns}\kappa `$ splittings are given by

$`m_{\sigma _{ns}}^2m_a^2`$ $`=`$ $`4\beta b+4\lambda ^{}a^2`$ (35)
$`m_{\sigma _{ns}}^2m_\kappa ^2`$ $`=`$ $`2\beta \left(a+b\right)\lambda \left(2a+b\right)\left(ba\right)+4\lambda ^{}a^2.`$ (36)
Notice that in the absence of the anomaly and the OZI forbidden $`\lambda ^{}`$ coupling the $`a`$ and $`\sigma `$ fields have the same mass and the $`\kappa `$ meson is slightly heavier as expected from their constituent quark content. The individual effects of the anomaly mentioned above (the anomaly pushes the sigma mass down and the $`a`$ and kappa masses up) causes the $`a\sigma `$ splitting making the $`a`$ field heavy and the sigma field light. In other words, the $`\sigma `$ meson is light as a consequence of the “anomaly-vacuum” and the “sign” effects. In the case of the $`\kappa \sigma `$ splitting, these mesons are pushed in opposite directions by the anomaly thus reinforcing the “normal” pattern. In addition we have to take into account $`\lambda ^{}`$ contribution in both cases.

The corresponding relations for $`\sigma _s`$ are

$`m_{\sigma _s}^2m_a^2`$ $`=`$ $`2\beta b+3\lambda \left(a+b\right)\left(ba\right)+2\lambda ^{}b^2`$ (37)
$`m_{\sigma _s}^2m_\kappa ^2`$ $`=`$ $`2\beta a+\lambda \left(a+2b\right)\left(ba\right)+2\lambda ^{}b^2.`$ (38)
Turning off the anomaly (and the $`\lambda ^{}`$ contribution ) yields a pattern as expected from the naive constituent quark picture. The strange scalar is heavier than the non-strange one and the kappa meson. It worth noticing that the strange scalar field is not affected by the “anomaly-vacuum” effect (see Eqs.(25)), but the chiral partner of the pion, the $`a`$ field, is pushed up by the anomaly rendering it almost degenerate with the physical (mostly strange) $`f_0`$ meson. Thus the $`a_0f_0`$ degeneracy is accidental and has its origin in the instanton induced quark interaction. In both equations above, the “normal” effect and the $`\lambda ^{}`$ contributions to the splittings have the same sign. These “normal” splittings are largely canceled by the “anomaly-vacuum”effect. This explains the $`a_0f_0`$ degeneracy and the close value of the $`\kappa `$ meson mass.

The $`a_0\kappa `$ splitting has the same structure as the corresponding pseudoscalar relation in Eq. (33), namely

$$m_\kappa ^2m_a^2=\left(2\beta +\lambda \left(2a+b\right)\right)\left(ba\right)$$
(39)
but in this case the individual“anomaly-vacuum” effects go in the opposite direction due to the “sign”effect. Both meson masses are pushed up by the anomaly with different strengths as the anomaly is coupled to different v.e.v.’s . The result is a not so strong anomaly effect in the splitting which nevertheless is strong enough to largely cancel the “normal” splittings due to quark masses making these mesons roughly equally heavy.

Finally, the $`\sigma _s\sigma _{ns}`$ splitting can be read from Eqs.(25) as

$$m_{\sigma _s}^2m_{\sigma _{ns}}^2=2\beta b+3\lambda \left(b+a\right)\left(ba\right)+2\lambda ^{}\left(b^22a^2\right).$$
(40)
In case of a vanishing anomaly, we obtain the “normal”pattern again. Switching on the anomaly affect the non-strange $`\sigma `$ mass only ( see Eq.(25 )) making this meson light. In other words, the “anomaly-vacuum” and “normal”contributions interfere constructively in the splitting. This makes the non-strange $`\sigma `$ much lighter than the strange sigma. The naively expected pattern is reinforced by the “anomaly-vacuum” effect in contrast to the pseudoscalar analogous case in Eq.(34) where the anomaly overcome the “normal” effect.

The qualitative analysis above can be graphically described as in Figs. 1-2, where the splittings due to the “anomaly-vacuum” and $`SU(3)`$ breaking effects are shown.

Fig.1. Effects of $`U_A(1)`$ symmetry breaking in the pseudoscalar sector.

Fig.2. Effects of $`U_A(1)`$ symmetry breaking in the scalar sector.

## V  Summary

Summarizing, we have obtained a qualitative (and quantitative) understanding of the lowest lying scalar (and pseudoscalar) meson spectrum on the basis of a simple model which incorporates the most relevant properties of QCD in the corresponding energy region, namely, the appropriate degrees of freedom, spontaneous breaking of chiral symmetry and $`U(1)_A`$ symmetry breaking. Most of the unusual properties of both sectors are explained by the coupling of the $`U(1)_A`$ anomaly to the vacuum expectation values of scalars by the spontaneous breaking of chiral symmetry (“anomaly-vacuum” effect ) which supplies mesons with not yet considered mass terms.

Certainly this is only a first step in the elucidation of the quark structure of scalar mesons, and we have considered the $`U_A(1)`$ anomaly as the only striking effect (beyond chiral symmetry breaking). It still remain to conciliate this point of view with the 1/N perspective  which is a main task. The fact that so many phenomena are explained by a single effect makes this model attractive as a starting point for more elaborate and systematic analysis. It is also reassuring that quantitative estimates of the model predictions for the scalar meson decays where intermediate vector mesons can not appear or their contribution is known to be small are successful. Recently, it has been shown that the puzzling $`a_0(980)\gamma \gamma `$ and $`f_0(980)\gamma \gamma `$ decays are properly described within the model  and the recent measurements for the $`\varphi \pi ^0\pi ^0\gamma `$ decay  are consistent with model predictions .

The improvement of the model is necessary in order to reach a more complete description of meson properties. This is shown by its failure to satisfactorily describe the pseudoscalar mixing angle. Other effects such as the hadronic loops considered in  could remedy this. Another possibility is to include degrees of freedom which could be relevant to the physics in this energy region and so far have not been considered. In particular, glueballs degrees of freedom could have some effects via their mixing with pure quarkonium states. A step in this direction is the work in  where glueball degrees of freedom are considered in a framework close to the one used here and its worthy to analyze its consequences in detail. A different line of action which is worth exploring is the relation of the instanton induced quark interaction with the recently raised group theoretical arguments in favor of the unitary symmetry as an accidental symmetry due to the gluon anomaly .


text:
# Star Formation Rate from Dust Infrared Emission

## 1 Introduction

Dust infrared (IR) luminosity111In this Letter, the term ’IR’ means the whole wavelength range of dust emission. The ’IR luminosity’ dose not include the component emitted by stars directly. is a familiar indicator of the current star formation rate (SFR). Especially, the most active star-forming galaxies (e.g., Arp 220) emit almost all their radiation energy in the IR range (Soifer, Houck, & Neugebauer, 1987; Sanders & Mirabel, 1996). For such galaxies, Kennicutt (1998a) derived a conversion formula from the IR luminosity to the SFR by assuming the IR luminosity equal to the bolometric luminosity and 10–100 Myr continuous star formation (hereafter K98 IR formula). More recently, Inoue, Hirashita, & Kamaya (2000; hereafter IHK00) theoretically constructed a new conversion law from IR to SFR, which we can apply not only to dusty starbursts but also to normal (or quiescent) galaxies.

Although the SFRs determined from various indicators should agree with each other, we often find discrepancies among them. Thus, many authors try to resolve such discrepancies (Sullivan et al., 2000; Bell & Kennicutt, 2001; Hopkins et al., 2001; Charlot et al., 2002; Buat et al., 2002; Rosa-González et al., 2002). In this context, the SFR from the IR luminosity via K98 IR formula is assumed to be the correct SFR in spite of the fact that K98 IR formula can be applied to only dusty starbursts in principle (Hopkins et al., 2001; Rosa-González et al., 2002). Indeed, the SFR via K98 IR formula shows a good agreement with that from the H$`\alpha `$ luminosity even for samples of normal galaxies if we correct the H$`\alpha `$ luminosity for the dust extinction in the interstellar medium (ISM), that in H ii regions, and the stellar Balmer absorption (Charlot et al., 2002; Rosa-González et al., 2002). Also, an empirical conversion factor into the SFR from the IR luminosity for normal late-type spiral galaxies agrees with that of K98 IR formula within a factor of 2 (Buat & Xu, 1996). Why can K98 IR formula present the correct SFR even for normal galaxies?

In this Letter, we examine what types of galaxies K98 IR formula is applicable to, and how much correction factor against the formula is required to obtain a more precise SFR. In §2, we formulate the relation between the IR luminosity and the star formation (SF) luminosity. In §3, we describe a model of dust extinction of galaxies. In §4, we show that the ratio of the IR luminosity to the SF luminosity becomes of order unity within a factor of 2 for a very wide parameter area, and many nearby galaxies really fall in this area.

## 2 Infrared Luminosity of Galaxies

Following IHK00, we formulate the relation between the observed IR luminosity of a galaxy, $`L_{\mathrm{IR}}`$, and the intrinsic bolometric luminosity of the newly ($``$ 10 Myr) formed stellar population (hereafter SF population), $`L_{\mathrm{SF}}`$.

The IR luminosity of a galaxy is expressed by the following:

$$L_{\mathrm{IR}}=L_{\mathrm{Ly}\alpha }+(1f)L_{\mathrm{LC}}+ϵL_{\mathrm{UV}}+\eta L_{\mathrm{under}},$$
(1)
where $`L_{\mathrm{Ly}\alpha }`$, $`L_{\mathrm{LC}}`$, $`L_{\mathrm{UV}}`$, and $`L_{\mathrm{under}}`$ are luminosities of Lyman $`\alpha `$ photons generated in H ii regions, Lyman continuum (LC) photons and nonionizing ($`\lambda >912`$ Å) photons emitted by the SF population, and photons from underlying stellar population, respectively. The parameters, $`f`$, $`ϵ`$, and $`\eta `$ are a fraction of LC photons contributing to hydrogen ionization, an efficiency of the dust extinction for nonionizing photons from the SF population, and that for photons from underlying population. Here, we assume that all Lyman $`\alpha `$ photons are eventually absorbed by dust during a number of resonance scatterings in the ISM. We neglect the escape of LC photons from a galaxy, the consumption of LC photons by Helium ionization, any extinctions for hydrogen’s recombination lines within H ii regions except for Lyman $`\alpha `$, and the contribution to the IR luminosity by AGN activity.

According to the Kurucz ATLAS 9 stellar spectra with Solar metallicity and a turbulence speed of 2 km s-1, we obtain $`L_{\mathrm{LC}}=0.37L_{\mathrm{SF}}`$ and $`L_{\mathrm{UV}}=0.63L_{\mathrm{SF}}`$ if the Salpeter IMF (0.1–100 $`M_{\mathrm{}}`$) is assumed. Since about two-thirds of the ionization-recombination processes produce Lyman $`\alpha `$ photons under case B (Spitzer, 1978), we also obtain $`L_{\mathrm{Ly}\alpha }=0.12fL_{\mathrm{SF}}`$. Here, we note that only a fraction $`f`$ of LC photons ionize hydrogen. Moreover, we introduce a new parameter, $`\gamma `$, to represent the intrinsic stellar luminosity fraction of the SF population to all stellar population; $`\gamma =L_{\mathrm{SF}}/L_{\mathrm{total}}`$, where $`L_{\mathrm{total}}=L_{\mathrm{SF}}+L_{\mathrm{under}}`$. This parameter is an indicator of the star-forming activity. Therefore, equation (1) is reduced to

$$\frac{L_{\mathrm{IR}}}{L_{\mathrm{SF}}}=0.370.25f+0.63ϵ+\eta \frac{1\gamma }{\gamma }.$$
(2)
For dusty starbursts ($`\gamma 1`$), the right hand side of equation (2) becomes unity since the dust opacity is enough high to absorb all photons emitted by stars (i.e. $`f=0`$ and $`ϵ=1`$). In principle, K98 IR formula can be applied to only such galaxies. On the other hand, if we determine all four parameters ($`f`$, $`ϵ`$, $`\eta `$, and $`\gamma `$), we can convert the observed IR luminosity of all kinds of galaxies to the SF luminosity, and then, to the current SFR. This is the algorithm of IHK00.

## 3 Dust Opacity of Galaxies

Here, we describe a model of dust extinction to determine the values of $`f`$, $`ϵ`$, and $`\eta `$, in equation (2). We will discuss in §4.2 the effect of the choice of a specific extinction law on our result.

Following Charlot & Fall (2000) (especially their Fig. 1), we consider, first, two types of dust extinction. One is the extinction by diffuse ISM dust ($`\tau _\lambda ^{\mathrm{ISM}}`$), and the other is by dust accompanying the birth clouds of young stars ($`\tau _\lambda ^{\mathrm{BC}}`$). Then, we assume $`\tau _\lambda ^{\mathrm{ISM}}=\tau _V^{\mathrm{ISM}}(\lambda /\lambda _V)^n`$ and $`\tau _\lambda ^{\mathrm{BC}}=\mu \tau _V^{\mathrm{ISM}}(\lambda /\lambda _V)^n`$ for $`\lambda >912`$ Å, where $`\mu =2`$, $`n=0.7`$, and $`\lambda _V=5500`$ Å are adopted, so as to reproduce the observed relation between UV flux slope and luminosity ratio of IR to UV for nearby star-forming galaxies (Charlot & Fall, 2000). We define the dust optical depth for the SF population as $`\tau _\lambda ^{\mathrm{SF}}=\tau _\lambda ^{\mathrm{BC}}+\tau _\lambda ^{\mathrm{ISM}}`$.

For simplicity, we assume that the dust geometry is an uniform screen (or the optical depth defined above is an effective one). Then, we define $`ϵ`$ as an average of $`ϵ(t)`$ over the age range of 0 – 10 Myr, where $`ϵ(t)=1L_\lambda (t)\mathrm{exp}(\tau _\lambda ^{\mathrm{SF}})𝑑\lambda /L_\lambda (t)𝑑\lambda `$ and $`L_\lambda (t)`$ is the luminosity density of a simple stellar population (SSP) with an age $`t`$. We use SSPs produced by PÉGASE 2.0 (Fioc & Rocca-Volmerange, 1997) throughout this Letter. In the same way, we define $`\eta `$ as an average of $`\eta (t)`$ over the age of 10 Myr – 15 Gyr but $`\tau _\lambda ^{\mathrm{SF}}`$ is replaced with $`\tau _\lambda ^{\mathrm{ISM}}`$. The obtained $`ϵ`$ and $`\eta `$ are a function of only $`\tau _V^{\mathrm{ISM}}`$, the normalization of the adopted extinction law.

Let us introduce one more type of dust extinction to determine the parameter $`f`$, Lyman continuum extinction (Inoue, Hirashita, & Kamaya 2001, Inoue 2001), which is the extinction for LC photons by dust in H ii regions before LC photons ionize neutral hydrogens. The parameter $`f`$ is determined from the dust optical depth for LC photons in H ii regions, $`\tau _{LC}^{\mathrm{HII}}`$, via eq. (8) in Petrosian, Silk, & Field (1972). Here, we assume a scaling law between $`\tau _{LC}^{\mathrm{HII}}`$ and $`\tau _V^{\mathrm{ISM}}`$, i.e., $`\tau _{LC}^{\mathrm{HII}}/\tau _V^{\mathrm{ISM}}\xi `$.

Now, we attempt to estimate $`\xi `$ from real data. We have shown that $`\tau _{LC}^{\mathrm{HII}}`$ is proportional to the dust-to-gas mass ratio, $`𝒟`$, in Inoue et al. (2001). Here, we adopt the relation, $`\tau _{LC}^{\mathrm{HII}}(𝒟/𝒟_{\mathrm{MW}})`$ (Inoue, 2001), where $`𝒟_{\mathrm{MW}}=6\times 10^3`$ is a typical Galactic value (Spitzer, 1978). We estimated $`𝒟`$ for seven nearby spiral galaxies observed by ISO (Alton et al., 1998) in Inoue et al. (2001). The mean value is $`1.8\times 10^3`$.222The mean value is different from that in Inoue et al. (2001) because we correct it for He contribution and for the calibration error of the flux density at 200 µm (-30%; Alton et al. 1998). In addition, we should note that $`𝒟_{\mathrm{MW}}`$ by Spitzer (1978) is estimated from the analysis of the extinction properties for a number of line of sights, which is different from our method. On the other hand, we can estimate the dust optical depth from the IR surface brightness. By neglecting the self-absorption of the IR photons and assuming the obtained optical depth to correspond to that in the diffuse ISM, we obtain $`\tau _V^{\mathrm{ISM}}=0.36`$ as a mean dust optical depth for these seven galaxies, where we adopt the dust emissivity determined by Bianchi, Davies, & Alton (1999). Therefore, we obtain $`\xi 1`$.

Although the real relation between $`\tau _{LC}^{\mathrm{HII}}`$ and $`\tau _V^{\mathrm{ISM}}`$ is quite uncertain, we find in equation (2) that the effect of $`f`$ on the ratio $`L_{\mathrm{IR}}/L_{\mathrm{SF}}`$ is small. Indeed, the effect of the deviation from $`f=0.5`$ on the $`L_{\mathrm{IR}}/L_{\mathrm{SF}}`$ is at most $`\pm 0.125`$.

## 4 Luminosity Ratio of Dust Infrared Emission to Star-Forming Population

Now, we can determine $`L_{\mathrm{IR}}/L_{\mathrm{SF}}`$ from only two parameters; a mean dust opacity ($`\tau _V^{\mathrm{ISM}}`$) and a star-forming activity ($`\gamma `$). We survey these parameters space to clarify when we can determine the SFR from $`L_{\mathrm{IR}}`$ within a factor of 2 uncertainty via K98 IR formula, i.e., when $`0.5L_{\mathrm{IR}}/L_{\mathrm{SF}}2.0`$ is satisfied.

In Figure 1, we show some constant $`L_{\mathrm{IR}}/L_{\mathrm{SF}}`$ lines in the $`\tau _V^{\mathrm{ISM}}`$-$`\gamma `$ plane. We define the ’applicable area’ as the parameter area satisfying the condition $`0.5L_{\mathrm{IR}}/L_{\mathrm{SF}}2.0`$. The ’applicable area’ is between two thick lines in Figure 1 since the upper thick line is the case of $`L_{\mathrm{IR}}/L_{\mathrm{SF}}=0.5`$ and the lower one is the case of $`L_{\mathrm{IR}}/L_{\mathrm{SF}}=2.0`$. Also, the thin line in Figure 1 represents the case of $`L_{\mathrm{IR}}/L_{\mathrm{SF}}=1.0`$. We find that the ’applicable area’ is very wide in the $`\tau _V^{\mathrm{ISM}}`$-$`\gamma `$ plane. This is a main conclusion in this Letter.

For a galaxy with an exponentially declining star formation history (SFH) of $``$ 10 Gyr time scale, age of $``$ 1 Gyr, and Salpeter IMF, we find that $`\gamma 0.5`$. Thus, many spiral and irregular galaxies are likely to fall in the ’applicable area’. In order to examine whether real galaxies fall in the ’applicable area’ or not, we will estimate their $`\gamma `$ and $`\tau _V^{\mathrm{ISM}}`$ in the following. Then, we will discuss implications from our results.

### 4.1 Estimation of $`\gamma `$ and $`\tau _V^{\mathrm{ISM}}`$

First, we estimate $`\gamma `$ of galaxies from the observed equivalent width of H $`\alpha `$ emission line corrected for dust and stellar absorption effects. By the definition of $`\gamma `$, we can express $`\gamma =EW_{\mathrm{H}\alpha }(C_{\mathrm{H}\alpha }/C_R)`$, where $`EW_{\mathrm{H}\alpha }`$ is the intrinsic equivalent width of H $`\alpha `$ emission line, $`C_{\mathrm{H}\alpha }`$ is the intrinsic ratio of the bolometric luminosity to the H $`\alpha `$ emission line luminosity of the SF population, and $`C_R`$ is the intrinsic ratio of the bolometric luminosity to the luminosity density at $`R`$-band of all stellar population.

If we specify the metallicity and IMF of galaxies, the value of $`C_{\mathrm{H}\alpha }`$ is determined from only the electron temperature and density in H ii regions under case B approximation. On the contrary, the value of $`C_R`$ depends on the SFH and age of galaxies as well as the metallicity and IMF. We examine the evolution of $`C_R`$ with the galactic age, so that we find that the factor $`C_{\mathrm{H}\alpha }/C_R`$ falls in the range of 0.002 – 0.004 unless we consider a galaxy younger than a few Gyr or a galaxy with a short ($``$ 1 Gyr) exponential time scale. Thus, we adopt a constant 0.003 as $`C_{\mathrm{H}\alpha }/C_R`$, which corresponds to the case of Solar metallicity, exponentially declining SFH with 10 Gyr time scale, and about 13 Gyr age. Here, we do not take into account the chemical enrichment effect.

Therefore, we estimate $`\gamma `$ from the following equation approximately;

$$\gamma =\{\begin{array}{cc}0.003EW_{\mathrm{H}\alpha }\hfill & \text{for }EW_{\mathrm{H}\alpha }333\text{ Å}\hfill \\ 1\hfill & \text{for }EW_{\mathrm{H}\alpha }>333\text{ Å}\hfill \end{array},$$
(3)
where the second case is needed due to $`\gamma 1`$ by its definition. In the real case, the value of $`\gamma `$ asymptotically approaches unity by decreasing the factor $`C_{\mathrm{H}\alpha }/C_R`$, when $`EW_{\mathrm{H}\alpha }`$ increases. Thus, we overestimate $`\gamma `$ for a galaxy with a very large $`EW_{\mathrm{H}\alpha }`$. However, the real value of $`\gamma `$ is about unity for such galaxies. On the other hand, we underestimate $`\gamma `$ for a galaxy with a very small $`EW_{\mathrm{H}\alpha }`$. We confirm that this effect is also not so large (at most factor of 2).

Although the effect of the change of the metallicity on $`C_{\mathrm{H}\alpha }/C_R`$ is small, that of the IMF slope is large. If we choose the IMF slope of 1.35 (Salpeter’s slope is 2.35), a typical $`C_{\mathrm{H}\alpha }/C_R`$ decreases an order of magnitude, about 0.0003. In the case of the slope of 3.35, we find $`C_{\mathrm{H}\alpha }/C_R0.007`$ typically. This is because a top-heavy IMF (1.35 slope) results in a smaller fraction of the $`R`$-band luminosity density in the bolometric luminosity since massive stars do not contribute to the $`R`$-band light so effectively. A steep IMF provides an opposite result.

Now, we must estimate the intrinsic equivalent width, $`EW_{\mathrm{H}\alpha }`$, from the observed one, $`EW_{\mathrm{H}\alpha }^{\mathrm{obs}}`$. Assuming that the underlying stellar population (age $``$ 10 Myr) dominates the $`R`$-band continuum and the stellar Balmer absorption, we can express

$$EW_{\mathrm{H}\alpha }=(EW_{\mathrm{H}\alpha }^{\mathrm{obs}}+EW_{\mathrm{H}\alpha }^{,\mathrm{abs}})\frac{1}{f}\mathrm{exp}(\tau _{\lambda _{\mathrm{H}\alpha }}^{\mathrm{SF}}\tau _{\lambda _{\mathrm{H}\alpha }}^{\mathrm{ISM}}),$$
(4)
where $`EW_{\mathrm{H}\alpha }^{,\mathrm{abs}}`$ is the equivalent width of the stellar absorption, and $`\lambda _{\mathrm{H}\alpha }=6563`$ Å. The factor $`1/f`$ means that only a fraction $`f`$ of LC photons is used to ionize neutral hydrogen atoms.

Next, we estimate $`\tau _V^{\mathrm{SF}}`$ from the Balmer decrement corrected for the stellar absorption by assuming $`EW_{\mathrm{H}\alpha }^{,\mathrm{abs}}=EW_{\mathrm{H}\beta }^{,\mathrm{abs}}=5`$ Å, which noted by Kennicutt (1992) as a typical value, and then, we convert the obtained $`\tau _V^{\mathrm{SF}}`$ into $`\tau _V^{\mathrm{ISM}}`$ by using the parameter $`\mu =2`$ (See §3). Here, we also assume the electron temperature of $`10^4`$ K and case B. Once the optical depth $`\tau _V^{\mathrm{ISM}}`$ is obtained, we can estimate $`f`$ parameter, $`EW_{\mathrm{H}\alpha }`$ from $`EW_{\mathrm{H}\alpha }^{\mathrm{obs}}`$, and then, $`\gamma `$. The obtained values of $`\tau _V^{\mathrm{ISM}}`$ are systematically smaller than but roughly consistent with those determined from the UV slope method (Meurer, Heckman, & Calzetti, 1999) for sample galaxies described below.

The uncertainty of $`\tau _V^{\mathrm{ISM}}`$ from the Balmer decrement is large, $`\pm 0.2`$. As a result, the uncertainty of the obtained $`\gamma `$ is $`\pm 0.2`$, which is estimated from a simple Monte Carlo simulation if we do not take account of the uncertainties of the adopted parameters of the dust extinction law (i.e., $`\eta `$, $`\xi `$, and $`n`$ in §3). Also, we find that a random error of $`\tau _V^{\mathrm{ISM}}`$ causes a somewhat overestimation of $`\gamma `$. To display this point, we show a typical error bar in Figure 1. In addition, we select only galaxies with enough strong H $`\beta `$ ’emission’ line to obtain a reliable Balmer decrement. Thus, our sample galaxies are biased towards larger $`\gamma `$.

One might think that the parameter $`\gamma `$ can be estimated from the ratio of the UV luminosity and the $`R`$ or $`K`$-band luminosity. However, we can do it reasonably only if we specify the SFH and age of galaxies. That is, the coefficient between $`\gamma `$ and $`L_{UV}/L_R`$ or $`L_{UV}/L_K`$ varies together with the galactic age significantly. This is because a significant fraction of the UV light originates in the underlying stellar population with its age of $``$ 10 Myr. On the contrary, the ionizing photons producing the H $`\alpha `$ photons originate in almost only the SF population (its age $``$ 10 Myr). Therefore, we estimated $`\gamma `$ from the equivalent width of the H $`\alpha `$ line.

### 4.2 Results & Discussions

To plot real galaxies in Figure 1, we compile all morphological types of galaxies except for elliptical galaxies from the literature. In Figure 1, the filled circles are early type spiral (Sa–Sab) galaxies observed by Usui, Saitō, & Tomita (2001), who selected active star-forming early type spiral galaxies with $`\mathrm{log}(L_{\mathrm{FIR}}/L_B)0.5`$. The filled squares, open squares, and crosses are early type spiral (S0–Sbc), late type spiral (Sc–Sm), and irregular (including peculiar and merger) galaxies, respectively, observed by Kennicutt (1992), who selected galaxies covering all morphological types. The filled triangles, open triangles, and stars are also early spiral, late spiral, and irregular (also peculiar and merger) galaxies, respectively, observed by Storchi-Bergmann, Kinney, & Challis (1995), who selected active star-forming galaxies observed by $`IUE`$ satellite. Any galaxies containing an AGN are excluded.

We find that almost all (43/52) galaxies fall in the ’applicable area’. For 7 out of 9 galaxies in the outside of the ’applicable area’, we fail to determine $`\tau _V^{\mathrm{ISM}}`$ because they become negative values. Although the uncertainties of the determined $`\tau _V^{\mathrm{ISM}}`$ and $`\gamma `$ are large and compiled galaxies do not make a complete sample, this result indicate that many nearby normal galaxies are really in the ’applicable area’.

We emphasize the following two points. First, it is a just coincidence that K98 IR formula is applicable to normal galaxies. It results from offsetting two effects of a small dust opacity and a large cirrus contribution of normal galaxies relative to starburst galaxies on the conversion coefficient of K98 IR formula (Kennicutt 1998b, IHK00). For an example case of a galaxy with $`\tau _V^{\mathrm{ISM}}=0.3`$ and $`EW_{\mathrm{H}\alpha }^{\mathrm{obs}}=50`$ Å, while only about 70% luminosity of the SF population is converted into the IR luminosity, it is almost the same amount as the rest 30% SF luminosity that is supplied with the IR luminosity by the underlying population.

Second, the result above is robust for changing the adopted parameter sets of the dust extinction. One might think that our result depends on the dust extinction law adopted. However, we note that it dose not change when we change the parameters in the range of $`0\mu 3`$, $`0.4n1.0`$, and $`0\xi 2`$ as shown in Figure 2 (an extra figure). This point will be discussed in detail in the forthcoming paper (Inoue 2002, in preparation).333Visit the following web cite where you will find figures showing the cases of other parameter choice; http://www.kusastro.kyoto-u.ac.jp/~inoue/irsfr.html.

We discuss the condition required to fall in the ’applicable area’. From the lower bound of $`\gamma `$ in Figure 1 ($`L_{\mathrm{IR}}/L_{\mathrm{SF}}=2.0`$), we find that a quiescent galaxy with a large dust amount (a small $`\gamma `$ and a large $`\tau _V^{\mathrm{ISM}}`$) is in the outside of the ’applicable area’. Since galaxies showing a lower activity of the star formation tend to be more transparent (Wang & Heckman, 1996; Hopkins et al., 2001), the number of galaxies falling in the area of large $`\tau _V^{\mathrm{ISM}}`$ and small $`\gamma `$ may be small. Moreover, we find that the lower bound of $`\gamma `$ is $`0.1`$ for a galaxy with $`\tau _V^{\mathrm{ISM}}0.3`$. This corresponds to $`EW_{\mathrm{H}\alpha }^{\mathrm{obs}}10`$ Å. This boundary can be exceeded not only by active star-forming galaxies but also by many nearby normal spiral and irregular galaxies (Kennicutt, 1983).

The upper bound of $`\gamma `$ in Figure 1 ($`L_{\mathrm{IR}}/L_{\mathrm{SF}}=0.5`$) reaches unity if $`\tau _V^{\mathrm{ISM}}`$ is larger than about 0.1. Thus, a starburst galaxy having very small dust content of $`\tau _V^{\mathrm{ISM}}0.1`$ is in the outside of the ’applicable area’. For example, I Zw 18, which is an active star-forming and the most metal deficient nearby dwarf galaxies, has typically $`A_V=0.2`$ mag in the southeast H ii region (Cannon et al., 2002). It corresponds to $`\tau _V^{\mathrm{ISM}}0.06`$ if we adopt $`\mu =2`$. By taking into account $`\gamma 1`$ estimated from very large $`EW_{\mathrm{H}\alpha }^{\mathrm{obs}}1000`$ Å (Cannon et al., 2002), I Zw 18 is in the outside of the ’applicable area’.

In conclusion, we find that the range of a dust opacity and a star-forming activity to be $`0.5L_{\mathrm{IR}}/L_{\mathrm{SF}}2.0`$ is very large. Furthermore, we confirm that these parameters of many nearby normal and active star-forming galaxies really fall in this range. Therefore, we can apply K98 IR formula assuming $`L_{\mathrm{IR}}=L_{\mathrm{SF}}`$ to many nearby normal galaxies as well as dusty starbursts without any correction factor within a factor of 2 uncertainty. However, we cannot apply it to dust rich but quiescent galaxies and extremely dust poor galaxies. Therefore, we need to divide the coefficient of K98 IR formula by a correction factor, $`L_{\mathrm{IR}}/L_{\mathrm{SF}}`$, when we apply it to such galaxies. Following the model in this Letter, we can estimate the correction factor for any galaxies if we have all data required, although we must more refine the model.

The author would like to thank H. Kamaya, T. T. Takeuchi, and H. Hirashita for their suggestions to improve the quality of this work, and also, to thank an anonymous referee for his/her many helpful comments to this paper. The author has made extensive use of NASA’s Astrophysics Data System Abstract Service (ADS).


text:
# ARE FOREST FIRES PREDICTABLE?

## 1 Introduction

It would be of obvious interest if one could predict a forest fire. An attractive aim would be also to evaluate the risk of the forest fires in a given year, having the statistics of the fires in the past. Here we argue that the time dependence of the forest area burned in a given year is inherently chaotic. This means that such predictions are not possible in longer time scale. The method used is the mean field equation based on an appropriate cellular automaton.

Cellular automata (CA) is a modern tool for qualitative simulations of numerous problems in statistical mechanics and theory of dynamical systems . In general, the rules of CA are probabilistic — they may depend on a random variable. When this dependence is absent, the automaton is termed deterministic. These automata are most exciting but least predictable. The behavior of a probabilistic cellular automaton can be predicted successfully  by means of a nonlinear mean field equation. The solution is the time-dependent probability distribution of cell states. The approach is termed “mean field”, because the correlations between the states of neighboring cells are lost within one time step. However, they are lost also in the original probabilistic automaton after several time steps. Similar technique can be applied also to investigate deterministic CA , but in this case the results of the mean field approximation are less reliable: the rules of an original deterministic automaton can preserve the correlations forever.

Our aim here is to apply this mean field technique to the problem of the forest fires. The list of references of the problem is vast and reveals its numerous aspects, from purely computational to inherently biological ones . In particular, some controversy has been mentioned  on the applicability of the percolation theory , and on the size dependence of results of possible experiments. The advantage of the mean field approximation is that the finite-size fluctuations are eliminated. Then, our results can be seen as a limit case both for simulations and experiments. Below we demonstrate that the critical concentration for the percolation problem is fairly reproduced by the forest fire simulation, if the lattice is large enough.

The approach can be formulated as follows: The model forest is a two-dimensional regular lattice. Each cell at the lattice may be occupied by a tree with probability $`p`$. During a year (one time step), trees appear at empty cells with a probability $`rp(1p)`$, where $`p`$ is the previous concentration of trees per cell, and $`0r1`$ is a model control parameter. This nonlinear character of the probability reflects the obvious tendency that a new tree appears less likely if there is almost no free space around the cell. This kind of dependence is known as the Verhulst term . All trees, those existing previously and those just grown, form clusters of different sizes. A cluster is a set of trees which are connected by nearest-neighbor bonds in the von Neumann neighborhood. Now, a probability that a cluster is ignited is proportional to its size. Then, a fraction of the forest burned each summer is just the mean square of the cluster size. We will show that below the percolation threshold $`p_c`$, this fraction is negligible. Above this concentration, most of the forest disappears — only small clusters survive.

Summarizing our method, the time evolution equation for the concentration of trees in the forest is a superposition of two maps, which describe a sequence of growing (“spring”) and fire (“summer”)

$$p_{n+1}=W(p_n+rp_n(1p_n)),$$
(1a)
where the map $`W(p)`$ is
$$W(p)=pA(p)$$
(1b)
and $`n`$ enumerates years, i.e. the time steps. The function $`A(p)`$ is the burned fraction of the total area of the forest, i.e. the mean square of the cluster size. This fraction is evaluated by means of the Hoshen–Kopelman algorithm .

The purpose of this paper is to demonstrate, that the map $`p_np_{n+1}`$ leads to chaotic behavior. We note that such a conclusion is not entirely new. In 1990, Chen, Bak and Jensen  proposed a deterministic coupled map lattice  for a model description of forest fires. In this model, fire was introduced as a continuous process. The Lyapunov exponent calculated there described the time evolution of the Hamming distance between two nearby trajectories. Initial distance was set by a random perturbation, added to each site at time $`t_0`$. As the distance was found to increase with the power law, the value of the Lyapunov exponent was zero. Later on, this model was generalized by Socolar, Grinstein and Jayaprakash  by adding a parameter, which allowed to get a positive Lyapunov exponent.

Both the above mentioned approach and this work have in common that they are computational rather than oriented to experiment. However, within this computational stream there are serious differences. The algorithm of Chen et al.  is deterministic. The forest fire is continuous; a tree burns spontaneously as soon as it is large enough. Our map can be treated as an extrapolation of the output of a probabilistic automaton to the area of infinite size. Then, a sequential iteration is applied, and the season of growing and the season of fires are treated as two successive parts of one step of the iteration. From Ref.  we know that there are many uncontrolled quantities relevant to the problem: moisture, humidity, landscape, etc. As it is impossible to take them into account within a simple scheme, our probabilistic approach seems to be closer to reality. Last but not least, our approach gives a sequence of yearly burned areas which, although expressed in arbitrary units, can be directly compared to available statistical data.

The paper is organized as follows: In Section 2 we describe the way how the Hoshen–Kopelman is used here. There, two different kinds of statistics are assigned to casual (lightnings) and intentional (“human factor”) forest fires. In Section 3 the results on the Lyapunov exponent and the bifurcation diagram are given. A remark is also added on the fire size distribution. The comparison of our results to the statistical data on forest fires in Canada in years 1970–2000  is included in Section 4. The last section is devoted to discussion.

## 2 Numerical approach

With the Hoshen–Kopelman algorithm  we are able to label all occupied sites on a $`L\times L`$ large square lattice in such a way that the sites with the same label belong to the same cluster and different labels are associated to different clusters. In this way we get the cluster size distribution, which allows us to obtain the function $`A(p)`$ and, in the next step, to formulate the mean field map given in Eq. (1).

The beauty of the Hoshen–Kopelman algorithm for the cluster characterization is that it goes through the lattice only once, and stores only the current line of the length $`L`$. Thus it is particularly useful for the bus-bar percolation problem investigation and it allows to check whether a given occupied site at distance $`x`$ from an arbitrary chosen (usually the first and full occupied) line is still connected to this line through the occupied sites in lines $`2,3,\mathrm{},x2,x1`$ .

However, to construct the map (1b) it is necessary to evaluate the mean cluster size (or its average mass) $`L^2A(p)`$. This requires an information on all bonds to sites at the distance $`x`$ from the first line also through the occupied sites at the lines $`x+1,x+2,\mathrm{},L1,L`$, as it is presented on Fig. 1.

This in turn requires storing of the whole lattice for all the time of simulation and passing through the lattice twice — what wastes computer memory and machine time, but ensures proper labeling of the clusters.

The simulation is carried out on $`10^3\times 10^3`$ square lattice, with a fraction of $`p+rp(1p)`$ occupied sites for different values of the initial trees concentration $`0p_01`$. An average concentration of trees in a “green” (but soon “red” and consequently “black”) cluster

$$A(p)=\frac{1}{L^2}\underset{i}{}w_is_i$$

may be evaluated with two kinds of weights

$$w_i^{}=s_i/L^2\text{ and }w_i^{\prime \prime }=s_i/\underset{i}{}s_i,$$

which describe the probability of finding a cluster of size $`s_i`$ on the lattice and the probability of choosing an $`i`$-th cluster among all other clusters, respectively. These statistics can be interpreted as reflecting the probability of a casual fire started from a lightning, which can strike at an empty cell ($`w^{}`$), or of a fire started intentionally at a randomly selected tree ($`w^{\prime \prime }`$). The obtained map (1) for each $`p`$ is the average over a hundred of independent lattice realizations.

Below we concentrate on the case of the fires of natural origin. The results for the other case will be only listed. The details of the latter calculation will be discussed more thoroughly elsewhere.

## 3 Results

### 3.1 The map

The map $`p_np_{n+1}`$ is given in Fig. 2. Two plots refer to two kinds of statistics, described in Section 2.

We see that the plots are the same below the maximum of the curve, which is near the percolation limit for the square lattice $`p_c=0.59273`$ . Actually, below this limit both plots are identical to the curve $`p+rp(1p)`$ from Eq. (1a). This means, that for $`p<p_c`$ the fires are negligible. Above the percolation limit, the curves remain similar for $`p`$ close to $`p_c`$, but for higher $`p`$ the curve for the weight $`w^{}`$ is remarkably higher, than the curve for the weight $`w^{\prime \prime }`$. It is clear that intentional fires ($`w^{\prime \prime }`$) are more devastating, because in this case it is sure that a tree is ignited.

### 3.2 The bifurcation diagram

The rest of the results refer to the case of casual fires, i.e. the weight $`w^{}`$. In Fig. 3 we show the bifurcation diagram. This diagram is obtained with the map $`p_np_{n+1}`$, supplemented by straight lines joining neighboring $`10^3`$ numerical points. We can see the windows of stable cycles of lengths four, eight, five, ten, six, twelve, seven and eight again (Fig. 3), five, three, five and two in Fig. 3, when the parameter $`r`$ decreases. The presence of the cycle three proves that, by means of the Sharkovskii theorem, all possible lengths of the limit cycles are present .

A question arises, if our bifurcation diagram is not an artificial consequence of the finite grid of the map and/or the finiteness of the lattice. To check this point, we calculated also the bifurcation diagram for a piecewisely analytical map $`G(H(p))`$, where

$$\begin{array}{cc}\hfill H(p)& =p+rp(1p),\hfill \\ \hfill G(p)& =\{\begin{array}{cc}0.1\mathrm{exp}[100(0.6p)]+0.4\sqrt{1p}\hfill & p>0.6,\hfill \\ p\hfill & p0.6.\hfill \end{array}\hfill \end{array}$$
(2)
This map is also shown in Fig. 2. It is designed not as to fit the plot obtained within the Hoshen–Kopelman algorithm, but rather to investigate the consequences of the particular shape of this plot. The obtained bifurcation diagram (Fig. 3) shows the windows of stability, with the width which increases with the parameter $`r`$. The lengths of the cycles decreases with $`r`$ in the same way as presented in Fig. 3. We deduce that this part of the bifurcation diagram is generic.

To each cycle, an admissible word can be assigned according to the rules of symbolic dynamics. We have only to check the values of $`p`$ in the sequence of a given limit cycle. The rule is as follows: The value of $`p`$ preceding the largest $`p`$ (close to $`p_c`$) in the cycle is labeled $`C`$. Let us call this value $`p^{}`$. The next value (close to $`p_c`$) is labeled $`R`$. Any other $`p`$ is labeled $`R`$ if it is larger than $`p^{}`$, and it is labeled $`L`$ if its value is smaller than $`p^{}`$ . We use this rule to assign symbolic words to each cycle in Fig. 3, and we compared their sequence to the standard one for the unimodal maps . The sequences of the cycles presented in Figs. 3 and 3 are respectively $`RL^2`$, $`RL^3RL^2`$, $`RL^3`$, $`RL^4RL^3`$, $`RL^4`$, $`RL^5RL^4`$, $`RL^5`$, $`RL^6`$ and $`RL^2R`$, $`RL`$, $`RLR^2`$, $`R`$ with decreasing parameter $`r`$. Note that $`C`$ is omitted for convenience at the beginning of each word. We have found that these sequences indicate opposite directions of the increase of the parameter $`r`$. According to the results discussed in the preceding paragraph, the left part of the bifurcation diagram (Fig. 3) is supposed to vanish in the thermodynamic limit. Still, the conclusion on all possible lengths of the limit cycles holds, as it follows from the presence of cycle five (all lengths but three) and cycle seven (all but three and five) . For the bifurcation diagram in Fig. 3 the above prescription gives words $`RL^n`$. However, in this case the maximum of the map does not appear in the limit cycle, except at the bifurcation points where the cycle length changes. In this sense the letter $`C`$ should not be used.

We have checked that the bifurcation diagram for weigths $`w^{\prime \prime }`$ seems to be a homogeneous spot. Still, we have found the cycle of length three by the symbolic dynamics method . The details will be given elsewhere.

### 3.3 The Lyapunov exponent

In Fig. 4 we show the largest Lyapunov exponent as dependent on the parameter $`r`$. Only positive values are shown, and the spaces between data agree with the windows of stability in the bifurcation diagram. We have found that the negative values of the Lyapunov exponent are much more difficult to be calculated when we have a limit cycle and not a fixed point. In this case, the time dependence of the difference between trajectories in not monotonic. A trajectory can be trapped by a limit cycle with a shifted phase, which leads to a fixing of the difference between trajectories even in the case of a stable limit cycle. Moreover, sometimes this trapping can be only temporal. As the result, in the presence of a stable limit cycle the difference between trajectories can depend on the initial point. These complexities are revealed in Fig. 5.

The Lyapunov exponent for weights $`w^{\prime \prime }`$ seems to be positive in almost the whole range of the parameter $`r`$. This confirms that the windows of stability are much more narrow.

### 3.4 Fire size distribution

We have calculated the histogram of the size of fires, given by the function $`A(p_n)A_n`$. Typical results for the map calculated by Eq. (1) for the weights $`w^{}`$ are given in Fig. 6. As we see, the “fires of size zero” are most frequent. However, larger fires can be more likely that smaller ones. This kind of curves do not show any self-similarity.

A serious difference between the results for $`w^{}`$ and $`w^{\prime \prime }`$ is found for the histogram of the forest fires. An example for $`r=0.3`$ is presented in Fig. 7. For large fires, the curve seems to follow the exponential dependence. To capture this tendency, we have fitted the map for $`w^{\prime \prime }`$ by a curve similar to Eq. (2), where the function $`G(p)`$ is chosen as

$$G(p)=p_c\mathrm{exp}(12.1735\sqrt{pp_c})$$
(3)
if $`p>p_c=0.59237`$. The obtained histogram is added to Fig. 7. Note that the latter curve is monotonic except the point at $`A=0`$, i.e. the events without a fire. As we see, small finite fires are neglected by the above analytic approximation.

## 4 Comparison with experimental data

Our map $`p_np_{n+1}`$ is obtained for the lattice of million sites, as an average over hundred realizations. The accessible statistical data  relate to areas larger than million trees each, but to thirty “iterations” for one initial state of the forest. On the other hand, several uncontrolled factors (moisture, etc.) cannot be taken into account within our simplified approach. We expect that these agents lead to some noise, which is absent in our deterministic map. This noise can be reproduced here by decreasing the model lattice. That is why we are going to compare the experimental data  with the results of our simulation performed for a small lattice.

We have checked that the time dependence of the size of the fires of both experiment and calculations show the same kind of oscillations. Examples are given in Fig. 8. In Fig. 8 the data are presented on the forest fire in Yukon, Canada, in years 1971–1998 .

In Fig. 8 the same data are shown for the simulation performed on the lattice of $`50\times 50`$ cells. As we see, the qualitative features of the data are reflected by the simulation.

We also checked that the data on other provinces in Canada in years 1970–2000  do not differ in any essential way. In Fig. 6 we show a collective histogram for experimental data. For each province or territory, the fires are rescaled to the maximal fire in this region within the years 1970–2000. The statistics does not allow to speculate on the size distribution of large fires. However, it is obvious that very small fires are most frequent. This particular effect is reproduced in our numerical results (Fig. 6 and 7).

## 5 Discussion

The bifurcation diagram for weights $`w^{}`$ in Figs. 33 reveals an unusual structure, which is a conglomerate of two diagrams. On the right side of the plot we observe some period doubling when the parameter $`r`$ decreases. On the left, clear period doublings are observed from a fixed point to a limit cycle of length two, and then to a cycle of length four, when the parameter $`r`$ increases. This left part of the bifurcation diagram is due to the size effect of finite lattice. On the other hand, the map for weights $`w^{\prime \prime }`$ constructed in the same way does not reveal wide windows of stability. This means that the period doublings and the windows are sensitive to the shape of the map on the right side of the percolation threshold. We could only prove that the cycles of all lengths are present there, by means of the presence of the cycle of period three. The map is similar to the tent map, where no windows are observed.

Concluding this point, the bifurcation diagram reveals that the map is not unimodal, what is not a surprise. Its structure may be caused by an intermittent character of the trajectories for small $`r`$. In this case, if the initial value of $`p`$ is small, forest fires do not occur in several years and the evolution is governed by its parabolic part $`p+rp(1p)`$, which remains almost constant through a long time. The diagrams in Figs. 3 and 3 are similar to the diagram obtained in . It was produced using a discontinuous map, where the time evolution was also intermittent. On the other hand, for larger $`r`$ the flat shape of the map above $`p_c`$ stabilizes the limit cycles. We note that the intermittent character of the problem is to some extent reflected also by the experimental data, where the majority of fires is of small sizes.

The histogram in Fig. 7 for weights $`w^{\prime \prime }`$ reveals an exponential curve of the number of fires of size $`A`$ in the range of intermediate and large $`A`$. The number of fires of size $`A`$ increases with $`A`$. This result is in contradiction to what is known from literature on the self-organized criticality . We interpret it as a particular consequence of our assumption that fires are initiated more likely in large clusters. In the approach presented in , a tree is ignited spontaneously if it is old enough.

Coming back to the problem of predictability, we note that when a deterministic mean field theory in the thermodynamic limit indicates the chaotic behavior of the system, it is hard to imagine that any other approach can produce a controllable solution. The exception is the the presence of the windows of stability, which allows us to predict the size of a forest fire for given values of $`r`$. Still, each limit cycle contains large fires, indicated by the letter $`R`$ in the respective word. Obviously, this kind of trajectory is not recommended by forest rangers. Concluding, both the bifurcation diagram and the Lyapunov exponent reveal that chaos is indeed present in the time evolution of the forest fire.

It seems also that the mean field approach can be a reference method in investigating other problems where probabilistic cellular automata are useful. A list of such problems includes immunology  and magnetic systems . To apply a map or a differential equation? The answer must be specific for a given problem.

## Acknowledgments

The authors thank to Professor Dietrich Stauffer for valuable comments. The simulations were carried out in ACK-CYFRONET-AGH. The time on SGI 2800 machine is financed by the Polish Committee for Scientific Research (KBN) with grant No. KBN/SGI2800/022/2002.