text:
# Classical Gluodynamics of High Energy Nuclear Collisions: an Erratum and an Update

## I Introduction

In a series of recent papers, we computed numerically the classical gluodynamics of the early stages of very high energy nuclear collisions. These computations were first performed for an SU(2) Yang-Mills gauge theory AR99 ; AR00 ; AR01  and later for the SU(3) Yang-Mills gauge theory AYR01 ; AYR02a ; AYR02b . In a very recent paper, Lappi has independently performed the same computation for the SU(3) gauge theory Lappi . We will show that the apparent discrepancies in the two computations can be understood as arising from an incorrect overall normalization in the non-perturbative formula for the energy in our previous work. Since the problem arises from a normalization error, all of our previous results are still useful if interpreted appropriately. The purpose of this note is primarily to clarify the sources of discrepancy in the two computations and to present corrected results wherever necessary. Now that two independent computations can be shown to agree, we will show that they significantly constrain phenomenological interpretations of the RHIC data.

This note is organized as follows. In Section II, we discuss the different conventions for the small x classical effective Lagrangean density used in the literature and isolate the source of the normalization error in our previous computations. A reader uninterested in the details can skip right ahead to Section III where we discuss the corrected results for various quantities. The implications of these results for RHIC phenomenology are discussed in Section IV.

## II Normalization of Classical Effective Lagrangean

There are two conventions for the classical effective Lagrangean density MV  that are followed in the literature. This is occasionally confusing so we shall discuss these explicitly here. The first of these is

$`L_I[A,\rho ]={\displaystyle \frac{1}{4}}F^aF^a\rho ^aA^a+i{\displaystyle \frac{\rho ^a\rho ^a}{2\mu _A^2}}.`$ (1)
For convenience, all Lorentz indices and space–time integrals have been suppressed here. For example, $`F^a`$ denotes the field strength tensor. In this form of the effective Lagrangean density, the coupling constant $`g`$ is absorbed in the definition of the charge density $`\rho ^a`$. The color charge squared per unit area $`\mu _A^2`$ contains a factor $`g^2`$. The correlator of color charges is then

$`\rho ^a(x_{})\rho ^a(y_{})=\mu _A^2\delta ^{ab}\delta ^{(2)}(x_{}y_{}).`$ (2)
This form of the effective Lagrangean density is used in several works – for instance, Refs. KovMuell ; IancuMcLerran .

The other form of the effective Lagrangean density found in the literature, for instance in Refs. MV ; JKW ; GyulassyMcLerran  is

$`L_{II}[A,\rho ]={\displaystyle \frac{1}{4}}F^aF^ag\stackrel{~}{\rho }^aA^a+i{\displaystyle \frac{\stackrel{~}{\rho }^a\stackrel{~}{\rho }^a}{2\mu _N^2}}.`$ (3)
Here $`g`$ is not absorbed in the definition of the charge density $`\rho ^a`$ but appears explicitly in the second term of the classical effective Lagrangean density. The relation between the color charge squared values in the two Lagrangean densities is therefore

$`\mu _A^2=g^2\mu _N^2.`$
The form of the effective Lagrangean density in Eq. (3) is the form used in our papers AR99 ; AR00 ; AR01 ; AYR01 ; AYR02a ; AYR02b , and the scale we use, $`\mathrm{\Lambda }_s^2`$ is defined to be $`\mathrm{\Lambda }_s^2=g^4\mu _N^2`$. Eq. (1) of Ref. Lappi  has the form

$`\rho ^a(x_{})\rho ^a(y_{})=g^2\mu _N^2\delta ^{ab}\delta ^{(2)}(x_{}y_{}).`$ (4)
Thus from Eqs. (2) and (II), it is clear that the form of the effective Lagrangean density used by Lappi is Eq. (1), albeit with the color charge squared defined to be identical to ours.

In Eq. (1) and (3), the field strength tensor has the canonical form

$`F_{\mu \nu }^a=_\mu A_\nu ^a_\nu A_\mu ^a+gf^{abc}A_\mu ^bA_\nu ^c.`$ (5)
In our numerical work – see for instance Ref. AR99  – the field strength tensor is defined as

$`\overline{F}_{\mu \nu }^a=_\mu \overline{A}_\nu ^a_\nu \overline{A}_\mu ^a+2f^{abc}\overline{A}_\mu ^b\overline{A}_\nu ^c.`$ (6)
Performing the field re-definition $`\overline{A}_\mu ^a=\frac{g}{2}A_\mu ^a`$, we obtain

$`\overline{F}_{\mu \nu }^a={\displaystyle \frac{g}{2}}F_{\mu \nu }^a.`$ (7)
The “canonical” effective Lagrangean density in Eq. (3) can then be re-written as

$`L_{II}`$ $`=`$ $`\left({\displaystyle \frac{2}{g}}\right)^2\left({\displaystyle \frac{1}{4}}\overline{F}^a\overline{F}^a+{\displaystyle \frac{g^2}{2}}\stackrel{~}{\rho }^a\overline{A}^a\right)`$ (8)
$`+`$ $`i{\displaystyle \frac{\stackrel{~}{\rho }^a\stackrel{~}{\rho }^a}{2\mu _N^2}}.`$
With a further re-definition $`\overline{\rho }^a=g^2\stackrel{~}{\rho }^a/2`$, we obtain

$`L_{(II)}`$ $`=`$ $`{\displaystyle \frac{4}{g^2}}\left({\displaystyle \frac{1}{4}}\overline{F}^a\overline{F}^a+\overline{\rho }^a\overline{A}^a\right)`$ (9)
$`+`$ $`i{\displaystyle \frac{\overline{\rho }^a\overline{\rho }^a}{2g^4\overline{\mu }_N^2}},`$
where $`\overline{\mu }_N=\mu _N/2`$. With the field re-definition, $`A^aA^a/g`$, both $`L_I`$ and $`L_{II}`$, as defined in Eqs. (1) and (3) respectively, can be expressed in the same form as the right hand side of Eq. (9), albeit without the factor of $`4`$ in front of the first two terms, and with $`\mu _N`$ instead of $`\overline{\mu }_N`$. Since the continuum limit of the lattice Hamiltonian in Ref. Lappi  corresponds to $`L_I`$ with $`A^aA^a/g`$ and $`\rho ^a=\overline{\rho }^a/g`$, there is this overall difference of a factor of 4 in the normalization between our lattice Hamiltonian and that of Lappi plus the requirement that the charge densities are related by $`\overline{\mu }_N=\mu _N/2`$.

Now consider the expression for the energy per unit rapidity per unit area,

$`\epsilon \tau ={\displaystyle \frac{1}{2}}\left(E^aE^a+B^aB^a\right),`$ (10)
again suppressing Lorentz indices. From Eq. (7),

$`\overline{B}^a={\displaystyle \frac{g}{2}}B^a;\overline{E}^a={\displaystyle \frac{g}{2}}E^a,`$ (11)
which gives

$`\epsilon \tau ={\displaystyle \frac{4}{g^2}}\overline{\epsilon }\tau .`$ (12)
In Refs. AR00  and AYR01 , we computed

$`\overline{\epsilon }\tau =\overline{f}_E\left(g^2\overline{\mu }_N\right)^3,`$ (13)
where $`\overline{f}_E`$ was determined by an extrapolation of the lattice results to the continuum limit. Then from Eq. (10) and (12), we obtain

$`\epsilon \tau ={\displaystyle \frac{4}{g^2}}\overline{f}_E\left(g^2\overline{\mu }_N\right)^3.`$ (14)
However, $`\overline{\mu }_N=\mu _N/2`$. Performing this substitution, we obtain

$`\epsilon \tau ={\displaystyle \frac{1}{g^2}}f_E\mathrm{\Lambda }_s^3,`$ (15)
where $`f_E=\overline{f}_E/2`$. Thus the physical value of the energy liberation coefficient is one half that computed previously. It accounts precisely for the result found by Lappi.

The number liberation coefficient remains the same after the field re-definitions. This is because the number per unit rapidity per unit area is defined to be

$`N={\displaystyle \frac{4}{g^2}}\overline{f}_N\left(g^2\overline{\mu }_N\right)^2.`$ (16)
Since $`\overline{\mu }_N=\mu _N/2`$, this replacement here exactly cancels the extra factor of $`4`$ appearing in the normalization. Hence $`\overline{f}_N=f_N`$ – our previous result is in agreement with the result of Lappi.

The gluon number distributions in the two papers coincide if we replace $`\mathrm{\Lambda }_s\mathrm{\Lambda }_s/2`$. The resulting comparison of the result of Lappi to our prior result in Ref. AYR01  is excellent and is shown in Fig. 1.

## III Corrected Results

Now that we understand the origin of the apparent discrepancies in the numerical results of Ref. Lappi  and our prior results, we will discuss the corrected results here. As discussed previously AR00 , all dimensional quantities in the classical effective theory can be expressed in terms of the appropriate power of $`\mathrm{\Lambda }_s`$ times a non-perturbative function of $`\mathrm{\Lambda }_sR`$. This is because $`\mathrm{\Lambda }_s`$ and $`R`$ are the only dimensional scales in the problem. $`\mathrm{\Lambda }_s`$ is, in principle, a function of the energy, the centrality and the atomic number. For RHIC energies, one can broadly estimate that $`\mathrm{\Lambda }_s1`$-2 GeV (more on this in the next section). The radius of a gold nucleus is approximately $`6.5`$ fm, so $`\mathrm{\Lambda }_sR`$ lies in the range $`32.5`$$`65`$ for Au-Au collisions at RHIC. We obtain, for the transverse energy,

$`{\displaystyle \frac{1}{\pi R^2}}{\displaystyle \frac{dE_T}{d\eta }}={\displaystyle \frac{f_E(\mathrm{\Lambda }_sR)}{g^2}}\mathrm{\Lambda }_s^3.`$ (17)
Previously, in Ref. AYR01 , we had $`f_E(25)=0.537`$ and $`f_E(83.7)=0.497`$. Following Eq. (15), this should be corrected to read $`f_E(50)=0.27`$ and $`f_E(167.4)=0.25`$. The corresponding expression for the gluon number is

$`{\displaystyle \frac{1}{\pi R^2}}{\displaystyle \frac{dN}{d\eta }}={\displaystyle \frac{f_N(\mathrm{\Lambda }_sR)}{g^2}}\mathrm{\Lambda }_s^2,`$ (18)
with $`f_N(50)=0.315`$ and $`f_N(167.4)=0.3`$. The variation in the number with $`\mathrm{\Lambda }_sR`$ is very small, on the order of 10%. These results are in good agreement with those of Lappi.

As we have explained in detail in our earlier papers, there are a variety of ways of defining the particle number for a generic field $`\varphi `$ and its conjugate momentum $`\pi `$. Let us compare the results for the following three definitions (here and in the following $`k_T`$ is the transverse momentum of the field):

$`n(\stackrel{}{k}_T)`$ $`=`$ $`\sqrt{|\varphi (\stackrel{}{k}_T)|^2|\pi (\stackrel{}{k}_T)|^2},`$ (19)
$`=`$ $`\omega (\stackrel{}{k}_T)|\varphi (\stackrel{}{k}_T)|^2`$ (20)
$`=`$ $`|\pi (\stackrel{}{k}_T)|^2/\omega (\stackrel{}{k}_T).`$ (21)
Here $`\omega (\stackrel{}{k}_T)`$ is the eigenfrequency for a free massless plane wave of the wave number $`\stackrel{}{k}_T`$ on a square lattice. In our work, as in Lappi’s, these prescriptions are applied to field configurations in the Coulomb gauge. In Fig. 1 our transverse momentum distributions for gluons are compared with Lappi’s result. KNV I (circles) in Fig. 1 is the same distribution as in our previous papers, but with the correct value of $`\mathrm{\Lambda }_s`$ used to set the horizontal scale. KNV I corresponds to the gluon number distribution from the definition Eq. (19),and $`k_T`$ dependence is extracted only for the transverse momenta along the principal lattice directions. The solid line is from Lappi’s numerical result, which is computed according to Eq. (21) using the entire Brillouin zone. KNV II is obtained by the same definition as Lappi’s and it is consistent. The deviation at large $`k_T`$ is considerd to be a consequence of how we define $`k_T`$ on the lattice. It should be noted that, with the lattice spacings used, the ultraviolet portion of the spectrum, $`k_T/\mathrm{\Lambda }_s3`$, cannot reliably reproduce continuum physics.

The initial transverse energy per particle changes (because $`f_E`$ does) and is now

$`{\displaystyle \frac{E_T}{N}}=0.9\mathrm{\Lambda }_s.`$ (22)
This value is nearly a factor of 2 lower than previously 111It is not exactly a factor of 2 because the value of $`\mathrm{\Lambda }_sR`$ at which $`f_E`$ is evaluated changes as well.. As we will discuss in the following section, this lower value of the transverse energy per particle makes the connection to RHIC phenomenology simpler to interpret.

The “formation time” $`\tau _i`$ (defined as $`\tau _i=1/\gamma (\mathrm{\Lambda }_sR)/\mathrm{\Lambda }_s`$) of gluons extracted in Ref. AR00  (it is the same for SU(2) and SU(3) as confirmed in Ref. AYR01 ) is unaffected, albeit the $`\mathrm{\Lambda }_sR`$ at which it is evaluated should be a factor of 2 larger than stated. One has approximately $`\gamma =0.3`$ in the range of interest. The initial energy density is then

$`\epsilon ={\displaystyle \frac{1}{\tau _i}}{\displaystyle \frac{dE_T}{\pi R^2d\eta }}={\displaystyle \frac{0.08}{g^2}}\mathrm{\Lambda }_s^4.`$ (23)
The number distributions, with coefficients appropriately scaled with $`\mathrm{\Lambda }_s`$, have the same form as previously:

$$\frac{1}{\pi R^2}\frac{dN}{d\eta d^2k_T}=\frac{1}{g^2}\overline{f}_n(k_T/\mathrm{\Lambda }_s),$$
(24)
where $`\overline{f}_n(k_T/\mathrm{\Lambda }_s)`$ is

$$\overline{f}_n=\{\begin{array}{cc}a_1\left[\mathrm{exp}\left(\sqrt{k_T^2+m^2}/T_{\mathrm{eff}}\right)1\right]^1\hfill & (k_T/\mathrm{\Lambda }_s1.5)\hfill \\ & \\ a_2\mathrm{\Lambda }_s^4\mathrm{log}(4\pi k_T/\mathrm{\Lambda }_s)k_T^4\hfill & (k_T/\mathrm{\Lambda }_s>1.5)\hfill \end{array}$$
(25)
with $`a_1=0.137`$, $`m=0.0358\mathrm{\Lambda }_s`$, $`T_{\mathrm{eff}}=0.465\mathrm{\Lambda }_s`$, and $`a_2=0.0087`$.

With this expression for the number distributions, we can now compute the occupation number for gluons after the collision. The relevant quantity for the validity of the three dimensional classical approximation is the three dimensional occupation number 222We thank Al Mueller for a discussion on this point. which can be estimated from the two dimensional boost-invariant number distribution computed in our simulation by the following relation MuellerQM2002 :

$`f^{\mathrm{glue}}={\displaystyle \frac{(2\pi )^3}{2(N_c^21)}}{\displaystyle \frac{1}{\pi R^2}}{\displaystyle \frac{dN}{d\eta d^2k_T}}.`$ (26)
Substituting our result for the number distribution from Eq. (19) here, one can compute the occupation number of gluons at the early stages of the collision. The results are shown in Fig. 2 for $`\mathrm{\Lambda }_sR=66`$.

Strictly speaking, the classical description is valid for $`f>1`$. It is not known however what the value of $`f`$ is at which the classical description breaks down completely and quantum corrections are dominant.

All of our other results can be understood by appropriate re-scaling. For instance, the $`p_t`$ distribution of the azimuthal asymmetry, $`v_2(p_t)`$, plotted versus $`p_t/\mathrm{\Lambda }_s`$ is the same after replacing $`\mathrm{\Lambda }_s\mathrm{\Lambda }_s/2`$ on the x-axis of the plot AYR02a . The peak in $`v_2(p_t)`$ is therefore at even lower momenta then previously – at $`p_t\mathrm{\Lambda }_s/8`$.

## IV Numerical Gluodynamics and RHIC Phenomenology

Now that we have shown that the results of two independent computations converge, it is useful to consider its implications for RHIC phenomenology. The results obtained in this classical approach are valid only at the very early stages of the nuclear collision and the connection to the final state distributions of measured hadrons might be considered remote. We will argue that, on the contrary, very general arguments can be applied to the RHIC data which will help constrain the properties of the initial state. In addition, they also significantly limit the freedom of final state interaction models which attempt to describe thermalization and the subsequent evolution of a Quark Gluon Plasma (QGP) in very high energy heavy ion collisions.

We know that the hadronic multiplicity at $`\eta 0`$ for central Au-Au collisions with center of mass energy $`\sqrt{s}=130`$ GeV/nucleon at RHIC is $`1000`$ while the hadronic transverse energy is $`500`$ GeV experiment . One expects, on general grounds, the following relations, respectively, to hold for the transverse energy and the multiplicity:

$`E_T^{\mathrm{glue}}`$ $`>`$ $`E_T^{\mathrm{hadrons}}`$ (27)
$`N^{\mathrm{glue}}`$ $``$ $`N^{\mathrm{hadrons}}.`$ (28)
The inequality for the transverse energy follows from the expectation that some of the initial transverse energy is converted into work due to the expansion of the system MuellerQM2002 . Moreover, if the system were to approach thermalization due to re-scattering, this process would again cause the initial transverse energy to be transferred into longitudinal energy Mueller2000 . Even if there were no $`pdV`$ work done and no re-scattering, it is difficult to envision a scenario where the final transverse energy is larger than the initial one. The second constraint is also very plausible. For a thermal system, entropy is simply proportional to the multiplicity and the final entropy must be equal to or greater than the initial one. Even in the other extreme scenario, independent fragmentation with no re-scattering, parton-hadron duality suggests that the parton multiplicity prior to hadronization is equal to the multiplicity of hadrons after phd .

Combining the first condition in Eq. (28) on the transverse energy with Eq. (17) and the empirical result for the RHIC data, we find that

$`\mathrm{\Lambda }_s^{\mathrm{RHIC}}>1.3\mathrm{GeV}.`$ (29)
Here we have assumed that the average transverse area for the range of impact parameters corresponding to the most central collisions is $`\pi R^2=130`$ fm2 and $`g=2`$ 333There is clearly an ambiguity in this bound due to our choice of the strong coupling constant $`g`$, which cannot be resolved at this order of the computation. If we assume that $`\alpha _s\alpha _s(\mathrm{\Lambda }_s)`$, the one loop value of $`\alpha _s`$ (from the bound we will set) is close to this estimate. Combining the second condition in Eq. (28) on the particle multiplicity (at $`\eta =0`$) with Eq. (18), we find

$`\mathrm{\Lambda }_s^{\mathrm{RHIC}}2\mathrm{GeV}.`$ (30)
Also, from Eq. (22), and Eqs. (29)-(30), we obtain the following constraint for the initial transverse energy per particle:

$`1.14\mathrm{GeV}{\displaystyle \frac{E_T}{N}}_{\mathrm{initial}}^{\mathrm{RHIC}}1.76\mathrm{GeV},`$ (31)
and, for the initial energy density, (from Eq. (23)),

$`7.1{\displaystyle \frac{\mathrm{GeV}}{\mathrm{fm}^3}}\epsilon _{\mathrm{initial}}^{\mathrm{RHIC}}40{\displaystyle \frac{\mathrm{GeV}}{\mathrm{fm}^3}}.`$ (32)
Thus the experimental data, the non-perturbative formulae in Eqs. (17)-(18) obtained from the classical numerical gluodunamics, and very general constraints can all be put together to significantly restrict the allowable range of $`\mathrm{\Lambda }_s^{\mathrm{RHIC}}`$. These in turn restrict the initial transverse energy per particle and the initial energy density. How these quantities evolve tell us a great deal about the space-time evolution of partonic matter in a heavy ion collision.

In our analysis, $`\mathrm{\Lambda }_s`$ is an external parameter – $`\mathrm{\Lambda }_s^2g^4\mu ^2`$, where $`\mu ^2`$, as discussed previously, has the physical interpretation of being the average color charge squared per unit area per unit rapidity of color sources at higher values of $`x`$ than those of interest. This scale is not quite the same as the gluon saturation scale $`Q_s^2`$, which denotes the scale at which the gluon distribution changes qualitatively due to saturation effects KovMuell . The relation between the two can be obtained by comparing the analytical expressions for the gluon distributions in the two approaches JKMW ; Kovchegov ; KovMuell . One obtains AYR02b

$`Q_s^2\mathrm{\Lambda }_s^2{\displaystyle \frac{N_c}{4\pi }}\mathrm{ln}\left({\displaystyle \frac{\mathrm{\Lambda }_s^2}{\mathrm{\Lambda }_{QCD}^2}}\right).`$ (33)
For the values of $`\mathrm{\Lambda }_s`$ in the range of interest, $`Q_s\mathrm{\Lambda }_s`$. Extrapolating the Golec-Biernat–Wusthoff parametrization to RHIC (with the appropriate $`A^{1/3}`$ dependence), we find $`\mathrm{\Lambda }_s^2Q_s^2=2`$ GeV2. This value for $`Q_s`$ is also obtained by Kharzeev and Nardi KN . $`\mathrm{\Lambda }_s`$ thus obtained lies within the range required by Eqs. (29) and Eq. (30). Nevertheless, there is an ambiguity in Eq. (33). For smaller values of $`x`$, one expects the infrared scale in this equation to be of order $`Q_s`$ and not $`\mathrm{\Lambda }_{QCD}`$ KIIM ; Mueller ; the relation between $`Q_s`$ and $`\mathrm{\Lambda }_s`$ will then be modified. Secondly, even at the classical level, the requirement that color neutrality be maintained Lam  will alter the relation between $`\mathrm{\Lambda }_s`$ and $`Q_s`$. Indeed, in this case, one can directly relate the two scales with an appropriate definition of $`Q_s`$. This point will be discussed further at a later date. Understanding this relation is important because the constraint on $`\mathrm{\Lambda }_s`$ discussed here will translate into a constraint on the gluon content in the nuclei at RHIC energies.

The bounds in Eq. (31) and Eq. (32) are very important for understanding the “late” stage dynamics of high energy heavy ion collisions, namely, when the classical approach breaks down as it must. One possibility is that the system undergoes hydrodynamic expansion Rasanen . At the end of the classical stage 444It is estimated in Ref. BMSS  that the average occupation number $`f`$, for $`\alpha _s1`$, is less than unity at $`\tau >1/Q_s/\alpha _{s}^{}{}_{}{}^{3/2}`$., the system is completely out of equilibrium. This is because the momentum distribution is very anisotropic: $`p_t\mathrm{\Lambda }_s`$ and $`p_z0`$. If the system is to thermalize, there must be particle production after the classical stage. Particle number conserving $`22`$ interactions are too inefficient to thermalize the system within a reasonable time scale Mueller2000 ; BMSS . If $`\mathrm{\Lambda }_s`$ is $`2`$ GeV, particles cannot be produced subsequent to the classical stage because this value saturates our bound for the particle number. If $`\mathrm{\Lambda }_s`$ is $`1.3`$ GeV (at the lower end of our allowed range), one can certainly increase the particle number by approximately a factor of 2 after the classical stage. The transverse energy per particle (see Eq. (31)) will then also decrease by a commensurate amount and approach the experimental value. However, one cannot then have any hydrodynamic expansion (and therefore thermalization) because for $`\mathrm{\Lambda }_s=1.3`$, the transverse energy was already as low as it could be: hydrodynamic expansion would lower it still further.

The problem with this “free streaming” scenario is that one does not have a satisfactory mechanism to explain the momentum dependence of the azimuthal anisotropy (especially at low $`p_t`$) observed at RHIC. The classical numerical simulations of $`v_2`$ do not generate enough $`v_2`$ AYR02a . There is an interesting attempt to explain the RHIC data as resulting from “non-flow” correlations but at low $`p_t`$ ($`p_t<Q_s`$) its energy dependence is different from the trend in the measured data from RHIC at $`\sqrt{s}=130`$ GeV/nucleon and $`\sqrt{s}=200`$ GeV/nucleon KovTuchin . A consistent picture of all the global features of the RHIC data might still be feasible within the straitjacket provided by the classical simulations and the data but there is limited room for models of final state interactions to maneuver.

###### Acknowledgements.

We would like to thank Rolf Baier, Keijo Kajantie, Dima Kharzeev, Tuomas Lappi, Larry McLerran and Al Mueller for very useful discussions and Tuomas Lappi and Keijo Kajantie for correspondence on their work. R. V.’s research was supported by DOE Contract No. DE-AC02-98CH10886 and the RIKEN-BNL Research Center. Y. N.’s research is supported by the DOE under Contracts DE-FG03-93ER40792. A. K. and R. V. acknowledge support from the Portuguese FCT, under grants POCTI/FNU/49529/2002 and CERN/FIS/43717/2001 and support in part of NSF Grant No. PHY99-07949. A. K. and R. V. also wish to thank the INT at the University of Washington for partial support during the completion of this work.


format: text

text:
# Neutrino-12C scattering in the ab initio shell model with a realistic three-body interaction

## Abstract

We investigate cross sections for neutrino-12C exclusive scattering and for muon capture on 12C using wave functions obtained in the ab initio no-core shell model. In our parameter-free calculations with basis spaces up to the $`6\mathrm{}\mathrm{\Omega }`$ we show that realistic nucleon-nucleon interactions, like e.g. the CD-Bonn, under predict the experimental cross sections by more than a factor of two. By including a realistic three-body interaction, Tucson-Melbourne TM(99), the cross sections are enhanced significantly and a much better agreement with experiment is achieved. At the same time, the TM(99) interaction improves the calculated level ordering in 12C. The comparison between the CD-Bonn and the three-body calculations provides strong confirmation for the need to include a realistic three-body interaction to account for the spin-orbit strength in $`p`$-shell nuclei.

The Gamow-Teller (GT) transition from the ground state of 12C to the 1+ T=1 isobar triplet (12Bg.s., 12C(15.11 MeV), 12Ng.s.) is a very sensitive test of nuclear structure models for mass 12 and, particularly, of the strength of the spin-orbit interaction. The two most common $`p`$-shell approximations for the structure of the ground state of 12C ((a) the $`p`$-shell equivalent of a L=0 S=0 three alpha-cluster structure and (b) the closed $`p_{3/2}`$ shell structure) give very different (indeed opposite) predictions for the B(GT) strength to T=1 1+ triplet. In the $`p`$-shell alpha-cluster limit the ground state of carbon has good SU(4) symmetry  and the Gamow-Teller transition is forbidden because there does not exist a 1+ T=1 state with  symmetry and $`\sigma \tau `$ operator cannot change SU(4) symmetry. This translates into an exact cancellation between the different $`p_{1/2}`$ and $`p_{3/2}`$ transition amplitudes. The observed transition strength requires the inclusion of higher SU(4) components in the wave functions and the breaking of the cancellation is quite sensitive to the assumed spin-orbit interaction. In the the $`jj`$-coupling limit, where one assumes that the ground state of 12C is described by a closed $`p_{3/2}`$ shell, the transition to the T=1 1+ state is pure $`p_{3/2}p_{1/2}`$. No cancellations between different transition amplitudes are allowed and the transition strength is over estimated by almost a factor of 6. When RPA correlations are included in the initial and final states the situation improves somewhat, but the transition remains over-estimated by about a factor of 4 kolbe ; nu\_C12\_SM . The strong contrast between the predictions of the pure $`jj`$-coupling and the pure SU(4) limits makes this Gamow-Teller transition an ideal test case for the strength of the spin-orbit interaction and for model wave functions of mass 12.

In this letter we present the predictions of no-core shell model (NCSM) C12\_NCSM  calculations of 12C for the T=1 1+ transition in 12C. We examine inelastic electron scattering to the 15.11 MeV state of 12C, muon capture to the ground state of 12B, and neutrino scattering to the ground state of 12N. These different electroweak reactions probe different momentum transfers and comparisons between theory and experiment allow us to test the convergence of the no-core shell model with increasing basis size, up to 6$`\mathrm{}\mathrm{\Omega }`$. We also investigate the contributions of a three-nucleon force since it is now well established pieper ; carlson ; wiringa  that realistic nucleon-nucleon interactions alone account for only half the observed $`p`$-shell splitting, while the rest arises form two-pion exchange between three or more nucleons. In the present calculations we include a realistic chiral-symmetry-based three-nucleon interaction (TNI), Tucson-Melbourne TM(99) TMprime99 .

A detailed description of the NCSM approach was presented, e.g. in Refs. C12\_NCSM ; v3eff . Here, we simply present extensions and modifications needed when a genuine TNI is included. The starting Hamiltonian is $`H_A=\frac{1}{A}_{i<j}\frac{(\stackrel{}{p}_i\stackrel{}{p}_j)^2}{2m}+_{i<j}^AV_{\mathrm{NN},ij}+_{i<j<k}^AV_{\mathrm{NNN},ijk}`$, where $`V_{\mathrm{NN},ij}`$ is the nucleon-nucleon (NN) interaction and $`V_{\mathrm{NNN},ijk}`$ is the TNI. We employ a large but finite harmonic-oscillator (HO) basis. Due to properties of the realistic nuclear interaction we have to derive an effective interaction appropriate for the selected finite basis space. To facilitate this, we modify the Hamiltonian by adding to it the center-of-mass (CM) HO Hamiltonian $`H_{\mathrm{CM}}=T_{\mathrm{CM}}+U_{\mathrm{CM}}`$, where $`U_{\mathrm{CM}}=\frac{1}{2}Am\mathrm{\Omega }^2\stackrel{}{R}^2`$, $`\stackrel{}{R}=\frac{1}{A}_{i=1}^A\stackrel{}{r}_i`$. The effect of the HO CM Hamiltonian will later be subtracted in the final many-body calculation. The modified Hamiltonian can be cast into the form $`H_A^\mathrm{\Omega }=H_A+H_{\mathrm{CM}}=_{i=1}^Ah_i+_{i<j}^AV_{ij}^{\mathrm{\Omega },A}+_{i<j<k}^AV_{\mathrm{NNN},ijk}`$,where $`h_i=\frac{\stackrel{}{p}_i^2}{2m}+\frac{1}{2}m\mathrm{\Omega }^2\stackrel{}{r}_i^2`$ and $`V_{ij}^{\mathrm{\Omega },A}=V_{\mathrm{NN},ij}\frac{m\mathrm{\Omega }^2}{2A}(\stackrel{}{r}_i\stackrel{}{r}_j)^2`$. Next we divide the $`A`$-nucleon infinite HO basis space into the finite active space ($`P`$) comprising all states up to $`N_{\mathrm{max}}`$ HO excitations above the unperturbed ground state and the excluded spaces ($`Q=1P`$). The basic idea of the NCSM approach is to apply a unitary transformation on the modified Hamiltonian, $`e^SH_A^\mathrm{\Omega }e^S`$ such that $`Qe^SH_A^\mathrm{\Omega }e^SP=0`$. If such a transformation is found, the effective Hamiltonian that exactly reproduces a subset of eigenstates of the full space Hamiltonian is given by $`H_{\mathrm{eff}}=Pe^SH_A^\mathrm{\Omega }e^SP`$. This effective Hamiltonian contains up to $`A`$-body terms and it is essentially as difficult to construct it as to solve the full problem. Therefore, we apply this approach with a cluster approximation. When a genuine TNI is considered, the simplest cluster approximation produces a three-body effective interaction. The NCSM calculation is then performed in four steps:

(i) We solve a three-nucleon system for all possible three-nucleon channels with the Hamiltonian $`H_A^\mathrm{\Omega }`$,i.e., using $`h_1+h_2+h_3+V_{12}^{\mathrm{\Omega },A}+V_{13}^{\mathrm{\Omega },A}+V_{23}^{\mathrm{\Omega },A}+V_{\mathrm{NNN},123}`$. Consequently, the three nucleons feel a pseudo-mean field of the spectator nucleons generated by the HO CM potential. It is necessary to separate the three-body effective interaction contributions from the TNI and from the two-nucleon interaction. Therefore, we need to find three-nucleon solutions for the Hamiltonian with and without the $`V_{\mathrm{NNN},123}`$ TNI term. The three-nucleon solutions are obtained by procedures described in Refs. Jacobi\_NCSM  (without TNI) and NCSM\_TM  (with TNI).

(ii) We construct the unitary transformation corresponding to the choice of the active basis space $`P`$ from the three-nucleon solutions using the Lee-Suzuki procedure LS1 ; UMOA . The three-body effective interaction is then obtained as $`V_{3\mathrm{e}\mathrm{f}\mathrm{f},123}^{\mathrm{NN}+\mathrm{NNN}}=P[e^{S_{\mathrm{NN}+\mathrm{NNN}}}(h_1+h_2+h_3+V_{12}^{\mathrm{\Omega },A}+V_{13}^{\mathrm{\Omega },A}+V_{23}^{\mathrm{\Omega },A}+V_{\mathrm{NNN},123})e^{S_{\mathrm{NN}+\mathrm{NNN}}}(h_1+h_2+h_3)]P`$ and $`V_{3\mathrm{e}\mathrm{f}\mathrm{f},123}^{\mathrm{NN}}=P[e^{S_{\mathrm{NN}}}(h_1+h_2+h_3+V_{12}^{\mathrm{\Omega },A}+V_{13}^{\mathrm{\Omega },A}+V_{23}^{\mathrm{\Omega },A})e^{S_{\mathrm{NN}}}(h_1+h_2+h_3)]P`$. The three-body effective interaction contribution from the TNI is then defined as $`V_{3\mathrm{e}\mathrm{f}\mathrm{f},123}^{\mathrm{NNN}}V_{3\mathrm{e}\mathrm{f}\mathrm{f},123}^{\mathrm{NN}+\mathrm{NNN}}V_{3\mathrm{e}\mathrm{f}\mathrm{f},123}^{\mathrm{NN}}`$.

(iii) As the three-body effective interactions are derived in the Jacobi-coordinate HO basis but the $`A=12`$ calculations will be performed in a Cartesian-coordinate single-particle Slater-determinant m-scheme basis, we need to perform a suitable transformation of the interactions. This transformation is a generalization of the well-known transformation on the two-body level that depends on HO Brody-Moshinsky brackets.

(iv) We solve the Schrödinger equation for the $`A=12`$ nucleon system using the Hamiltonian $`H_{A,\mathrm{eff}}^\mathrm{\Omega }=_{i=1}^Ah_i+\frac{1}{A2}_{i<j<k}^AV_{3\mathrm{e}\mathrm{f}\mathrm{f},ijk}^{\mathrm{NN}}+_{i<j<k}^AV_{3\mathrm{e}\mathrm{f}\mathrm{f},ijk}^{\mathrm{NNN}}`$, where the $`\frac{1}{A2}`$ factor takes care of overcounting the contribution from the two-nucleon interaction. At this point we also subtract the $`H_{\mathrm{CM}}`$. The $`A=12`$ nucleon calculation is then performed using the Many-Fermion Dynamics shell model code MFD  generalized to handle three-body interactions. Eventually, the transition densities are computed that serve as an input for evaluating our selected observables.

Detailed 12C NCSM calculations using realistic two-nucleon interactions were reported in Ref. C12\_NCSM . Here we extend those calculations by including the TNI and reach the $`4\mathrm{}\mathrm{\Omega }`$ ($`6\mathrm{}\mathrm{\Omega }`$) basis in calculations with (without) the TNI. In Table 1 we summarize some of our results. In general, in addition to increase of binding energy, we observe a substantial sensitivity of the low-lying spectra to the presence of the TNI and a trend toward level-ordering and level-spacing improvement in comparison to experiment. The sensitivity is the largest for states where the spin-orbit interaction strength is known to play a role. Note the correct ordering of the $`1^+04^+0`$ states and ordering and spacing improvement of the lowest T=1 states.

The significant increase in the spin-orbit splitting obtained from the inclusion of the TNI is seen most strikingly in the predicted cross sections to the T=1 1+ states. In all our electron scattering and weak interactions results here, only one-body currents are included and the bare operators are used. Dubach and Haxton dubach  have shown that at high-momentum-transfers it is necessary to include two-body meson-exchange currents to describe the transverse magnetic electron scattering form factor for excitation of the 15.11 MeV state. However, a reasonable description of the form factor up to momentum transfers of about 200 MeV/c can be obtained with a one-body current.

Figure 1 shows a comparison between the form factors predicted by the NCSM and experiment. The experimental data are represented by the black circles, which represent a fit to the data assuming only a one-body current obtained by Dubach and Haxton dubach . The theoretical curves shown are NCSM results for $`2\mathrm{}\mathrm{\Omega }`$, $`4\mathrm{}\mathrm{\Omega }`$, $`6\mathrm{}\mathrm{\Omega }`$ using the CD-Bonn NN interaction and a $`4\mathrm{}\mathrm{\Omega }`$ calculation using the AV8 plus the TM(99) realistic TNI. The qualitative features of our results are seen by looking at the height of the first maximum and the position of the minimum. With two-body interactions alone, the change in transition form factor from $`4\mathrm{}\mathrm{\Omega }`$ to $`6\mathrm{}\mathrm{\Omega }`$ is small compared to the differences between theory and experiment. The magnitude of the form factor is too low and minimum occurs too far out in momentum. When the TNI is included a significant improvement is seen in both the shape and magnitude of the theoretical form factor. The magnitude of the form factor up to the first maximum is close but somewhat lower than experiment. The shape of the form factor is also improved but it is still stretched out too far in momentum. Comparing the $`4\mathrm{}\mathrm{\Omega }`$ 3-body calculation with the $`4\mathrm{}\mathrm{\Omega }`$ CD-Bonn calculation the magnitude of the form factor at the peak has increased by about 75$`\%`$, and the position of the first minimum has shifted from $`q_{\mathrm{min}}`$ 400 MeV/c to 360 MeV/c. This improvement is almost entirely due to the improved strength of the spin-orbit splitting when the TNI is included.

The 15.11 MeV state was included in the fit to the $`p`$-shell interaction by Cohen and Kurath CK  (CK), and the CK interaction probably represents the best description of this transition using a globally fitted $`p`$-shell interaction. Our AV8+TM(99) form factor agrees well with the CK prediction, when for consistency between the two calculations, we use b=1.663 fm.

The conclusion drawn from the transverse magnetic form factor results is further supported by our B(M1; $`1^+10^+0`$) results presented in Table 1 and Fig. 2. The calculations with 2-body forces show saturation and underpredict the experiment by almost a factor of three. By including the TNI, the B(M1) value increases dramatically. We fully expect that further increases in the basis size will produce results with TNI close to experiment. For smaller basis sizes, effective transition operators may be important and work in this direction is underway.

Table 2 shows the comparison between the theoretical and experimental neutrino scattering cross section for the same selections of Hamiltonians and basis spaces. These results show a similar trend to the electron scattering results above. In this case, the neutrino spectrum for electron neutrinos from decay-at-rest (DAR) of the pion peaks around 30 MeV and the average momentum transfer is about 40 MeV/c. The CD-Bonn interaction (without TNI) results indicate an approach to convergence by $`6\mathrm{}\mathrm{\Omega }`$ but experiment is under-predicted by about a factor of 2.4. When the TNI is included with the AV8 interaction the predicted cross section is only 30$`\%`$ lower than experiment. Based on the similarity of trends with the electron scattering results, we anticipate that when the model space is eventually expanded to $`6\mathrm{}\mathrm{\Omega }`$ theory would be within 15$`\%`$ of experiment. The substitution of AV8 for CD-Bonn in the calculations with TNI is expected to be of minor consequence.

Muon capture involves a higher momentum transfer than the $`(\nu _e,e^{})`$ reaction and the average momentum transfer is $`q`$ 100 MeV/c. By $`6\mathrm{}\mathrm{\Omega }`$ the CD-Bonn calculations show signs of converging yet experiment is underestimated a factor of 2.6. The inclusion of the TNI shows a significant improvement and, for $`4\mathrm{}\mathrm{\Omega }`$, theory is 34$`\%`$ lower than experiment. Again extrapolating using the trends of the inelastic electron scattering results suggests that a 6$`\mathrm{}\mathrm{\Omega }`$ calculation that included a realistic TNI would come within 20$`\%`$ of experiment.

The $`(\nu _\mu ,\mu ^{})`$ neutrino cross section to 12Ng.s. corresponds to the LSND muon neutrinos from decay-in-flight (DIF) of the pion. This spectrum involves neutrinos up to about 250 MeV, with a average neutrino energy of about 150 MeV and an average momentum transfer of about 200 MeV/c. In this case the $`6\mathrm{}\mathrm{\Omega }`$ CD-Bonn calculation is off by a factor of 1.8 compared with experiment. The $`4\mathrm{}\mathrm{\Omega }`$ calculation that includes the 3-body TM(99) interaction is, in fact, in agreement with experiment. However, based on the trends established above, this suggests that a larger model space may over-predict experiment. Examining the elastic scattering form factor suggests that the problem lies in the fact that at 200 MeV/c the predicted form factor is too large. Of course, as the model space is increased we expect the form factor to be shifted down in momentum.

In conclusion, the transition from the 12Cg.s. to the T=1 1+ states in mass 12 is very sensitive to the strength of the spin orbit interaction. We have investigated neutrino-12C exclusive cross sections and muon capture on 12C as well as inelastic electron scattering using wave functions obtained in the ab initio NCSM. In our parameter-free calculations with basis spaces up to $`6\mathrm{}\mathrm{\Omega }`$ we show that realistic NN interactions under predict the experimental weak interaction cross sections by more than a factor of two. At high momentum transfers around $`q`$ 200 MeV/c the electron scattering form factor is over-predicted and the position of the predicted minimum to close to 400 MeV/c, compared to the experimental minimum at $`q260`$ MeV/c. By including a realistic TNI the weak interaction cross sections are enhanced significantly which considerably improves agreement with experiment. The shape of the electron scattering form factor is also significantly improved, but the predicted form factor still peaks at too large a momentum transfer and is too large at the momentum transfers relevant to the LSND DIF cross section. The difference between the observed and predicted shape for the (e,e’) form factors and the very different momentum transfers involved in the three weak processes examined here imply that a single experiment/theory scale factor cannot be defined for all three. The comparison between the CD-Bonn and the three-body calculations discussed here provide a strong confirmation of the need to include a realistic three-body interaction to account for the spin-orbit strength in $`p`$-shell nuclei.

This work was performed in part under the auspices of the U. S. Department of Energy by the University of California, Lawrence Livermore National Laboratory under contract No. W-7405-Eng-48. This work was also supported in part by USDOE DE-FG02-87ER40371.


text:
# Chapter 1 A Fast Hybrid Approach to Air Shower Simulations and Applications

## Chapter 1 A Fast Hybrid Approach to Air Shower Simulations and Applications

### Abstract

The SENECA model, a new hybrid approach to air shower simulations, is presented. It combines the use of efficient cascade equations in the energy range where a shower can be treated as one-dimensional, with a traditional Monte Carlo method which traces individual particles. This allows one to reproduce natural fluctuations of individual showers as well as the lateral spread of low energy particles. The model is quite efficient in computation time.

As an application of the new approach, the influence of the low energy hadronic models on shower properties for AUGER energies is studied. We conclude that these models have a significant impact on the tails of lateral distribution functions, and deserve therefore more attention.

### 1. Introduction

In the field of air shower simulations, large computation times limit the applicability of traditional Monte-Carlo approaches at very high energies. The thinning algorithm tries solve this problem, at the price of introducing artificial fluctuations. Here, we present a hybrid approach to air shower simulation, as introduced in : A system of transport equations is employed in the energetic region where an air shower can be considered as one dimensional. The lateral spread of low energy particles is calculated with traditional Monte-Carlo methods, where one traces each particle individually.

### 2. Hadronic Cascade Equations

If one considers an air shower as a one-dimensional system, it can be described with the following set of cascade equations (CE) :

$`{\displaystyle \frac{h_n(E,X)}{X}}`$ $`=`$ $`h_n(E,X)\left[{\displaystyle \frac{1}{\lambda _n(E)}}+{\displaystyle \frac{B_n}{EX}}\right]`$
$`+{\displaystyle \underset{m}{}}{\displaystyle _E^{E_{\mathrm{max}}^{\mathrm{had}}}}h_m(E^{},X)\left[{\displaystyle \frac{W_{mn}(E^{},E)}{\lambda _m(E^{})}}+{\displaystyle \frac{B_mD_{mn}(E^{},E)}{E^{}X}}\right]𝑑E^{},`$
where $`h_n(E,X)dE`$ is the number of particles of type $`n`$ (basically nucleons, pions and kaons) with an energy $`[E,E+dE]`$ at an atmospheric depth $`X`$. The first line with the minus-sign accounts for particles falling out of the system by either inelastic collisions with air-nuclei or by decay. The second line accounts for producing particles by the same two mechanisms at higher energies, $`W_{mn}(E,E^{})=\frac{dN}{dE^{}}`$ is the differential energy spectrum for collisions of particle type $`m`$ at energy $`E`$ and $`D_{mn}`$ are the corresponding decay functions. The technique to solve these equations is described in detail in references .

### 3. Electromagnetic Cascading

Electromagnetic cascades can in principle be described in a similar way. However, due to the fact that no decays are involved it is sufficient to compute the energy spectra of shower particles for a reference layer of air of a given thickness, e.g. $`2.5\mathrm{g}/\mathrm{cm}^2`$. We use discretized energies as $`E_i=10^{(i1)/n_d}/rmGeV`$, with $`n_d`$ bins per decade. If $`V_{ji}^{mn}`$ describes the energy of electrons/positrons and photons, for a particle of type $`m`$ traversing this layer of air, then a simple iteration

$$g_i^n(X+\mathrm{\Delta }X)=\underset{m,ji}{}g_j^m(X)V_{ji}^{mn}(\mathrm{\Delta }X).$$
(2)
describes the shower at any depth $`X`$.

### 4. Source Functions

The system of hadronic and electromagnetic cascade equations is solved down to energies of $`1000`$ GeV and $`10`$ GeV, respectively. Below these energies, the three-dimensional character of the air shower becomes important. Thus, instead of evolving further, we calculate directly the energy spectra of low energy particles in the following way:

$`{\displaystyle \frac{h_n^{\mathrm{source}}(E,X)}{X}}`$ $`=`$ $`{\displaystyle \underset{m}{}}{\displaystyle _{E_{\mathrm{min}}^{\mathrm{had}}}^{E_{\mathrm{max}}^{\mathrm{had}}}}h_m(E^{},X)\left[{\displaystyle \frac{W_{mn}(E^{},E)}{\lambda _m(E^{})}}+{\displaystyle \frac{B_mD_{mn}(E^{},E)}{E^{}X}}\right]dE^{}.`$
The source function is then used to generate low energy particles which are then followed in a standard Monte-Carlo approach.

### 5. Tests of Results

The use of cascade equations should be of technical nature only and not change the physics results. Therefore we compare results to the pure Monte Carlo approach with thinning. The average of 500 proton induced $`10^{19}`$ eV showers with thinning level $`10^5`$ has been computed. Fig. 1 shows comparisons for longitudinal and lateral profiles as well as the energy and arrival time spectra at 1000 m distance from the shower axis. The agreement between the CE and the MC approach is satisfactory.

### 6. Low Energy Model Dependence on Lateral Distribution Functions

In Ref.  it was pointed out that the low-energy hadronic model is of crucial importance for muons in general, and for the tails the lateral distribution functions (LDF) of electrons and photons. The reason was that electromagnetic subshowers contributing to large distances are emitted by low energy particles at a large distance from the shower axis.

Here we want to analyse the influence of different low energy hadronic models on the LDF of electrons, muons and photons. We use the following models: GHEISHA  and G-FLUKA  as found in GEANT3.21, and the UrQMD  model developed at the University of Frankfurt. For high energies we employ the QGSJET01  model. In Fig. 2 we show the ratio of LDFs of $`e^\pm ,\mu ^\pm `$ and $`\gamma `$ for the different models. We see how the tails of the distributions change due to different choices of the low-energy hadronic model. GHEISHA produces more muons than UrQMD, since it gives a larger yield for charged pion production as seen in Fig. 2. The differences even influence the tails of the LDF of electrons and photons, though less pronounced. The results for UrQMD and G-FLUKA are more similar to each other.

### 7. Acknowledgments

This works has been supported by NASA-Grant NAG-9246 and by the German BMBF/Desy. The computational resources were provided by the Center for Scientific Computing Frankfurt (CSC).

### 8. References

1. Bossard et al.,Phys.Rev.D63:054030,2001.

2. H.J. Drescher, G.R. Farrar, \[astro-ph/0212018\], in print Phys.Rev.D

3. H.J. Drescher, G.R. Farrar, Astropart.Phys.19:235-244,2003

4. N.N. Kalmykov and S.S. Ostapchenko and A.I. Pavlov ,Nucl. Phys. B 17,1997

5. S.A.Bass et al. Prog.Part.Nucl.Phys. 41 (1998) 225

6. M.Bleicher et al. J.Phys. G25 (1999) 1859

7. H.Fesefeldt PITHA 85/02, Aachen 1985

8. A.Fasso et al. Proceedings of the Workshop on Simulating Accelerator Radiation Environments, Santa Fe, 1993


text:
# 1 Introduction

## 1 Introduction

The investigation in phase space of the properties of quantum systems contributes with a remarkable picture of those systems in which the states are described by areas in such a way that the emergence of diverse quantum phenomena  appears as the result of the interference of the portions of these areas which overlap. The area which represent each state may be fixed from considerations in the context of the old Bohr-Sommerfeld approach , with the aid of phase-space distributions  or using the WKB approximate wave functions . This last approach has the advantage of being more systematic but is constrained by the limitations of the approximation. In this paper we illustrate both terms in this sentence discussing the photon distribution of squeezed number states. The photon distribution of squeezed number states is parity sensible and present characteristic oscillations with increasing photon number until it decay and vanish for sufficiently large photon number. These oscillations which may be computed exactly  are similar to the ones which appear for squeezed states which been discussed thoroughly using the phase space description . It is then interesting to investigate also how this approach may also be applied to the case of the squeezed number states. A first step in this direction was taken in Ref.  where it was shown that the main characteristics of the distribution are consistent with interference in phase space picture. Here we develop with detail the WKB approach of the phase space picture to investigate how the parity restrictions may be enforced and under which conditions the approximation breaks down. The outcome is that the WKB approximation is not able to represent the last maximum of the distribution in the the high photon number region. This failure may be traced to to the limitations of the WKB approximation and not to the phase space picture which as we show below is still realized (at least for high squeezing) using the more direct scheme to assign the areas to the states which relies on the Wigner-Cohen distributions.

## 2 The WKB approximation

Let us begin discussing the WKB approximation for squeezed number states and its relation with the phase space picture. We work with the normalized one mode Hamiltonian

$$H=\frac{1}{2}p^2+\frac{1}{2}x^2=a^{}a+\frac{1}{2}$$
(1)
and denote by $`|n>`$ the eigenstates of $`H`$ and of the number operator $`N=a^{}a`$ with eigenvalues $`n+1/2`$ and $`n`$.

The squeezed number states  are defined by

$`|n,r=S(r)|n`$ (2)
with the squeeze operator $`S(r)`$ given by

$$S(r)=\mathrm{exp}\left(\frac{1}{2}r(a^{})^2\frac{1}{2}ra^2\right).$$
(3)
(We are taking the phase parameter of the squeezing equal to zero since the photon distribution does not depends on it). These states are eigenstates with eigenvalues $`(n+1/2)`$ of the transformed Hamiltonian $`H^{}=S(r)HS(r)^1`$ which is shown to be equal to ,

$$H^{}=\frac{1}{2}p^2e^{2r}+\frac{1}{2}x^2e^{2r}.$$
(4)
Alternatively the squeezed number states may be interpreted as the Fock states of a related system with mass $`e^{2r}`$ and the same frequency. For a more detailed description of these states see . The WKB wave functions for these states are given by ,

$$\mathrm{\Phi }_n^{(r)}(x)=\frac{2}{C_n}(T_n^{(r)}p_n^{(r)}(x))^{1/2}\mathrm{cos}\left(S_n^{(r)}(x)\pi /4\right),$$
(5)
where $`C_n`$ is a normalization constant. The value of $`p_n^{(r)}`$ is obtained from the equation

$$n+1/2=\frac{1}{2}(p_n^{(r)}(x))^2e^{2r}+\frac{1}{2}x^2e^{2r}$$
(6)
and is given in terms of the classical turning point value $`\epsilon _n^{(r)}=e^r\sqrt{2m+1}`$ by

$`p_n^{(r)}(x)=`$ $`e^{2r}((\epsilon _n^{(r)})^2x^2)^{1/2}`$
$`=`$ $`e^rp_n^{(0)}(xe^r)`$ (7)
The phase function is given by,

$`S_n^{(r)}(x)=`$ $`{\displaystyle _x^{\epsilon _n^{(r)}}}p_n^{(r)}(x^{})𝑑x^{}`$
$`=`$ $`S_n^{(0)}(xe^r),`$ (8)
and the normalization constant has been written in terms of the period of the related classical motion,

$`T_n^{(r)}=`$ $`{\displaystyle _{\epsilon _n^{(r)}}^{\epsilon _n^{(r)}}}{\displaystyle \frac{dx^{}}{p_n^{(r)}(x^{})}}=`$ $`2\pi e^{2r}`$ (9)
For $`r=0`$ we recover the description of the number states. The WKB wave functions are compressed by the squeezing in the same form as the exact states,

$$\mathrm{\Phi }_n^{(r)}(x)=e^{r/2}\mathrm{\Phi }_n^{(0)}(e^rx).$$
(10)
We note that $`p_n^{(r)}`$ is a symmetric function and that

$$_{\epsilon _n^{(r)}}^{\epsilon _n^{(r)}}p_n^{(r)}(x^{})𝑑x^{}=\pi /2+n\pi .$$
(11)
Hence the phase function $`S_n(x)`$ satisfies the relation

$$S_n^{(r)}(x)\pi /4=\left(S_n(x)\pi /4\right)+n\pi $$
(12)
and one can show that the WKB wave functions satisfy the parity conditions of the exact states.

$`\mathrm{\Phi }_n^{(r)}(x)`$ $`=`$ $`(1)^n\mathrm{\Phi }_n^{(r)}(x)`$ (13)
Again for $`r=0`$ we are dealing with the number states.

The WKB wave functions are meaningless near the classical return points where they diverge. Also we note that the approximation holds better for higher values of $`n`$. The normalization of the wave functions is such that

$`|C_n|^2=1+{\displaystyle \frac{1}{\pi }}{\displaystyle _{\epsilon _n^{(r)}}^{\epsilon _n^{(r)}}}{\displaystyle \frac{\mathrm{sin}(2S_m^{(r)}(x))}{p_m^{(r)}(x)}}𝑑x`$ (14)
with the integral in the right term vanishing for large n.

## 3 Photon statistics in the WKB approximation

Now we turn to the computation of the photon distribution of the squeezed number states within this approximation. It is given by

$$P_{mn}=|W_{mn}|^2$$
(15)
with

$$W_{mn}=<m|n,r>=_{\mathrm{}}^{\mathrm{}}\mathrm{\Phi }_m^{(0)}(x)\mathrm{\Phi }_n^{(r)}(x)𝑑x.$$
(16)
Using the WKB wave functions is shown that, as was suggested in References  the probability amplitude $`W_{mn}`$ may be expressed in terms of a phase $`\phi _{mn}`$ and an overlapping area $`A_{mn}`$ in the form,

$`W_{mn}`$ $`=`$ $`(A_{mn})^{1/2}e^{i\phi _{mn}}+(A_{mn})^{1/2}e^{i\phi _{mn}}`$ (17)
$`+`$ $`(A_{mn})^{1/2}e^{i\phi _{mn}^{}}+(A_{mn})^{1/2}e^{i\phi _{mn}^{}}`$
In this expression $`A_{mn}`$ is given by

$$A_{mn}=2\pi \left(T_mT_n^{(r)}p_m^2(X_c)\right)^1|\frac{d}{dx}p_m(X_c)\frac{d}{dx}p_n^{(r)}(X_c)|^1$$
(18)
and is shown to be the overlapping area between a circular ring of interior radius $`\sqrt{2m}`$ and exterior radius $`\sqrt{2m+2}`$ of total area $`2\pi `$ representing the $`m`$ Fock state in phase space and a deformed elliptical ring of equal area which corresponds to the $`n`$ squeezed number state. The phases are given by

$`\phi _{mn}=`$ $`S_m^{(0)}(X_c)S_n^{(r)}(X_c)\pi /4,`$
$`\phi _{mn}^{}=`$ $`S_m^{(0)}(X_c)S_n^{(r)}(X_c)+\pi /4`$ (19)
where $`X_c`$ is the point where,

$$p_m^{(0)}(X_c)=p_n^{(r)}(X_c)$$
(20)
as can be shown in Figure (1). We have then,

$$X_c=\left(\frac{e^{2r}(2n+1)(2m+1)}{(e^{4r}1)}\right)^{1/2}.$$
(21)
From the properties of $`S_m^{(r)}`$ (Ec. 12) we obtain:

$`\phi _{mn}^{}`$ $`=`$ $`\phi _{mn}+(mn)\pi `$ (22)
Using (22) and (17) the probability amplitude is given by

$$W_{mn}=2(A_{mn})^{1/2}\mathrm{cos}(\phi _{mn})\left(1+\mathrm{cos}\left((mn)\pi \right)\right).$$
(23)
This result is in agreement with the form of the amplitude proposed in  and in particular 22 gives the correct parity behavior of the distribution.

The connection with the geometrical construction is obtained noting that $`S_m^{(0)}(X_c)`$ is given by the shadowed area in Fig.(2). Then,

$$S_m^{(0)}(X_c)=\frac{1}{2}\theta (\epsilon _m^{(0)})^2\frac{1}{2}X_c\left((\epsilon _m^{(0)})^2X_c^2\right)^{1/2}$$
(24)
where

$$\theta =\mathrm{arccos}\{\frac{X_c}{\epsilon _m^{(0)}}\}$$
(25)
In the same way we have,

$$S_n^{(r)}(X_c)=\frac{1}{2}\theta (\epsilon _n^{(r)})^2\frac{1}{2}X_c\left((\epsilon _n^{(r)})^2X_c^2\right)^{1/2}$$
(26)
where

$$\theta =\mathrm{arccos}\{\frac{X_c}{\epsilon _n^{(r)}}\}$$
(27)
In Fig (3) we show the WKB approximation for the photon distribution of a squeezed number state with $`n=5`$ and $`r=2`$ compared with the exact computation. The approximation is faithful for photon numbers below the last maximum.

## 4 Wigner-Cohen distributions and interference in phase space

The failure of the approximation presented in the last section for large photon numbers is related to the the breakdown of the validity of the WKB approximation. Nevertheless as we show below it does not exclude the description of this part of the distribution using the area overlapping in phase space. To see this let us introduce $`F_n^{(r)}(x,p)`$ to denote any of the phase space Wigner-Cohen distributions  associated to the squeezed number states which satisfies,

$$F_n^{(r)}(x,p)𝑑p=|\psi _n^{(r)}(p)|^2$$
(28)
and

$$F_n^{(r)}(x,p)𝑑x=|\overline{\psi }_n^{(r)}(p)|^2$$
(29)
where $`\psi _n^{(r)}(x)`$ y $`\overline{\psi }_n^{(r)}(p)`$ are the wave functions of the squeezed number states in the position and momentum representations respectively.

The photon number distribution may be obtained by integrating $`F_n^{(r)}(x,p)`$ over the $`m`$-th ring $`\mathrm{\Omega }_m`$ associated to the number state $`|m>`$ which is the area between to circles of radii $`\sqrt{2m}`$ and $`\sqrt{2m+2}`$ respectively,

$$A_{mn}=\frac{1}{2}_{\mathrm{\Omega }_m}𝑑p𝑑xF_n^{(r)}(x,p)$$
(30)
For a highly squeezed number state the integral in the configuration space variable $`q`$ can approximately taken over the whole line

$$A_{mn}_{\sqrt{2m}}^{\sqrt{2m+2}}|\overline{\psi }_n^{(r)}(p)|^2𝑑p$$
(31)
and by supposing a smooth behavior in the ring we get,

$$A_{mn}|\overline{\psi }_n^{(r)}(\sqrt{2m+1})|^2(\sqrt{2m+2}\sqrt{2m})$$
(32)
Taking into account the parity of the functions we finally have,

$$P_{mn}\left[1+\mathrm{cos}((mn)\pi )\right]^2A_{mn}$$
(33)
Substituting the values of the waves functions we obtain,

$`P_{mn}`$ $``$ $`\left(1+\mathrm{cos}((mn)\pi )\right)^2\left(\sqrt{2m+2}\sqrt{2m}\right)`$ (34)
$`e^r{\displaystyle \frac{e^{(2m+1)e^{2r}}(H_n(\sqrt{2m+1}e^r))^2}{\sqrt{\pi }2^nn!}}`$
which is in an excellent agreement with the exact value as shown in Fig. (4) in the particular case of $`n=5`$ and $`r=2`$. The agreement stills holds for photon numbers beyond the last oscillation where the WKB approximation breaks down. This expression also shows that for large squeezing the photon number oscillation are driven by the oscillations of the wave function in momentum space.

## 5 Conclusion

We have computed the oscillations photon number distribution of squeezed number states within the WKB approximation making connection with the phase space description given in Ref.. In this case there are four sectors where the areas overlap and give a contribution to the probability amplitude . We show that for high photon number the WKB approximation fails to describe the last oscillation. Finally we note that a description in terms of interference of the overlapping areas in phase space may be recovered, but not based in the WKB approximation but in the Wigner-Cohen distributions. For large squeezing this prescription reproduces the photon number distribution for the whole range.


text:
# Phase separation in coupled chaotic maps on fractal networks

## Abstract

The phase ordering dynamics of coupled chaotic maps on fractal networks are investigated. The statistical properties of the systems are characterized by means of the persistence probability of equivalent spin variables that define the phases. The persistence saturates and phase domains freeze for all values of the coupling parameter as a consequence of the fractal structure of the networks, in contrast to the phase transition behavior previously observed in regular Euclidean lattices. Several discontinuities and other features found in the saturation persistence curve as a function of the coupling are explained in terms of changes of stability of local phase configurations on the fractals.

Coupled map lattices have provided fruitful and computationally efficient models for the study of a variety of dynamical processes in spatially distributed systems Chaos . The discrete-space character of coupled map systems makes them specially appropriate for the investigation of spatiotemporal dynamics on nonuniform or complex networks. Phenomena such as pattern formation, spatiotemporal intermittency, nontrivial collective behavior, synchronization, etc., have been extensively studied in coupled map systems defined on fractal lattices CK2 , hierarchical structures Gade , trees Tree , random graphs Vol , small-world networks Kay , and scale-free networks Amri .

Recently, there has been much interest in the study of the phase-ordering properties of systems of coupled chaotic maps and their relationship with Ising models in statistical physics Chate ; Chate2 ; Wei ; Stra ; Just . These works have invariably assumed the phase competition dynamics taking place on a uniform Euclidean space; however, in many physical situations the medium that supports the dynamics can be nonuniform on some length scales. The nonuniformity may be due to the intrinsic heterogeneous nature of the substratum such as porous or fractured media, or it may arise from random fluctuations in the medium. This article investigates the phenomenon of phase ordering in coupled chaotic maps on fractal networks as a model for studying this phenomenon on nonuniform media. The class of fractal networks being considered corresponds to generalized Sierpinski gaskets (GSG) embedded in Euclidean spaces of arbitrary dimension $`d`$ CK2 . In particular, this model of coupled maps on fractal networks yields a situation to explore the role that the connectivity of the underlying lattice plays on the statistical properties of phase ordering processes in nonlinear coupled systems.

Deterministic fractal networks, such as GSG, can be generated in any $`d`$-dimensional Euclidean space as follows CK2 . At the $`n`$th level of construction, the fractal consists of $`N=(d+1)^n`$ $`d`$-dimensional hypertetrahedral cells whose coordinates can be specified by a sequence $`(\alpha _1\alpha _2\mathrm{}\alpha _n)`$, where $`\alpha _m`$ can take any value in a set of $`(d+1)`$ different symbols which can be chosen to be $`\{0,1,\mathrm{},d\}`$. At level $`n+1`$, each cell $`(\alpha _1\alpha _2\mathrm{}\alpha _n)`$ splits into $`d+1`$ cells scaled down by a longitudinal factor of two, and which are now labeled by $`(\alpha _1\alpha _2\mathrm{}\alpha _n\alpha _{n+1})`$, where the first $`n`$ symbols of the sequence are the same as the parent cell. Given this construction rule, the fractal dimension of the GSG is $`d_f=\mathrm{log}(d+1)/\mathrm{log}2`$. A label $`(\alpha _1\alpha _2\mathrm{}\alpha _n)`$ can be written as $`(\alpha _1\alpha _2\mathrm{}\alpha _{ns}\alpha _{ns+1}^s)`$ for some $`s\{1,2,\mathrm{},n\}`$, where $`\alpha _i^s`$ means the sequence of $`s`$ identical symbols $`\alpha _i`$. The cell with this label has a neighborhood $`𝒩_{(\alpha _1\alpha _2\mathrm{}\alpha _n)}`$ with $`d+1`$ elements labeled by $`(\alpha _1\mathrm{}\alpha _{n1}(\alpha _n+j))`$ $`(j=1,2,\mathrm{},d)`$ and $`(\alpha _1\mathrm{}\alpha _{ns+1}\alpha _{ns}^s)`$, where the addition $`\alpha _i+j`$ is defined modulo $`(d+1)`$. If $`s=n`$, then the cell is one of the $`d+1`$ vertices of the gasket, labeled by $`(\alpha _1^n)`$, and it has only $`d`$ neighbors belonging to the same parent cell. An integer index $`i=0,1,\mathrm{},(d+1)^n1`$, can be assigned to each cell of the lattice at the level of construction $`n`$ by the rule $`i=_{m=1}^n\alpha _m(d+1)^{nm}`$.

The equations describing the dynamics of the diffusively coupled map system defined on these fractal networks at level of construction $`n`$, embedded in a $`d`$-dimensional Euclidean space, are

$$\begin{array}{cc}x_{t+1}(\alpha _1\mathrm{}\alpha _n)=(1ϵ)f\left(x_t(\alpha _1\mathrm{}\alpha _n)\right)+\hfill & \\ & \\ \frac{ϵ}{d+1}_{(\beta _1\mathrm{}\beta _n)𝒩_{(\alpha _1\mathrm{}\alpha _n)}}f\left(x_t(\beta _1\mathrm{}\beta _n)\right),\hfill & \end{array}$$
(1)
where $`x_t(\alpha _1\mathrm{}\alpha _n)`$ gives the state of the cell $`(\alpha _1\mathrm{}\alpha _n)`$ at discrete time $`t`$; $`(\alpha _1\mathrm{}\alpha _n)`$ and $`(\beta _1\mathrm{}\beta _n)`$ label the $`(d+1)^n`$ cells on the gasket; $`ϵ`$ is a parameter measuring the coupling strength between neighboring sites, and $`f(x)`$ is a nonlinear function that expresses the local dynamics. Equation (1) also applies to the $`(d+1)`$ vertex cells of the fractal network, except that the coefficient of the sum is $`ϵ/d`$, since each of these cells have $`d`$ neighbors.

As local dynamics, we assume a piecewise linear, odd map Chate

$$f(x)=\{\begin{array}{cc}2\mu /3\mu x\hfill & \text{if }x[1,1/3]\hfill \\ \mu x\hfill & \text{if }x[1/3,1/3]\hfill \\ 2\mu /3\mu x\hfill & \text{if }x[1/3,1].\hfill \end{array}$$
(2)
When the parameter $`\mu [1,2]`$, the map possesses two symmetric chaotic attractors contained in the intervals $`I^\pm =[\pm \mu (2\mu )/3,\pm \mu /3]`$, and separated by a gap $`I^o=[\mu (2\mu )/3,\mu (2\mu )/3]`$. For values of $`\mu `$ close to two, the size of the chaotic intervals is larger than the gap. Then the local states have two well defined symmetric phases that can be characterized by spin variables defined as the sign of the state at time $`t`$, $`\sigma _t(\alpha _1\mathrm{}\alpha _n)=\text{sign}(x_t(\alpha _1\mathrm{}\alpha _n))`$.

To study the phase-ordering phenomenon of the coupled maps on fractal networks we fix the local map parameter at the value $`\mu =1.9`$ and set the initial condition as follows: if the number of cells $`(d+1)^n`$ in a lattice is even ($`d`$ odd), exactly one half of the sites are randomly chosen and assigned random values uniformly distributed on the interval $`I^+`$, while the other half are similarly assigned values on $`I^{}`$. If the number of cells in a lattice is odd ($`d`$ even), then the state of the remaining cell is assigned at random on either interval $`I^+`$ or $`I^{}`$.

The statistical properties of the phase-ordering process on the fractal networks can be characterized by using the persistence probability $`P_t`$, defined as the fraction of cells which have not changed sign up to time $`t`$ Derrida . Figure 1 shows $`P_t`$ as a function of time for the GSG embedded in $`d=3`$, for several values of the coupling parameter. For some ranges of the coupling, the persistence saturates in a few iterations, while for some other ranges of $`ϵ`$, $`P_t`$ reaches its saturation value more slowly. In contrast, in regular Euclidean lattices the persistence saturates for small couplings, while it decays algebraically in time for coupling strengths greater than some critical value Chate .

The domains formed by the two phases on fractal lattices reach a frozen configuration for all values of the coupling $`ϵ[0,1]`$. Figure 2 shows the asymptotic patterns of the phase separation process on a GSG embedded in $`d=3`$ at level of construction $`n=3`$, for several values of the coupling. Note that the configuration of the blocked phase domains changes as the coupling is varied. The domain configurations can be characterized by the fraction of sites in a given phase that have $`k`$ neighbors in that same phase at time $`t`$, denoted by $`F_t(k)`$, with $`k=0,1,\mathrm{},d+1`$. For example, consider the pattern displayed in Fig. 2(a) where there are several sites in one phase, as those indicated by arrows, having all of their four neighbors in the opposite phase. Thus the asymptotic fraction $`F_{\mathrm{}}(0)`$ is greater than zero in this case. In Fig. 2(b), as $`ϵ`$ increases, $`F_{\mathrm{}}(0)`$ becomes zero, but $`F_{\mathrm{}}(1)`$ is finite since there are sites in one phase, as those signaled by arrows, having just one neighbor in that same phase.

The relationship between the asymptotic behavior of the persistence and the frozen domain configurations on fractals becomes manifest in Fig. 3(a), which shows the saturation value of the persistence, $`P_{\mathrm{}}`$, as well as the fractions $`F_{\mathrm{}}(0)`$, $`F_{\mathrm{}}(1)`$ and $`F_{\mathrm{}}(2)`$, as functions of the coupling parameter for the GSG embedded in Euclidean dimension $`d=3`$. Several discontinuities are observed in the curve of $`P_{\mathrm{}}`$ in Fig. 3(a). The first discontinuity of $`P_{\mathrm{}}`$ occurs at the value of the coupling $`ϵ_o=0.24`$ where the fraction $`F_{\mathrm{}}(0)`$ vanishes. This implies that local blocked configurations where a site has no neighbors in its same phase, as those indicated in Fig. 2(a), become unstable at the value of coupling $`ϵ_o`$. The second discontinuity of $`P_{\mathrm{}}`$ takes place at the value of coupling $`ϵ_1=0.30`$ where $`F_{\mathrm{}}(1)`$ becomes zero, and it is related to the loss of stability of local frozen domain configurations as those signaled in Fig. 2(b). In addition, $`P_{\mathrm{}}`$ reaches a minimum at the value of coupling $`ϵ_2=0.36`$, where the fraction $`F_{\mathrm{}}(2)`$ decays to zero. For $`ϵ>ϵ_2`$ the domains of the two phases grow in size as seen in Fig. 2(c) and 2(d). The phase domains also form faster, reducing the number of phase switching of the elements and therefore producing an increment in the saturation values $`P_{\mathrm{}}`$ up to a maximum occurring at $`ϵ_{max}=0.72`$.

A similar behavior is observed for all the fractal networks embedded in different Euclidean dimensions $`d`$. Figure 3(b) shows the saturation persistence $`P_{\mathrm{}}`$ as a function of $`ϵ`$ for GSG associated to different $`d`$. The discontinuities on each curve are related to the loss of stability of the configurations where one local element has a majority of its neighbors in the opposite phase. The number of those configurations is $`J=\text{nint}((d+1)/2)`$, where $`\text{nint}(x)`$ rounds $`x`$ to its nearest integer. These $`J`$ discontinuities take place at increasing values of the coupling $`ϵ_k`$ for which the asymptotic fractions $`F_{\mathrm{}}(k)`$ vanish, with $`k<J`$. The first discontinuity is related to the vanishing of the fraction $`F_{\mathrm{}}(0)`$ at the value $`ϵ_o=0.24`$, independently of the embedding dimension $`d`$. This independence is due to the presence of the normalization factor $`(d+1)^1`$ in the coupling term in Eq. (1) for all the lattices. For embedding dimensions $`d`$ odd, there occur a minimum of $`P_{\mathrm{}}`$ when the fraction $`F_{\mathrm{}}(J)`$ becomes zero. In those cases the local configurations losing stability are those consisting of a site with half of its neighbors in one phase and the other half in the opposite phase. These local configurations are symmetric, in the sense that a change of phase of that site does not alter the phase composition of its neighborhood. These symmetric configurations can not happen in GSG associated to $`d`$ even, and for those networks the minimum of $`P_{\mathrm{}}`$ coincides with the the last discontinuity that takes place when $`F_{\mathrm{}}(J1)`$ vanishes.

Note that for each fractal lattice there is a value of the coupling $`ϵ_{max}`$ where a maximum of $`P_{\mathrm{}}`$ is observed. The origin of such maximum can be analyzed through the average fraction of neighboring pairs that have opposite phases at time $`t`$, defined as

$$G_t=1\frac{1}{N(d+1)}\underset{i=1}{\overset{N}{}}\left(\underset{j𝒩_i}{}\delta (\sigma _t(i),\sigma _t(j))\right),$$
(3)
where $`\delta (\sigma _t(i),\sigma _t(j))=1`$ if $`\sigma _t(i)=\sigma _t(j)`$, and $`\delta (\sigma _t(i),\sigma _t(j))=0`$ otherwise. In Fig. 4 we show $`G_1`$ as a function of the parameter $`ϵ`$ for different fractal lattices. The maximum of $`G_1`$ for each lattice takes place at the value $`ϵ_{max}`$ at which the corresponding curve of $`P_{\mathrm{}}`$ reaches a maximum. When $`G_1`$ is maximum the probability that there exist cells in one phase having all of their neighbors in the opposite phase at the first iteration is also maximum. Therefore, for the value of coupling $`ϵ_{max}`$ there is a greater chance that at the next iteration such cells change phase instead of their neighbors, yielding more stable domains with fewer elements of the network having to switch their initial phase. Consequently, the persistence probability, that measures the number of elements that have not changed phases, is maximum at that time. Since the domains that are being formed are the most stable, in successive times the persistence will sustain a maximum value for the value of coupling $`ϵ_{max}`$ corresponding to each lattice. For $`ϵ>ϵ_{max}`$ the coupling is strong enough to induce transient changes in the phase of elements having the majority of their neighbors in that same phase and therefore producing lower values of $`P_{\mathrm{}}`$.

The local effect captured by the quantity $`G_1`$ in fractal lattices also appears in regular Euclidean lattices, although the asymptotic behavior of the persistence is different in those two network topologies. Figure 4 includes the calculation of $`G_1`$ as a function of $`ϵ`$ for a two-dimensional regular lattice; the maximum of $`G_1`$ in this case occurs at $`ϵ_{max}=0.67`$. This is the critical value of the coupling parameter found in Ref. Chate , after proper normalization, for the phase ordering transition in the scaling behavior of the persistence in a two dimensional Euclidean lattice. At $`ϵ=0.67`$ the blocked states in the two-dimensional regular lattice give place to growing domains of the two phases, separated by a continuous interface. The interface motion is driven by curvature effects that cause changes of phase in many elements of the system and therefore a temporal decay in the persistence probability Chate2 .

In general, in regular Euclidean lattices the ratio between the length of the interface to the size of a domain decays as $`r^1`$, where $`r`$ is the average radius of the domain. On the other hand, in the fractal networks the interface consists of a few disconnected cells separating large domains and no curvature can be defined. Because of the self-similarity of the structure the number of elements in a domain grows as $`r^{d_f}=(d+1)^l`$, where $`l`$ is an integer, while the size of the interface is of the order of $`(d+1)`$, independently of the domain sizes. Thus the ratio of the interface to domain size decreases as $`r^{d_f(l1)/l}`$ in the fractals. This decay is much faster than in regular Euclidean lattices and accelerates with increasing domain size, forming stable separated phase domains. As a consequence, domains always freeze on the fractal networks and the phase transition observed in the temporal behavior of the persistence in regular Euclidean lattices does not occur in fractals.

In summary, we have found that phase domains in chaotic maps coupled on fractal networks always reach a frozen configuration, causing the saturation of the persistence in time for all values of the coupling parameter, in contrast to the phase-transition behavior of the persistence observed in Euclidean regular lattices. The fractal nature of the spatial support is also reflected in the discontinuities observed in the $`P_{\mathrm{}}`$ vs. $`ϵ`$ curves in Fig. 3. The phase configurations of the local neighborhoods have similar transient manifestations in fractal networks and in Euclidean regular lattices, as seen in the emergence of a maximum of $`G_1`$ at a value of coupling $`ϵ_{max}`$. However, the asymptotic and global properties of the phase ordering process on these two network topologies are quite different, even when the number of local connections in the neighborhood is the same, as it happens for the two-dimensional Euclidean lattice and the GSG embedded in $`d=3`$. The necessary correlations between elements for building complex collective dynamics are more likely to occur in the denser Euclidean lattices. These results suggest that the topology of the network and not the number of local connections or the dimensionality of the space determine the asymptotic collective behaviors that may emerge on networks of coupled chaotic elements.

This work was supported by Consejo de Desarrollo Científico, Humanístico y Tecnológico, Universidad de Los Andes, Mérida, Venezuela.


text:
# Hyperon Photoproduction in the Nucleon Resonance Region

The CLAS Collaboration

(May 25, 2003)

## Abstract

Cross-sections and recoil polarizations for the reactions $`\gamma +pK^++\mathrm{\Lambda }`$ and $`\gamma +pK^++\mathrm{\Sigma }^0`$ have been measured with high statistics and with good angular coverage for center-of-mass energies between 1.6 and 2.3 GeV. In the $`K^+\mathrm{\Lambda }`$ channel we confirm a structure near $`W=1.9`$ GeV at backward kaon angles, but our data shows a more complex $`s`$\- and $`u`$\- channel resonance structure than previously seen. This structure is present at forward and backward angles but not central angles, and its position and width change with angle, indicating that more than one resonance is playing a role. Rising back-angle cross sections at higher energies and large positive polarization at backward angles are consistent with sizable $`s`$\- or $`u`$-channel contributions. None of the model calculations we present can consistently explain these aspects of the data.

preprint: CMU/JLab/5-2003

Characterizing the non-strange baryon resonances is of fundamental interest in non-perturbative QCD. The masses, quantum numbers, and decay branches of the higher-mass baryon resonances have remained difficult to establish, both experimentally and theoretically. Experimentally, most information comes from the use of pion beams interacting with nucleon targets, combined with detection of one or more pions or the nucleon. For increasing masses, both the energy overlap of resonances and meson production (e.g. the $`\rho `$) make it more difficult to separate the resonance contributions. The long list of poorly-established higher-mass resonances pdg  illustrates this problem above the strangeness threshold near $`W`$=1600 MeV. Theoretically, there is an apparent over-supply of baryons predicted in quark models, the so-called “missing baryons” problem cap . Various ways have been suggested whereby dynamical effects such as di-quarks could reduce the number of states to something closer to what has been already observed klempt .

Photoproduction of non-strange resonances detected via decay into strange particles offers two benefits in this field. First, two-body $`KY`$ final states are easier to analyze than the three-body $`\pi \pi N`$ final states that dominate decays at higher masses. So, while the cross sections for strangeness production tend to be small (on the order of 1 or 2 $`\mu `$b in electromagnetic production), the energy and angular distributions are simpler. Also, the recoil polarization observables are readily accessible via hyperon decays. Second, couplings of nucleon resonances to $`KY`$ final states are expected to differ from coupling to $`\pi N`$ or $`\pi \pi N`$ final states cap . Therefore, looking in the strangeness sector casts a different light on the resonance excitation spectrum, and thus may emphasize resonances not revealed in $`\pi N`$ scattering. Some “missing resonances” may only be “hidden” by the character of the channels studied previously. To date, however, the PDG compilation pdg  gives poorly-known $`K\mathrm{\Lambda }`$ couplings for only five well-established resonances, and no $`K\mathrm{\Sigma }`$ couplings for any resonances. The most widely-available model calculation of the $`K\mathrm{\Lambda }`$ photoproduction, the Kaon-MAID code maid , includes a mere three well-established resonances: the $`S_{11}(1650)`$, the $`P_{11}(1710)`$, and the $`P_{13}(1720)`$. Thus it is timely and interesting to have additional good-quality photoproduction data of these channels to see what additional resonance formation and decay information can be obtained. Here we report the global features of our results john  which are new, and compare them to published reaction models.

Differential cross section and hyperon recoil polarization data were obtained with the CLAS clas0  system in Hall B at the Thomas Jefferson National Accelerator Facility. A beam of tagged photons from a bremsstrahlung beam spanned energies from threshold at $`E_\gamma =0.911`$ GeV ($`W`$ = 1.609 GeV) up to 2.325 GeV ($`W`$ = 2.290 GeV). The event trigger required an electron signal from the photon tagger, and at least one charged-track coincidence between the time-of-flight ‘start’ counters near the 18-cm liquid-hydrogen target and the time-of-flight ‘stop’ counters surrounding the drift chambers. Kaons were identified using momentum and time-of-flight measurements to compute their mass, and were the only particles detected in CLAS to obtain the cross sections. The $`\mathrm{\Lambda }`$ and $`\mathrm{\Sigma }^0`$ yields were separated from the background due to mis-identified pions using lineshape fits to missing-mass spectra in each of over 900 kinematic bins of photon energy and kaon angle. The results are binned in 25 MeV steps in $`E_\gamma `$ and in 18 bins in the center-of-mass angle of the kaons, $`0.9<\mathrm{cos}(\theta _K^{CM})<+0.9`$. Consistency among several variations of kaon selection cuts and background shapes was demanded in extracting the hyperon yields, with $`\chi _\nu ^2`$ always less than 1.75 and signal-to-background ratios of greater than 2.5. A hyperon missing-mass resolution of $`\sigma =6.1`$ MeV was obtained when averaged over all detection angles and photon energies. The estimated method-dependent yield uncertainties are included bin by bin in our results, and average $`6\%`$. A total of 427,000 $`K^+\mathrm{\Lambda }`$ events and 354,000 $`K^+\mathrm{\Sigma }^0`$ events were accumulated.

The acceptance and efficiency for CLAS were modeled twice, using two independent Monte Carlo models. One was a full GEANT-based simulation involving hit digitizations, while the other was a faster parametric simulation that modeled detector effects starting at the level of reconstructed tracks. The results were in good agreement overall, and analysis of the remaining variations led to an estimated global systematic uncertainty of $`7\%`$. This was the dominant systematic uncertainty in the experiment.

The photon flux was determined by integrating the tagger rate. The rate was sampled by counting hits from accidental photons in the tagger TDC’s. Photon losses due to beam collimation were determined using a separate total-absorption counter downstream of CLAS. As a check on our results, the $`p(\gamma ,\pi ^+)n`$ cross section was measured using the same analysis chain as the $`p(\gamma ,K^+)Y`$ data. The pion cross section was found to be in good agreement with the SAID said  parameterization of the world’s data between 0.6 and 1.6 GeV, albeit with an overall scale factor of 0.92. In energy and angle the variation in the ratio of CLAS to SAID pion cross sections was small, approximately $`\pm 3\%`$, compared to our overall kaon cross section uncertainties. This shows that the yield extractions, acceptance calculations, and photon flux determinations were all consistent, but the overall normalization was made to the world’s pion photoproduction data. The final global systematic uncertainties on our cross sections are $`8.2\%`$ for the $`\mathrm{\Lambda }`$ data and $`7.7\%`$ for the $`\mathrm{\Sigma }^0`$ data.

Figure 1 (top) shows the differential cross section for $`\mathrm{\Lambda }`$ hyperon photoproduction at $`W=2.0`$ GeV. It is forward peaked, as has been seen in previous experiments bon . However, we also see a backward rise in the cross section for this and similar high values of $`W`$. This can be due either to $`u`$-channel components of the reaction mechanism or to the interference of $`s`$-channel resonances. The agreement between CLAS and previous data from SAPHIR at Bonn bon  varies: generally the measurements agree within the estimated uncertainties at back angles and near threshold energies, but CLAS measures consistently larger $`K^+\mathrm{\Lambda }`$ cross sections at forward kaon angles.

The $`\mathrm{\Lambda }`$ and $`\mathrm{\Sigma }^0`$ hyperons have isospin 0 and 1, respectively, and so intermediate states leading to the production of $`\mathrm{\Lambda }`$’s can only have isospin 1/2 ($`N^{}`$ only), whereas for the $`\mathrm{\Sigma }^0`$’s intermediate states with both isospin 1/2 and 3/2 ($`N^{}`$ or $`\mathrm{\Delta }`$) can contribute. Figure 1 (bottom) shows the data for $`\mathrm{\Sigma }^0`$ production at the same $`W`$ as above, showing different trends induced by differing resonance structure. In general the $`\mathrm{\Sigma }^0`$ cross sections are less forward peaked, with less indication of a back-angle rise.

The Regge-model calculation lag1 ; lag2  shown in Figure 1 uses only $`K`$ and $`K^{}`$ exchanges, with no $`s`$-channel resonances. The prediction was made using a model that fit high-energy kaon electroproduction data well, and could be expected to reproduce the average behavior of the cross section in the resonance region. However, extrapolated down to the resonance region, the model overpredicts the size of the $`\mathrm{\Lambda }`$ cross section and underpredicts that of the $`\mathrm{\Sigma }^0`$. Since it is a $`t`$-channel reaction model, it cannot produce a rise at back angles as seen for the $`\mathrm{\Lambda }`$, and illustrates the need for $`s`$\- and $`u`$-channel contributions to understand that feature. Two hadrodynamic models based on similar effective Lagrangian approaches maid ; mart  jan  are also shown. Both emphasize the addition of a small set of $`s`$-channel resonances to the non-resonant Born terms, and differ in their treatment of hadronic form factors and gauge invariance restoration. Both were fit to the previous data from SAPHIRbon , and therefore do not agree well with our results.

Resonance structure in the $`s`$-channel should appear most clearly in the $`W`$ dependence of the cross sections. In Figure 2 (top) we show the $`K^+\mathrm{\Lambda }`$ cross section at our most forward kaon angle, showing a sharp rise from threshold up to 1.72 GeV, a slow decline, and then a structure at 1.95 GeV with a full-width of about 100 MeV. The peak in the threshold region is understood in model calculations as due to the known $`S_{11}(1650)`$, $`P_{11}(1710)`$, and $`P_{13}(1720)`$ resonances. At a moderate forward angle, shown in Figure 2 (middle), the higher mass structure near 1.95 GeV is not visible. At a moderate backward angle, shown in Figure 2 (bottom), we again see clear structure, but it is broader, centered near 1.90 GeV, and is about 200 MeV wide. These structures are prominent at forward and backward angles; for most intermediate angles the energy dependence near 1.9 GeV falls smoothly. In contrast, our corresponding measured $`\pi ^+n`$ cross sections are featureless throughout this angle and energy range john .

Indications of the structure near 1.9 GeV were first seen in data from SAPHIR bon , which was interpreted by some mart  as evidence for a “missing” resonance at this mass. Based on theoretical guidance from one particular quark model cap , an assignment of $`D_{13}(1895)`$ seemed consistent with the angular distributions. However, other groups sag ; jan  showed that the same data could be accommodated using $`u`$-channel hyperon exchanges, an extra $`P`$-wave resonance, or alternative hadronic form factors. CLAS data, which shows a structure that varies in width and position with kaon angle, suggests an interference phenomenon between several resonant states in this mass range, rather than a single well-separated resonance. This should be expected, since many $`s`$-resonances occupy this mass range. The best modeling of the backward-angle structure near 1.9 GeV is given in Ref. mart  by incorporating a $`D_{13}(1895)`$. We see, however, that the resulting fixed position and width in this model is not consistent with the variation with angle seen in the data.

The hyperon recoil polarization provides another test of reaction models. This observable is related to interferences of the imaginary parts of the resonant amplitudes with the real part of other amplitudes, including the non-resonant Born terms. Unpolarized photons on an unpolarized target can only produce hyperons that are polarized along the axis $`(\widehat{\gamma }\times \widehat{K^+})`$ normal to the production plane. The parity-violating weak decay asymmetry in hyperon decays enables us to determine this polarization by measuring the angular distribution of the decay protons. The large acceptance of CLAS made it straightforward to detect protons from the decay of hyperons in coincidence with the $`K^+`$ mesons.

Figure 3 shows the $`\mathrm{\Lambda }`$ recoil polarization as a function of $`W`$ for representative kaon angles in the backward and the forward directions. The data has been binned such that the statistical uncertainty on each data point is less than $`\pm 0.15`$. The error bars combine statistical and estimated systematic uncertainties arising from the yield extraction. Our results are generally consistent with a few older data points from SAPHIR bon , but our energy binning is finer and reveals more structure. The data shows negative polarization of the $`\mathrm{\Lambda }`$ hyperons when the kaons go forward in the center-of-mass frame and a comparably strong positive polarization when the kaons go backward.

Of the three models tested here, only the Janssen et al. jan  model (dashed line) predicts the large back-angle polarization seen in the data near 2.0 GeV. This prediction is strongly influenced by $`u`$-channel $`Y^{}`$ contributions in that model that are added to a $`D_{13}(1895)`$ $`s`$-channel component. At the forward angle, however, this model does not perform better than the other hadrodynamic calculation or the Regge-based model. The positive back-angle polarization arises in the Kaon-MAID maid  calculation from the presence of a $`D_{13}(1895)`$, but the strength is too small. The Regge-based model lag1 ; lag2  can provide only very weak back-angle polarization, since by construction it has only a $`t`$-channel (forward-angle) production mechanism, leading here to the wrong sign. Near threshold energy the hadrodynamic models show the most structure due to interference of the known resonances cited earlier; but here our data has limited precision and cannot distinguish among these models.

In summary, we present results from an experimental investigation of hyperon photoproduction from the proton in the energy range where nucleon resonance physics should dominate. Our $`K^+\mathrm{\Lambda }`$ cross section data reveals an interesting $`W`$-dependence: double-peaked at forward and backward angles, but not at central angles. The structure near 1.9 GeV shifts in position and shape from forward to backward angles. This finding cannot be explained by a $`t`$-channel Regge-based model or by the addition of a single new resonance in the $`s`$ or $`u`$ channel. Our polarization data shows large values of polarization that change from negative values at forward angles to positive for backward kaon angles. Since a $`t`$-channel Regge model is unable to explain the backward, positive polarization, it appears that additional $`s`$\- or $`u`$-channel resonances are needed to explain the data. Our results show that hyperon photoproduction can reveal resonance structure previously “hidden” from view. Comprehensive partial wave analysis and amplitude modeling for this data can therefore be hoped to firmly establish the mass and possibly the quantum numbers of these states.

We thank staff of the Accelerator and the Physics Divisions at JLab. Major support came from the U.S. DOE and NSF, INFN Italy, the French atomic energy agency, and the Korean S&E Foundation. The Southeastern Universities Research Association (SURA) operates the Jefferson Lab for the United States DOE under contract DE-AC05-84ER40150.