text:
# Long-lived oscillons from asymmetric bubbles: existence and stability

## I Introduction

The existence of classical field configurations exhibiting soliton-like properties is not only an interesting consequence of nonlinear effects in field theory, but also an important ingredient in the understanding of nonperturbative effects in particle physics rajaraman . Apart from one-dimensional kinks, examples of these structures in 3+1 dimensions include nontopological solitons (NTS) lee  and $`Q`$ balls coleman : in both cases, choosing a simple harmonic time-dependence for the scalar field allows one to obtain a static solution of the field equations, describing a spherically-symmetric configuration which, for a range of parameters, may be the lowest energy state. Such configurations may be found in extensions of the standard model, supersymmetric or not, as has been suggested recently nts-obs . They may be sufficiently stable as to allow for a quantization procedure and form a legitimate bound state (see e.g. lee ). NTSs and Q-balls have also been of great interest to applications of particle physics to the early universe, often being proposed as possible candidates for dark matter nontop-cosm . A time-dependent, long-lived, and localized configuration in 3+1 dimensional $`\varphi ^4`$ scalar field theory was re-discovered and thoroughly studied by one of us and collaborators some years ago osc1 ; osc2 . It was shown that these configurations, named in Ref. osc1  “oscillons,” naturally arise from collapsing unstable spherically-symmetric bubbles in models with symmetric and asymmetric double-well potentials, being mainly characterized by a rapid oscillation of the field at the bubble’s core. Their relevant feature is that, albeit not strictly periodic, they possess a very long lifetime, of order $`10^310^4m^1`$ ($`\mathrm{}=c=1`$ in this work). Oscillons may be thought of as the higher-dimensional cousins of one-dimensional breather states found from kink-antikink bound states breathers . Just as kink-antikink pairs may be thermally or quantum-mechanically nucleated through nonperturbative processes, so may oscillon states, although here the calculation must be done in real and not Euclidean time. More recently, Gleiser and Sornborger investigated whether oscillons are present in 2+1 dimensions, finding not only that they do exist, but also that their lifetime is at least of order $`10^7m^1`$ sornb . Motivated by this result, in the present paper we investigate two important related questions: first, if oscillons still appear during the collapse of asymmetric – as opposed to symmetric – initial configurations; and, second, if they are stable against small angular and radial fluctuations. Since this implies that we will be dealing with a higher-dimensional parameter space, we restrict ourselves here to 2+1 dimensional oscillons. Apart from being of interest in their own right, we expect that our results will be indicative of the behavior of 3+1 dimensional oscillons. We also note that it should be quite easy to build oscillons from more complicated field theories, including interactions between the “oscillon” field and other scalar or fermionic fields. The robustness of these configurations, as demonstrated here, should provide enough motivation for a careful search of such generalized oscillons (and possibly more realistic) in the near future.

The paper is organized as follows. By means of a numerical scheme suited to tackle long-lived configurations (described in the Appendix), in Sec. II we show that oscillons quickly appear during the collapse of most elliptically deformed bubbles and, moreover, that they are all circularly symmetric and extremely long-lived, leaving no trace of the initial asymmetry. This suggests that oscillons can be understood as attractors in field-configuration space, ordered spatio-temporal structures that emerge during the nonlinear evolution of a wide variety of initial configurations. In fact, in 2+1 dimensions, the attractor basin is quite deep, as was initially hinted in Ref. sornb  and will be further shown here. We then move on to study, through a dynamical approach, whether these symmetric configurations are stable against small asymmetric perturbations, finding no indication of spectral instability (Sec. III). We conclude in Sec. IV summarizing our results and pointing out future avenues of research.

## II Oscillons from asymmetrical bubbles

The Lagrangian density for our 2+1 dimensional scalar field theory is:

$$=\frac{1}{2}(_\mu \varphi )(^\mu \varphi )\frac{\lambda }{4}\left(\varphi ^2\frac{m^2}{\lambda }\right)^2,$$
(1)
with $`\mu =0,1,2`$. We introduce dimensionless variables by rescaling the coordinates and the field as $`x_\mu ^{}=x_\mu m`$ and $`\varphi ^{}=\varphi \sqrt{\lambda }/m`$ (henceforth we drop the primes). The energy and the equation of motion are

$$E[\varphi ]=\frac{m^2}{\lambda }\text{d}^2x\left[\frac{1}{2}(\varphi /t)^2+\frac{1}{2}(\varphi )^2+\frac{1}{4}(\varphi ^21)^2\right],$$
(2)
and

$$\frac{^2\varphi }{t^2}=^2\varphi (\varphi ^3\varphi ),$$
(3)
respectively. So far, all previous studies have obtained oscillons from symmetric initial configurations, with either thick or thin walls (Gaussian or tanh profiles, respectively). We will restrict our investigation to Gaussian initial profiles, as these proved to be the most interesting in 2+1 dimensions (cf. sornb ). For convenience, we will follow Ref. sornb  and restrict the initial field configuration to interpolate between the two minima of the potential. Of course, one could select different values for the initial value of the field at the core \[$`\varphi (t=0,r=0,\theta )`$\]: as it was argued in Ref. osc2 , as long as the value of the field at the core probes the nonlinearity of the potential, and the initial configuration has an energy above the “plateau” energy (the energy of the oscillon configuration), oscillons are bound to appear.

The asymmetry in the initial field configuration is introduced by means of an elliptical deformation:

$$\varphi (\rho ,\theta ,0)=2\mathrm{exp}\left[\frac{r^2(\rho ,\theta )}{R_0^2}\right]+1,$$
(4)
where:

$$r^2(\rho ,\theta )=\frac{\rho ^2}{1\epsilon ^2\mathrm{cos}^2\theta },$$
(5)
with $`\epsilon :[0,1)`$ the bubble eccentricity, $`R_0`$ the bubble “radius” and $`\rho ,\theta `$ polar coordinates (notice that this expression reduces to the usual symmetric ansatz when $`\epsilon =0`$). We note in passing that a similar parameterization was adopted in the study of eccentric pulsons in the sine-Gordon theory chris .

In order to measure the asymmetry of the field configuration, we take advantage of the effective radius $`R_{\text{eff}}`$ defined in Ref. osc2  and introduce an “effective radial dispersion” $`\mathrm{\Delta }R_{\text{eff}}`$ (see below). The effective radius could be recast in the form:

$$R_{\text{eff}}(t)R(\theta ,t)_\theta ,$$
(6)
where $`_\theta =(2\pi )^1_0^{2\pi }\text{d}\theta `$ is an angular average and $`R(\theta ,t)`$ is a $`\theta `$-dependent effective radius defined by (compare Eq. (28) of Ref. osc2 ):

$$R(\theta ,t)=\frac{_0^{R_\text{s}}\text{d}\rho \rho ^2\left[\frac{1}{2}(\varphi )^2+V(\varphi )\right]}{(2\pi )^1_0^{2\pi }\text{d}\theta _0^{R_\text{s}}\text{d}\rho \rho \left[\frac{1}{2}(\varphi )^2+V(\varphi )\right]},$$
(7)
with $`R_\text{s}`$ the “shell radius” (i.e. the radius within which we integrate all the quantities relative to the oscillon). This notation allows us to construct angular averages and dispersions in a manner analogous to time averages. We therefore define the effective radial dispersion as the relative root-mean-square deviation from the average radius:

$$\mathrm{\Delta }R_{\text{eff}}(t)\frac{\sqrt{R^2(\theta ,t)_\theta R(\theta ,t)_\theta ^2}}{R(\theta ,t)_\theta }.$$
(8)
As shown in Figure 1, the above quantity is indeed a good measure of asymmetry, i.e. it increases monotonically with $`\epsilon `$. It is approximately independent of the bubble size, being the limiting case $`\mathrm{\Delta }R_{\text{eff}}=0`$ an indication of a symmetrical state (though not exactly zero on a lattice due to its finite resolution). We now turn to the presentation of the main numerical experiments obtained by solving Eq. (3) with the initially asymmetrical bubbles (4) for eccentricities $`\epsilon `$ ranging from $`0.1`$ to $`0.8`$ and $`R_0`$ from $`2.0`$ to $`5.0`$. This investigation can easily be extended to greater values of $`R_0`$, although this would require much longer computational times without generating results of much physical interest. \[The computational time $`𝒞_{\text{time}}`$ is proportional to $`L^2`$, which in turn is proportional to $`(R_0)^2`$, see Appendix\].

Figures 2-5 show the time evolution of the total energy within a shell of radius $`R_\text{s}=5R_0`$ surrounding the initial configuration for different values of $`R_0`$ and $`\epsilon `$. It is seen that, in general, initially asymmetric bubbles tend to decay into coherent field configurations with an approximately constant energy plateau similar to those found in Ref. sornb , which focused on the evolution of symmetric configurations. With the help of the effective radial dispersion $`\mathrm{\Delta }R_{\text{eff}}`$, we can investigate whether these configurations correspond to “excited” states of an oscillon (i.e., non-spherically symmetric configurations analogous to an excited state of a hydrogen atom for $`\mathrm{}0`$) or if the bubble asymmetry is completely lost and the system decays into a “ground” (i.e. symmetric) state.

In Fig. 6 we show the time evolution of $`\mathrm{\Delta }R_{\text{eff}}`$ for the $`R_0=4.0`$ case. It is clearly seen that the initially asymmetric bubble decays into a $`\mathrm{\Delta }R_{\text{eff}}=0`$, symmetric configuration, after a brief asymmetric pulsation. A similar pattern was observed for all initial radii investigated here, suggesting that whenever an oscillon stage is set, the resulting configuration is circularly symmetric. The peaks in this figure also suggest that $`\mathrm{\Delta }R_{\text{eff}}`$ might follow a scaling law for different values of $`\epsilon `$. Indeed, a collapse of these curves using the scaling $`\mathrm{\Delta }R_{\text{eff}}(\epsilon ,t)=\epsilon ^\alpha f(t),`$ where $`\alpha `$ is a real constant and $`f(t)`$ is a function of time only, is shown in Figure 7.

As the reader must have noticed, an intriguing feature of these results is the presence of some “instability windows” for some values of $`\epsilon `$. These can be observed here in the cases $`R_0=3.0`$ and $`R_0=5.0`$. Thus, oscillons do not always appear as the asymmetric configurations decays. A finer investigation of the parameter space for the elliptical deformations, generalizing what was done in great detail for spherically-symmetric 3+1 dimensional oscillons honda , will quite possibly reveal a very rich and detailed substructure of stable and unstable windows. It is important to stress that once $`\mathrm{\Delta }R_{\mathrm{eff}}0`$ (cf. fig. 6), the field does settle into an oscillon, as the phase space portrait of Fig. 8 exemplifies. This justifies our earlier claim that oscillons are attractors in field-configuration space.

A crucial step not yet studied is the stability analysis against small but arbitrary asymmetric perturbations. The oscillon stability with regard to these perturbations is fundamental for the computation of quantum corrections around the classical solution rajaraman ; lee  and is the subject of the next session.

## III Stability analysis

Our task now is to investigate whether symmetric oscillons are stable against small radial and angular perturbations $`\delta \varphi (\rho ,\theta ,t)`$, i.e. to probe the “linear stability” of the oscillon. We are unaware of any previous study where the stability of a field configuration with a time-dependent amplitude has been tested against small fluctuations. So far, stability investigations have been restricted to either time-independent configurations, such as bounce solutions coleman2 , or to configurations with a linear time-dependence in the phase, such as $`Q`$-balls coleman ). The stability analysis in these two situations is greatly simplified by the fact that the dynamical equation dictating the behavior of the fluctuations is separable into its spatial and temporal parts; the resulting problem reduces to finding the eigenvalues $`\omega _n^2`$ of a time-independent operator (an alternative approach based on the so-called “Bogomolnyi bound” bogomolnyi  can also be applied in the time-independent case, see e.g. calberto ). The existence of at least one negative eigenvalue $`\omega _n^2`$ (and thus of a complex eigenfrequency $`\omega _n`$) signals the presence of an exponentially-growing instability jackiw  (a well-known example is the so-called “bounce” solution, which we will investigate further below coleman2 ). The present problem, however, is not amenable to such treatment due to the anharmonic time-dependence of the oscillon; we must consider both the space- and time-dependence of the background field, making the stability analysis of oscillon-type configurations considerably more challenging both analytically and numerically, as we will now discuss.

In order to appreciate these difficulties, let us write (in polar coordinates) the linearized equation of motion that follows from Eq. (3) through the substitution $`\varphi (\rho ,\theta ,t)\varphi _0(\rho ,t)+\delta \varphi (\rho ,\theta ,t)`$, where $`\varphi _0`$ is the symmetric oscillon solution and $`\delta \varphi \varphi _0`$ is the perturbation, i.e.,

$$\frac{^2}{t^2}\delta \varphi =\frac{^2}{\rho ^2}\delta \varphi +\frac{1}{\rho }\frac{}{\rho }\delta \varphi +\frac{1}{\rho ^2}\frac{^2}{\theta ^2}\delta \varphi (3\varphi _0^21)\delta \varphi .$$
(9)
Here one might be tempted to separate the variables as $`\delta \varphi R(\rho )\mathrm{\Theta }(\theta )T(t)`$. However, the resulting equations show that one cannot get rid of the simultaneous radial and temporal dependence of the background configuration, $`\varphi _0(\rho ,t)`$ 111In Ref. osc1  it was assumed that the time-scale of possible unstable fluctuations would be faster than typical oscillon time-scales. However, such “adiabatic” assumption is too restrictive, as the fluctuation time-scales are a priori unknown.. This situation should be contrasted to the usual case where $`\varphi _0`$ is a time-independent solution, e.g. the bounce, or to the case where the time-dependence of $`\varphi _0`$ is in a phase factor $`\mathrm{exp}(i\omega t)`$, and thus immediately eliminated in the full equation of motion coleman . A considerable simplification can nevertheless be accomplished by writing $`\delta \varphi \mathrm{\Phi }(\rho ,t)\mathrm{\Theta }(\theta )`$, isolating at least the angular part of the problem. Performing such substitution gives the pair of equations,

$$\ddot{\mathrm{\Phi }}\mathrm{\Phi }^{\prime \prime }\frac{1}{\rho }\mathrm{\Phi }^{}+\left(U_0(\rho ,t)+\frac{\mathrm{}^2}{\rho ^2}\right)\mathrm{\Phi }=0,$$
(10)
and

$$\frac{d^2\mathrm{\Theta }(\theta )}{d\theta ^2}=\mathrm{}^2\mathrm{\Theta }(\theta ),$$
(11)
with overdots and primes indicating time and radial derivatives, respectively. Here $`U_0(\rho ,t)3\varphi _0^21`$ and $`\mathrm{}`$ is a separation constant. The solutions for $`\mathrm{\Theta }`$ are trivial, viz. $`\mathrm{\Theta }\mathrm{exp}(\pm i\mathrm{}\theta )`$, and by requiring $`\delta \varphi `$ to be single-valued we have $`\mathrm{}=0,1,2,\mathrm{}`$. Our original (2+1)-dimensional problem, Eq. (9), reduces therefore to solving the above (1+1)-dimensional one, Eq. (10), since the time-dependence is present only in $`\mathrm{\Phi }`$.

Our goal is to probe the linear stability of the oscillon by solving Eq. (10) for arbitrary initial conditions. The strategy is to monitor the time evolution of the perturbations $`\delta \varphi `$, which should grow without bounds in the case of a linearly unstable configuration merkin ; morrison . An obvious limitation of this approach is that it is impossible to scan all initial values of perturbations, viz. $`\delta \varphi (\rho ,0)`$ with $`\mathrm{}=0,1,2,\mathrm{}`$. The method thus is only indicative of stability, not being able to provide conclusive proof. The more thorough the search, the more one is guaranteed to show stability, at least against most types of perturbation. This unavoidable limitation should be contrasted with the simpler case for time-independent background configurations based on a harmonic decomposition $`T(t)\mathrm{exp}(i\omega _nt)`$ (see e.g. coleman ; jackiw ), where the existence of exponentially unstable modes is clearly related to imaginary eigenvalues. However, we would like to point out a limitation of the spectral method that is often overlooked. By restricting the analysis to an exponential time-dependence, as in $`T(t)`$ above, one can obtain only spectral instabilities of a configuration, leaving aside other possible forms of instabilities, for example, linear (or power-law) ones. In other words, a system that is spectrally stable may still be unstable against slower growing modes morrison ; lichtenberg . Since we are here essentially watching the full time-dependence of $`\delta \varphi `$, we should be able to detect any sort of instability by observing its long-time behavior, although in practice the infinite-time limit or a complete scan of possible fluctuations cannot be achieved numerically. Fortunately, we shall shortly see that typical spectral instabilities (such as that of the bounce) do not require a long-time integration or a very wide search, being therefore bound to be observed through our method. Before we do so, it is worth testing the reliability of the numerical implementation itself.

### III.1 Linear test

The first step is to compare the numerical solution of Eq. (10) with a closed-form analytical one in order to prepare and test our numerical implementation, since the singular behavior at the origin requires a careful treatment. This can be done most easily by setting $`U_0(\rho ,t)=0`$, in which case Eq. (10) becomes linear and separable, with $`\mathrm{\Phi }=R(\rho )T(t)`$. This gives $`T\mathrm{exp}(\pm i\omega t)`$, where $`\omega `$ is a separation constant, and the equation for $`R`$,

$$R^{\prime \prime }+\frac{1}{\rho }R^{}+\left(\omega ^2\frac{\mathrm{}^2}{\rho ^2}\right)R=0,$$
(12)
which we recognize as Bessel’s equation. By requiring regularity at the origin and $`\dot{\mathrm{\Phi }}(\rho ,0)=0`$, the solution can be written as

$$\mathrm{\Phi }_{\mathrm{}}(\rho ,t)=_0^{\mathrm{}}𝑑\omega f(\omega )J_{\mathrm{}}(\omega \rho )\mathrm{cos}(\omega t),$$
(13)
where $`f(\omega )`$ is determined by the initial condition. Choosing $`\mathrm{\Phi }(\rho ,0)=AJ_{\mathrm{}}(a\rho )`$ we have

$$f(\omega )=A\omega _0^{\mathrm{}}𝑑\rho \rho J_{\mathrm{}}(a\rho )J_{\mathrm{}}(\omega \rho )=A\delta (\omega a),$$
(14)
and therefore

$`\mathrm{\Phi }_{\mathrm{}}(\rho ,t)`$ $`=`$ $`A{\displaystyle _0^{\mathrm{}}}𝑑\omega \delta (\omega a)J_{\mathrm{}}(\omega \rho )\mathrm{cos}(\omega t)`$ (15)
$`=`$ $`AJ_{\mathrm{}}(a\rho )\mathrm{cos}(at).`$
The above solution maintains its shape but oscillates harmonically with period $`2\pi /a`$. We have verified that our numerical implementation reproduces correctly this analytical solution for various values of $`a`$ and $`\mathrm{}`$.

### III.2 The bounce

As a first application of our method we investigate the stability of the so-called “bounce” solution coleman2 , which is guaranteed to be spectrally unstable in any dimension greater than (1+1) due to Derrick’s theorem derrick . (In fact, Coleman has showed that only one negative eigenvalue exists coleman3 ). Should our method be reliable, the solution $`\delta \varphi `$ for the case where $`\varphi _0`$ is a bounce solution will grow exponentially at late times, indicating the presence of an unstable mode.

The bounce solution $`\varphi _b(\rho )`$ is the $`O(2)`$-symmetric static configuration that solves the equation

$$\frac{d^2\varphi }{d\rho ^2}+\frac{1}{\rho }\frac{d\varphi }{d\rho }=V^{}(\varphi ),$$
(16)
with the asymmetric potential

$$V(\varphi )=\frac{\beta }{2}\varphi ^2\frac{1}{3}\varphi ^3+\frac{1}{4}\varphi ^4.$$
(17)
In order to detect the instability, we solved Eq. (10) with $`U_0(\rho )=\beta 2\varphi _b(\rho )+3\varphi _b^2(\rho )`$ and various initial conditions, sweeping the lattice at every time step to find the maximum value of the perturbation, $`\mathrm{\Phi }_{max}`$. In Fig. 9, we show our results for $`\beta =0.011`$ \[the initial conditions are Eq. (18) with $`n=0,1`$ and $`m=2`$, and Eq. (20) for $`m=2`$, both with $`\mathrm{}=0`$\]; one can clearly identify the exponential growth of $`\delta \varphi `$ even at early times $`t<100`$. Also shown is the slope of the curve, which should match the unstable eigenvalue $`\omega _n`$ obtained with the usual spectral method. \[We have attempted to obtain such eigenvalue by solving numerically the associated Schrödinger-like equation. However, in two spatial dimensions the severe singularity at the origin causes a numerical instability which we were unable to control even with sophisticated methods liu . Since this is not the focus of this paper, we will leave this question aside.\]

### III.3 The oscillon

We are now ready to apply our method to the stability of the oscillon, which was obtained here by solving Eq. (3) with the symmetrical version ($`\epsilon =0`$) of the ansatz Eq. (4). We have essentially followed the same procedure described above for the bounce, but now evolving both $`U_0`$ and $`\mathrm{\Phi }`$ in Eq. (10). Since the dimensionality of the configuration space is infinite, we chose arbitrarily the initial profiles of the fluctuations $`\mathrm{\Phi }(\rho ,t_0)`$, with the only constraint that they should vanish at $`\rho \mathrm{}`$ to ensure localization around the oscillon. The time $`t_0`$ was chosen to be about $`200`$, since that is roughly when the initial bubbles have just decayed into an oscillon (cf. Figs. 2-5). Some examples of the initial configurations investigated here are

$`\mathrm{\Phi }_{\mathrm{}}(\rho ,t_0)`$ $`=x^n\mathrm{exp}\left({\displaystyle \frac{\rho ^m}{R_0^m}}\right),`$ (18)
$`\mathrm{\Phi }_{\mathrm{}}(\rho ,t_0)`$ $`=\mathrm{sin}\left({\displaystyle \frac{2n\pi \rho }{R_0}}\right)\mathrm{exp}\left({\displaystyle \frac{\rho ^m}{R_0^m}}\right),`$ (19)
$`\mathrm{\Phi }_{\mathrm{}}(\rho ,t_0)`$ $`=J_{\mathrm{}}(\rho )\mathrm{exp}\left({\displaystyle \frac{\rho ^m}{R_0^m}}\right),`$ (20)
for various integer values of $`n,m`$ and $`\mathrm{}`$ (namely, $`n=0,1,2`$, $`m=2,3`$ and $`\mathrm{}=0,1,2,3`$).

In Figure 10 we show a typical outcome of our search. In all cases investigated we found that the fluctuations $`\mathrm{\Phi }_{max}`$ are bounded from above, as one would expect from a linearly stable configuration. We conclude that if, indeed, there are any unstable modes, they are sufficiently slow-growing to justify the use of the oscillon as a stable bound state. \[We have integrated the linearized equations of motion up to $`t=10^4`$; see also discussion below\]. Note that the large amplitude of $`\mathrm{\Phi }`$, e.g. $`1000`$ in Fig. 10, does not mean that the condition $`\delta \varphi \varphi _0`$ is violated: since the resulting equation is linear, any solution can always be rescaled without changing its shape by choosing a different constant prefactor for the initial conditions.

Although not as systematic and transparent as the investigation above, another approach to check the stability of oscillons is to superimpose the perturbation to the full (2+1)-dimensional oscillon dynamics discussed in Sec. II. One can then probe the oscillon stability simply by checking the persistence of the energy plateau: if the added energy from the perturbations is radiated away, the oscillon is stable. Due to the dimensionality of the problem, the numerical treatment is quite more challenging than the one use above within the linear method. Nevertheless, we have investigated the stability of oscillons against superimposed fluctuations for similar initial conditions. In Fig. 11 we present the outcome of a particular choice of initial condition for three different initial radii. The results are consistent with the previous stability analysis, as can be seen by the persistence of the energy plateau.

On the basis of our extensive search with many different initial conditions and long integration times, we find it very unlikely that an exponentially-growing mode exists. If it does, it would be either very small and/or related to a very “rare” excitation; the oscillon configuration would still be stable for large times and could be considered a legitimate (or at least a very long-lived) bound state in semi-classical quantization.

## IV Conclusions

We have investigated, in 2+1 dimensions, two key questions concerning the properties of oscillons – time-dependent, localized field configurations that emerge during the deterministic evolution of $`\varphi ^4`$ models. First, we have shown that initially asymmetric configurations evolve, for a wide range of elliptic deformations, into symmetric oscillons states. Thus, oscillons are not just particular to symmetric initial states. This result led us to propose that oscillons are attractors in field configuration space, with a very deep attractor basin, at least in 2+1 dimensions. Second, we have shown that oscillons are stable against a wide range of asymmetric small perturbations. This result was obtained by two distinct approaches, one solving the linearized equation for the perturbations and the other by superimposing the perturbations on the oscillons and evolving the perturbed configurations with the full equation of motion. Clearly, both methods are restricted to the choice of initial fluctuations. However, after an extensive search, we were unable to find any unstable fluctuation with either approach. To the best of our knowledge, this is the first dynamical investigation of the stability of explicitly time-dependent scalar field configurations. We expect that both of these results will carry on to 3+1 dimensions, although probably the attractor basin will be shallower in this case.

These results suggest the importance time-dependent spatio-temporal structures may have in a wide range of physical systems, from condensed matter to early universe cosmology. Although we have restricted our study to simple $`\varphi ^4`$ models, we expect, as suggested in Refs. osc2 ; sornb , that oscillons will be present whenever there is a bifurcation instability related to the negative curvature of the nonlinear potential. Oscillons will emerge in a wide variety of dynamical systems, possibly representing a bottleneck to equipartition of energy, thus delaying the approach to equilibrium.

One possible arena for oscillons in early universe cosmology is during the reheating supposed to occur after inflation. Oscillons may be thermally nucleated with a probability proportional to $`\mathrm{exp}[E_{\mathrm{osc}}/T]`$, where $`E_{\mathrm{osc}}`$ is the energy of the oscillon configuration. They will act as “entropy sinks”, confining several degrees of freedom to an ordered state, delaying the thermalization of the universe. Eventually, when they decay into radiation, they will dump more entropy to the early universe, possibly changing the final reheating temperature.

Finally, it would be interesting to compute the spectrum of quantum fluctuations around oscillon states, to investigate their effect on oscillon stability. Time-dependent bound states may have much to add to our knowledge of nonperturbative quantum field theory, which has traditionally focused on time-independent configurations, such as instantons.

###### Acknowledgements.

The authors are indebted to the referee for pointing out a limitation of the original stability analysis, which led to the present treatment. A. B. Adib acknowledges Rafael Howell for providing the numerical bounce solution for our stability analysis and for many fruitful discussions, the Physics Department at UFC for the kind hospitality during part of this work and Dartmouth College for financial support. M. Gleiser was partially supported by NSF grants PHY-0070554 and PHY-0099543, and by the “Mr. Tompkins Fund for Cosmology and Field Theory” at Dartmouth. C. A. S. Almeida was supported in part by CNPq (Brazil).

*

## Appendix A Numerical method

The integration scheme adopted here is a standard leapfrog algorithm, which ensures second-order precision in time numrec  (the spatial discretization used is fourth-order). We have adopted the “adiabatic damping method” (or simply the damping method) of Ref. sornb  together with Higdon’s first-order boundary conditions higdon . Their combined use turned out to be very effective and of easy implementation, allowing us to tackle this otherwise demanding numerical problem with current workstations.

Put briefly, the leapfrog equations with the damping method read:

$`\dot{\varphi }_{i,j}^{n+1/2}`$ $`=`$ $`{\displaystyle \frac{(1\eta _{i,j}\mathrm{\Delta }t/2)\dot{\varphi }_{i,j}^{n1/2}+\mathrm{\Delta }t[^2\varphi _{i,j}^nV^{}(\varphi _{i,j}^n)]}{1+\eta _{i,j}\mathrm{\Delta }t/2}},`$
$`\varphi _{i,j}^{n+1}`$ $`=`$ $`\varphi _{i,j}^n+\mathrm{\Delta }t\dot{\varphi }_{i,j}^{n+1/2},`$ (21)
where superscripts (subscripts) denote temporal (spatial) indices, $`\eta _{i,j}`$ is the damping function of Ref. sornb  and $`V^{}(\varphi )`$ is a first partial derivative of the potential with respect to the field. The second spatial derivatives in the Laplacian are discretized with a fourth-order scheme (to wit, $`_{xx}\varphi _{i,j}[16(\varphi _{i+1,j}+\varphi _{i1,j})\varphi _{i+2,j}\varphi _{i2,j}30\varphi _{i,j}]/12\mathrm{\Delta }x^2`$ and analogously for $`_{yy}\varphi _{i,j}`$), which gives an energy conservation of one part in $`10^3`$ for $`\mathrm{\Delta }x=0.1`$ and $`\mathrm{\Delta }t=0.06`$ (and, of course, with $`\eta _{i,j}0`$). A better energy conservation could be obtained with smaller $`\mathrm{\Delta }x`$ or $`\mathrm{\Delta }t`$, but this comes with a high price tag since, as remarked below, the computational time is inversely proportional to both $`\mathrm{\Delta }x^2`$ and $`\mathrm{\Delta }t`$. Despite this fact, with the above parameters we were able to reproduce quite accurately the results of Gleiser and Sornborger sornb . Even though the damping method is already a major improvement over more naive methods (such as huge lattices or even moving boundary conditions), for the problem at hand it is still demanding. As an example, for small oscillons of radius $`R_02.0`$, the required lattice of radial dimension $`R200`$ adopted in Ref. sornb  (and thus $`L400`$ in our square grid, where $`L`$ is the lattice edge), would already demand a total of $`N10^7`$ sites for $`\mathrm{\Delta }x=0.1`$, as opposed to the $`N10^3`$ used in the latter reference. Another aggravating fact comes from the large integration times $`\tau `$ involved in such problems \[notice that the required computational time for this problem goes roughly as $`𝒞_{\text{time}}(\tau /\mathrm{\Delta }t)N=(\tau /\mathrm{\Delta }t)(L/\mathrm{\Delta }x)^2`$\]. We note in passing that there has been some effort to find a more natural and efficient discretization for the $`\varphi ^4`$ theory which might reduce significantly the computational time of such problems speight . Motivated by this possibility, two of us have recently investigated these lattices and have found that, unfortunately, they are of limited practical use even for simple dynamical problems kinkdyn . It was seen, however, that if the above scheme is supplied with the boundary conditions of Ref. higdon :

$$\left(_t\pm _\alpha \right)\varphi |_\mathrm{\Omega }=0,$$
(22)
where $`\alpha `$ is either $`x`$ or $`y`$, then a significantly smaller lattice could be used, resulting in an energy error smaller than (or equal to) the error due to numerical energy fluctuations. \[These first-order “absorbing boundary conditions” were obtained for the rather simple (linear) wave equation. We expect, however, that the damping introduced before the boundaries could reduce the amplitude of the outgoing waves such that Eq. (3) is effectively linearized in that region, and thus that the boundary condition (22) becomes applicable\]. With regard to the example in the previous paragraph, we have found that the required lattice with this “mixed” method needs only $`L140`$ (in contrast to the former $`L400`$), such that $`N`$ (and thus $`𝒞_{\text{time}}`$) is roughly one order of magnitude smaller than the previous one (this trend is also found for greater $`R_0`$). For the sake of completeness, we quote here the parameters of the damping method used throughout our simulations (we use the same functional form for $`\eta (\rho )`$ as Ref. sornb ): $`k=0.005`$ (damping constant), $`\rho _0=10R_0`$ (initial radius of the damping) and $`\rho _{\mathrm{}}=50`$ (damping length), these latter two being defined such that $`L=2(\rho _0+\rho _{\mathrm{}})`$.

We expect that the method adopted here might be useful not only in higher dimensional systems (the two methods above do not really make any dimensional requirement), but also in other finite-domain problems not necessarily related to oscillons.


format: text

text:
# Unification of Gravity and Electro-Magnetism in Six Dimensions

## Abstract

Recently by us was proposed the model where Einstein’s equation on the brane was connected with Maxwell’s multi-dimensional equations in pseudo-Euclidean space. Based on this idea unification of 4-dimensional gravity and electromagnetism in (2+4)-space is found. In this picture photon is massless in four dimensions and obtains large mass in extra (1+1)-space normal to the brane.

PACS numbers: 04.50.+h, 03.50.De, 98.80.Cq

One of the main ideas of the brane models is that Plank’s scale $`M_{Pl}`$ can be constructed with the fundamental scale $`M`$ and the brane width $`ϵ`$ . In our previous paper  it was shown that in the brane approach Einstein’s equations can be effective as well. Gravitation was connected with the geometry of the brane embedded in pseudo-Euclidean space, where solutions of 4-dimensional Einstein’s equation could be constructed with the solutions of multidimensional Maxwell’s equations. This approach is the inverse to the standard Kaluza-Klein picture, since the vector field is fundamental and gravity is obtained from the brane geometry. In this direction geometrical unification of different interactions is easier, since Dirac equation also can be derived from the constrained Yang-Mills Lagrangian . The surface where Einstein’s equations appear is not necessary to be 4-brane. For example, it can be 5-dimensional space with exponential warp factor needed to cancel cosmological constant . Embedding of this 5-space into Euclidean (2+4)-space was found in .

The main idea of  was based on the fact that any $`n`$-dimensional Riemannian space can be embedded into $`N`$-dimensional pseudo-Euclidean space with $`nNn(n+1)/2`$ . This means that the multidimensional ’tetrad’ field $`h_\alpha ^A`$ can be expressed as a derivative of embedded functions

$$h_\alpha ^A=_\alpha \varphi ^A.$$
(1)
Capital Latin letters $`A,B,\mathrm{}`$ numerate coordinates of embedded space, while Greek indices $`\alpha ,\beta ,\mathrm{}`$ numerate coordinates in four dimensions.

The metric of bulk flat space-time (where the brane with arbitrary geometry is embedded) can be transformed to the Gaussian normal coordinates

$$ds^2=\eta _{ij}dx^idx^j+g_{\alpha \beta }(x^i,x^\nu )dx^\alpha dx^\beta .$$
(2)
Here $`\eta _{ij}`$ is the Minkowskian metric of the normal space to the brane and small Latin indices $`i,j,\mathrm{}`$ numerates extra $`(N4)`$ coordinates. The induced metric on the brane is $`g_{\alpha \beta }(0,x^\nu )`$.

Embedding theory allows us to rewrite Hilbert’s 4-dimensional action in terms of derivatives of the normal components $`\varphi ^i`$ of some multi-dimensional vector

$$S_g=M_{Pl}^2R\sqrt{g}d^4x=\eta _{ij}M_{Pl}^2\left(\mathrm{}\varphi ^i\mathrm{}\varphi ^j_\alpha _\beta \varphi ^i^\alpha ^\beta \varphi ^j\right)\sqrt{g}d^4x,$$
(3)
where $`\mathrm{}=_\alpha ^\alpha `$ is the 4-dimensional wave operator and $`g`$ is the determinant in Gaussian coordinates.

In  it was shown that to the same form (3) can be transformed Maxwell’s multi-dimensional action

$$S_A=\frac{1}{4}F_{AB}F^{AB}d^NX=\frac{1}{2}[_AA_B(^AA^B+^BA^A)2_AA^A_BA^B]d^NX,$$
(4)
where $`F_{AB}=_AA_B_BA_A`$ and $`X^A`$ are N-dimensional Cartesian coordinates. For this reduction we must transform $`X^A`$ to the Gaussian coordinates (2) and insert in (4) the exact solution of Maxwell’s equations for the brane

$`A^\alpha =\left({\displaystyle \frac{2}{\pi ϵ^2}}\right)^{(N4)/4}M{\displaystyle \underset{i}{\overset{(N4)}{}}}exp\left[{\displaystyle \frac{(x^i)^2}{ϵ^2}}\right]^\alpha \varphi ^i(x^\beta ),`$
$`A^i=\left({\displaystyle \frac{2}{\pi ϵ^2}}\right)^{(N4)/4}{\displaystyle \frac{2Mx^i}{ϵ^2}}exp\left[{\displaystyle \frac{(x^i)^2}{ϵ^2}}\right]\varphi ^i(x^\beta ).`$ (5)
Then we must integrate (4) by the extra coordinates $`x^i`$, remove boundary terms and as in  put

$$M_{Pl}^2=ϵ^{(N4)}M^2.$$
(6)
In the present paper we want to find 4-dimensional unification of gravity and Electro-Magnetism based on Maxwell’s equations in 6-dimensional pseudo-Euclidean space with the signature $`(2+4)`$.

The (2+4)-space is interest object for investigations. As it was found a long time ago Schwarzschild’s metric can be embedded in (2+4)-space ; also non-linear 15-parameter conformal group can be written as a linear Lorentz-type transformations there . For compact Kaluza-Klein models (2+4)-space was studied in  and in the context of brane models with non-factorizable geometry in .

In our previous paper  was considered only possible connection of Einstein’s equations with multi-dimensional Maxwell’s equation and was not received Maxwell’s equations on the brane. In (Unification of Gravity and Electro-Magnetism in Six Dimensions) there presents only derivatives of $`\varphi ^i`$ (which is connected with gravity) and it is clear that $`F_{\alpha \beta }`$ is zero. Now we want to insert in (Unification of Gravity and Electro-Magnetism in Six Dimensions) also non-integrable part $`a^\alpha (x^\beta )`$ and identify it with the 4-dimensional vector-potential. So, in (2+4)-space we wand to rewrite (Unification of Gravity and Electro-Magnetism in Six Dimensions) in the form

$`A^\alpha =\sqrt{{\displaystyle \frac{2}{\pi ϵ^2}}}\left[Mexp\left({\displaystyle \frac{\tau ^2}{ϵ^2}}\right)^\alpha \varphi ^\tau (x^\beta )+Mexp\left({\displaystyle \frac{\kappa ^2}{ϵ^2}}\right)^\alpha \varphi ^\kappa (x^\beta )+exp\left({\displaystyle \frac{\tau ^2+\kappa ^2}{ϵ^2}}\right)a^\alpha \right],`$
$`A^\tau =\sqrt{{\displaystyle \frac{2}{\pi ϵ^2}}}{\displaystyle \frac{2M\tau }{ϵ^2}}exp\left({\displaystyle \frac{\tau ^2}{ϵ^2}}\right)\varphi ^\tau (x^\beta ),`$ (7)
$`A^\kappa =\sqrt{{\displaystyle \frac{2}{\pi ϵ^2}}}{\displaystyle \frac{2M\kappa }{ϵ^2}}exp\left({\displaystyle \frac{\kappa ^2}{ϵ^2}}\right)\varphi ^\kappa (x^\beta ),`$
where $`\kappa `$ and $`\tau `$ are the normal to the brane space-like and time-like coordinates respectively.

Using the anzats (Unification of Gravity and Electro-Magnetism in Six Dimensions) Maxwell’s action (4) in (2+4)-space takes the form

$`S_A={\displaystyle \frac{1}{4}}{\displaystyle F_{AB}F^{AB}d^6X}=`$
$`=M^2ϵ^2{\displaystyle R\sqrt{g}d^4x}{\displaystyle \frac{1}{2\pi ϵ^2}}{\displaystyle exp\left[\frac{2(\tau ^2+\kappa ^2)}{ϵ^2}\right]\left[f_{\alpha \beta }f^{\alpha \beta }+\frac{4}{ϵ^4}(\tau ^2\kappa ^2)a_\alpha a^\alpha \right]\sqrt{g}d^4x𝑑\tau 𝑑\kappa }=`$ (8)
$`={\displaystyle \left(M_{Pl}^2R+\frac{1}{4}f_{\alpha \beta }f^{\alpha \beta }\right)\sqrt{g}d^4x},`$
where $`f_{\alpha \beta }=_\alpha a_\beta _\beta a_\alpha `$. Latin indices are summing with the 6-dimensional Minkowskian metric tensor $`\eta _{AB}`$ and Greek indices with the induced on the brane metric tensor $`g_{\alpha \beta }`$. So, 6-dimensional Maxwell’s action on the brane is reduced to the 4-dimensional Einstein-Maxwell action. As in  the brane is placed at $`\tau ^2=\kappa ^2`$ what means that it moves with the speed of light in transversal (1+1)-space.

In Gaussian coordinates (2) Maxwell’s 6-dimensional equations

$$_AF^{AB}=0,$$
(9)
after inserting there the anzats (Unification of Gravity and Electro-Magnetism in Six Dimensions) reduces to the system

$`{\displaystyle \frac{1}{\sqrt{g}}}_\alpha (\sqrt{g}f^{\alpha \beta }){\displaystyle \frac{4}{ϵ^4}}(\tau ^2\kappa ^2)a^\beta =0,`$
$`_\alpha a^\alpha =0.`$ (10)
One can see that the situation is similar to superconductors. On the brane ($`\tau ^2=\kappa ^2`$) photon is massless and obtains the large mass $`1/ϵ`$ in the direction normal to the brane.

At the end we want to note that even if we don’t want to connect $`\varphi ^i`$ with gravity the anzats (Unification of Gravity and Electro-Magnetism in Six Dimensions) can be used for the trapping of photons on the brane in the (2+4)-space.

Acknowledgements: Author would like to acknowledge the hospitality extended during his visits at the Abdus Salam International Centre for Theoretical Physics where this work was done.


text:
# Infrared Surveys for AGN

## 1.  Motivation: Active Galaxies in the Infrared

There are a number of motivations for infrared studies of AGN and for carrying out surveys for AGN at infrared wavelengths:

1. AGN unification models require a dusty, molecular obscuring screen (torus?) which is expected to emit in the infrared.
2. Cosmic X-Ray Background models require a population of highly obscured AGN which may be detectable in the infrared.
3. Infrared/Starburst emission is a frequent companion to AGN activity; IR emission originates from dust which is either associated with the active nucleus (e.g. the PG QSO sample; Sanders et al. 1989) or with circumnuclear starburst emission as in the classical Seyfert galaxies NGC 1068 & NGC 7469.
4. Luminous Infrared Galaxies frequently harbor AGN cores and may be a stage in the development of AGN from the merger of gas-rich galaxies.

### 1.1.  AGN Unification Models

The AGN Unification models (\[Antonucci 1993\]) which have been popular for some time, and which have garnered considerable observational support, require an optically-thick, dusty obscuring screen, frequently assumed to be toroidal (\[Pier & Krolik 1993\]; \[Granato & Danese 1994\]), which obscures the central engine and broad emission-line region from equatorial lines of sight. In this scenario broad-line objects (Type 1 AGN – Sy1 galaxies and classical QSOs) are objects viewed near the pole of the torus, while narrow-line objects (Type 2 AGN – Sy 2 galaxies) are viewed edge-on. Direct evidence for the existence of this obscuring screen is largely from studies of scattered broad-line emission in Sy 2 systems, but claims have been made that the mid-infrared spectra of some nearby galaxies match various torus models (\[Pier & Krolik 1993\]; \[Granato & Danese 1994\]). A critical point regarding this torus is that it will be optically thick at all wavelengths from soft x-rays through the mid-infrared. Even mid-infrared diagnostics may not reveal a highly obscured AGN. Hard x-rays and VLBI radio imaging are currently required to peer through the veil of obscuration.

### 1.2.  The Cosmic X-Ray and Infrared Backgrounds

Recent models for the Cosmic X-Ray Background (CXB) require a population of highly obscured, perhaps even compton-thick, AGN which emit primarily in the hard x-ray region in order to fit the x-ray background spectrum (\[Comastri et al. 1995\]; \[Gilli et al. 2001\]). This population increases in density with redshift. These AGN are expected to re-emit the absorbed radiation at infrared wavelengths providing a link between the CXB and the Cosmic InfraRed Background (CIRB). Recent measures of the CIRB with COBE (\[Hauser et al. 1998\]; \[Puget et al. 1996\]) show that more than half of the cosmic energy density (excluding the CMB) comes out in the infrared. An important question is then the relative importance of accretion energy due to AGN compared with stellar nucleosynthesis. Deep surveys with Chandra and XMM-Newton are resolving the CXB (\[Hornschemeier et al. 2001\]; \[Hasinger et al. 2001\]; \[Rosati et al. 2001\]) and ISO (\[Elbaz et al. 2002\]) has begun to resolve CIRB. New Surveys with SIRTF will convincingly determine relationship between the faint x-ray and infrared populations revealing a great deal about the history of our Universe.

### 1.3.  Luminous Infrared Galaxies

The most luminous galaxies in the Local Universe are Luminous Infrared Galaxies (LIGs) which emit the vast majority of their radiant power in the far-infrared between about 40–120$`\mu `$m. These are gas-rich systems which are in the late stages of collisions or mergers. Extrapolation from the properties of lower luminosity Starburst galaxies suggests that the LIGs should be active star-forming systems (see Sanders & Mirabel 1996 for a review). The LIGs also show many characterisatics of AGN and their luminosities reach values comparable to those of luminous QSOs. Much effort has been focused on whether LIGs are powered principally by Starburst or AGN activity, although both types of activity are almost certainly present. The discussion has been framed around a scenario proposed by Sanders et al. (1988) in which a merger of gas-rich disk galaxies stimulates a massive nuclear Starburst which in turn feeds a coalescing AGN core in the galaxy nucleus. As the AGN turns on, radiation pressure drives out the shroud of dust, revealing a nascent quasar. The goal must be not only to understand the dominant source of energy in LIGs, but to understand the relationship between Starburst and AGN activity and other galaxy characteristics, and to place them into an evolutionary context.

We have for some time been using VLBI techniques to attempt to understand the power sources of LIGs and to place them in an evolutionary context. Analysis of a complete LIG sample (\[Smith, Lonsdale & Lonsdale 1998\]) suggests that many LIGs may be interpreted as intense Starbursts like Arp 220 (Figure 1), but a nearly equal number must house AGN cores, as in the case of Mrk 231, whose nuclear structure implies a recent ignition ($`t<<10^6yr`$) for AGN activity. Evidence is strong that a significant number of LIGs house obscured AGN activity and schematics are consistent with this Starburst-to-AGN scenario.

Of particular interest are the extremely luminous high-z galaxies detected in recent SCUBA and other submillimeter surveys (\[Ivison et al. 2000\]). The SEDs of LIGs in the submillimeter have the unique characteristic that the positive K-correction offsets cosmological dimming for redshifts from $`z1`$ – 10 such that it is equally easy (or difficult) to detect infrared galaxies at 850$`\mu `$m over a range of high redshifts. Although only a handful of redshifts are available for submm sources owing to the large error circles for submm sources and the optical faintness of the small number of identified galaxies, existing redshifts confirm that the submm population lies at redshifts $`1<z<3`$. The existence of these sources implies a population of very luminous galaxies at very early epochs, with concommittant rapid evolution. The star-formation rates inferred from the luminosities of these systems are not easily produced in CDM galaxy formation models. In analogy with local LIGs, the submillimeter sources may be candidates for nascent QSOs at high-redshift.

## 2.  Active Galaxies with IRAS

One of the most striking results from IRAS is that the infrared galaxy population evolves more rapidly than the optical galaxy population (\[Lonsdale et al. 1990\]) with a form similar to that for QSOs. Selection of Luminous IR Galaxies associated with FIRST radio sources (\[Stanford et al. 2001\]) or by FIR/Optical ratio (\[Smith et al. 2002\]; \[Lonsdale et al. 2002\]) have led to the detection of IRAS Faint Sources up to $`z1`$ with a small number of “hyperluminous” or lensed sources with redshifts in excess of 1 (e.g. \[Rowan-Robinson et al. 1991\]).

IRAS showed that the infrared spectra of QSOs are essentially similar, with a broad ’infrared bump’ from about 2$`\mu `$m to 1 mm, with a significant fraction (10–50%) of the bolometric luminosity emitted at infrared wavelengths. This infrared emission is almost certainly thermal dust emission, but the location of this warm dust remains uncertain (\[Sanders et al. 1989\]) possibly associated with the outer edges of the accretion disk, the torus, or circumnuclear star formation.

Most of the known Sy 2 galaxies were discovered by IRAS, but the sensitivity of IRAS has limited limits our census of “Type 2 AGN” to $`z<0.2`$. As Padovani (1998) has previously stressed, at $`z<0.2`$ AGN Type 2 outnumber Type 1, suggesting that there are a large number of Sy 2 galaxies and Type 2 QSOs yet to be discovered at higher redshift. The expectation from our local census and from the CXB models suggests that an important task of future IR AGN Surveys will be detecting the high-redshift AGN 2 population, or explaining its absence in terms of AGN evolution.

## 3.  2MASS Quasars

The 2MASS AGN surveys have been reviewed by R. Cutri (2002; this volume) and the characteristics of that population will be discussed only briefly. The surface density of 2MASS AGN is of the order of 0.5 deg<sup>-2</sup> with a space density comparable to that found in low-redshift optical-to-x-ray selected samples. The colors and polarization of 2MASS QSOs suggests that these are reddened rather than intrinsically red QSOs and the ratio of Type 1 AGN to Type 2 is about 2:1 — the obscured QSO 2’s are not being found in 2MASS.

The color selection applied to detect 2MASS QSOs is $`(JK)>2`$ and Cutri et al. argue that these QSOs are sufficiently red to have been missed in previous optical surveys. Analysis by the SDSS and 2MASS teams (Ivezic et al. 2002; this volume) suggests, however, that the incompleteness of the SDSS to such red QSOs is less than 10%. The 2MASS AGN are faint in x-rays with a wide range in hardness-ratio (\[Wilkes et al. 2002\]) suggesting that the faintest and hardest sources may be highly obscured and could provide the missing x-ray population required to explain the CXB.

## 4.  ISO AGN studies and Surveys

### 4.1.  Mid-Infrared Diagnostics for AGN

ISO has shown that mid-infrared spectroscopy is a powerful tool for source classification. In a series of papers (\[Lutz et al. 1998\],\[Tran et al. 2001\]) the mid-infrared molecular band strengths haver been used to discriminate between AGN-dominated (7.7$`\mu `$m line to continuum $`<1`$) and Starburst-dominated (7.7$`\mu `$m to continuum $`>1`$) LIGs. These results suggest that the majority of LIGs with $`logL_{fir}<12(L_{})`$, the “ULIGs”, are Starburst-dominated, with an increasing AGN-fraction with increasing $`L_{fir}`$. Comparison with optical spectroscopy (\[Lutz, Veilleux & Genzel 1999\]) shows that the infrared classifications are consistent with classical optical excitation methods and suggests that infrared galaxies with LINER spectra are Starbursts (Veilleux 2002; this volume).

A somewhat different direction has been taken by Clavel et al. (2000) who compared the strengths of the mid-infrared bands among AGN types. In their sample there is a clear distinction between Sy 1 galaxies, showing low $`W_\lambda (7.7\mu m)`$, and Sy 1.5-2 galaxies with higher values of $`W_\lambda (7.7\mu m)`$. The 7.7$`\mu `$m luminosities of the Sy galaxy types are comparable, however. They interpret this result in terms of the classical AGN-torus model, suggesting that the warm AGN continuum from the inner torus is extinguished in the edge-on Sy 2 systems but visible in the Sy 1 systems, viewed face-on, whereas the extended, Starburst-related mid-infrared features are visible in both types of galaxies.

Taken together these results underscore the caveat discussed above, that AGN activity in compact, highly-obscured infrared galaxies may remain hidden even at mid-infrared wavelengths.

### 4.2.  ISO Deep Surveys

A number of deep surveys were undertaken with the ISO satellite at wavelengths of 7, 15, 90, and 170$`\mu `$m. These are reviewed in detail by Taniguchi (2002; this volume) and the results will only be summarized here. The principal result from these surveys is the continuing high surface densities of infrared galaxies as ISO pushed to lower flux-densities — as low as $`10\mu `$Jy at the lower wavelengths — requiring continued steep evolution with redshift. There are a number of models with varying prescriptions for the evolving population (\[Franceschini et al. 2001\]; \[Rowan-Robinson 2001\]; \[Xu et al. 2001\]). The models of Xu et al., based upon a local complete 24$`\mu `$m sample, require luminosity evolution rates as high as $`L(1+z)^{4.2}`$ combined with density evolution, $`\rho (1+z)^2`$, up to $`z1`$ for the Starburst population to match the $`logNlogS`$ relation; the number of LIGs at $`z1`$ is thus estimated to be approximately 40 times higher that in the local Universe (\[Franceschini et al. 2001\]). The Xu et al. model has pure luminosity evolution for the AGN population, $`L(1+z)^{3.5}`$, whereas the “normal” galaxy population evolves as $`L(1+z)^{1.5}`$.

With about 2000+ galaxies resolved by ISO and over 400 spectroscopic observations, AGN account for about 10% of the identified sources in the ISO Surveys. The ratio of AGN type 1 to type 2 in these surveys is about unity, at variance with the IRAS results at higher flux-density.

Elbaz et al. (2002) have analyzed the deep ISO 15$`\mu `$m counts and estimated that the IR galaxies detected by ISO to a 15$`\mu `$m flux-density, $`S_{15\mu m}>50\mu `$Jy, contribute over half of the CIRB and that AGN may contribute, at most about 20% of the IR background. Fadda et al. (2002) have combined the deep ISOCAM/XMM-Newton data from the HDFN and Lockman Hole with the brighter ELAIS S1/BeppoSAX data to estimate the AGN contribution th mid-IR surveys and the CIRB. Using x-ray emission as an indicator of AGN activity, Fadda et al. estimate that 15-20% of the mid-infrared emission in Lockman and HDFN originates from AGN. The detection (AGN) fraction of 15$`\mu `$m sources increases with x-ray energy from 30% below 2keV to over 60% above 5keV, as might be expected if the x-ray background is produced by obscured, high-column sources. Again, the fraction of the CIRB attributed to AGN-accretion energy is estimated to be less than 20%.

Further analysis of red x-ray luminous galaxies ($`L_x10^{43}`$$`10^{45}erg/s`$) in the same sample (\[Franceschini et al. 2002\]) suggests that these may be highly-obscured AGN as predicted by CXB background models. The mid-infrared SEDs of these sources are well reproduced by model spectra of obscured QSOs with $`\tau _{0.3\mu m}30`$–40 and the ratio of Type 1-to-Type 2 AGN is 1:3 in agreement with predictions.

The areas surveyed remain small, less than a few hundred square arcminutes to $`S_{15\mu m}<100\mu `$Jy, and statistics are restricted to small numbers with a few tens of confirmed AGN. The lower sensitivities at the longer wavelengths require the above analyses to employ template SEDs for estimating the contributions of ISO sources at the peak of CIRB near 140$`\mu `$m. With the wide range of mid-to-far-IR SEDS observed in the local Universe (e.g. the Far-IR/Mid-IR ratio may differ by an order of magnitude between a “typical Starburst” like M82 and a LIG such as Arp 220) these intriguing results remain tentative.

## 5.  SWIRE: The SIRTF Wide-area InfraRed Extragalactic Survey

For SIRTF, the last of its Great Observatories, NASA has selected a set of Legacy Programs, designed to be major surveys of general interest to the Astronomical community and to be carried out in the first year of the SIRTF mission. The Survey data will be distributed to the community in time for use in preparation of General Observer Proposals with no proprietary period for Legacy data.

The SIRTF Wide-area InfraRed Extragalactic Survey (SWIRE, Dr. Carol Lonsdale, P.I.) is the largest of the six SIRTF Legacy Surveys (851 hours), surveying approximately 67 square degrees in all 7 SIRTF imaging bands. A current description of the SWIRE Survey is given on the SWIRE WebPages: http://www.ipac.caltech.edu/SWIRE. Table 1 lists the Survey sensitivities.

The Survey will cover seven high-latitude fields, selected to be the most transparent, lowest background fields in the sky. The fields, covering between 5 and 15 sq. deg. include previously well-known IR extragalactic survey fields (e.g. Lockman and the ELAIS ISO Survey Fields) and x-ray fields (Chandra Deep South and XMM Large Scale Survey) are shown in Table 2.

The SWIRE science goal is to enable fundamental studies of galaxy evolution in the infrared for 0.5 $`<z<`$3:

* evolution of star-forming and passively evolving galaxies in the context of structure formation and environment.
* spatial distribution and clustering of evolved galaxies, Starbursts, & AGN.
* the evolutionary relationship between galaxies and AGN and the contribution of AGN accretion energy to the cosmic backgrounds.

Galaxy evolution models which match the IRAS/ISO galaxy counts at all wavelengths from 7–100$`\mu `$m as well as the CIRB (\[Xu et al. 2001\]) predict that SWIRE will detect of the order of 2 million galaxies — spheroids and evolved stellar systems with IRAC, and active star-forming systems with MIPS. SWIRE will also detect about 25,000 classical AGN, and an unknown number, perhaps several times as many, dust-enshrouded AGN.

Recent estimates of the “Universal Star-formation History (SFH)” (\[Steidel et al. 1999\]) suggest that the bulk of cosmic evolution occurs between redshifts, $`0.5<z<3`$, the redshift interval for which SWIRE is optimized. The median redshift is predicted to be, $`z1`$, where many estimates find a peak in the SFH; luminous infrared galaxies will be detected by SWIRE out to $`z3`$. Previous estimates of the SFH have varying, frequently large and uncertain corrections for extinction. SWIRE will directly measure the total star-formation rates as a function of redshift and environment over this critical range of time and redshift.

A key element in the SWIRE Survey design is to enable galaxy evolution studies in the context of large-scale structure/environment. One of the SWIRE Survey fields covers the deep survey areas of the XMM-LSS Survey (\[Pierre 2001\]) so that the infrared galaxy census may be directly tied to the presence of rich x-ray clusters to $`z>1`$. SWIRE will sample several hundred, 100 Mpc scale co-moving volume cells enabling a variety of large-scale structure measures from correlation functions, power spectra, and counts-in-cells to direct comparison with model calculations. SWIRE’s measures of the star-formation as a function of environment will be important input for CDM simulations which have been exceedingly successful in simulating the development of structure in the early Universe, but perhaps less so in simulating galaxy evolution within that structure owing to the complexity of the physics of star formation (e.g. Kay et al. 2002).

Of more direct importance to this Conference, the similarities of AGN SEDs in the mid-far Infrared suggests that SWIRE will be unbiased with respect to AGN types and ages, enabling a complete census of AGN out to redshifts greater than 1. Although the detection rates should be unbiased, the similarity between the SEDs of obscured AGN to those of Starbursts, and the extreme optical depths will make identifying the obscured AGN population very challenging. Low-frequency radio surveys will, of course, identify radio-loud AGN, but these make up only 10–15% of the AGN population. For this reason the XMM-LSS Survey, along with current and planned deeper surveys in hard x-rays will be vital to identifying SWIRE AGN.

### 5.1.  Supporting Observations

An aggressive program of ground-based optical, near-infrared and radio observations is planned in support of the SWIRE Survey and we are actively pursuing other programs with HST, Chandra, XMM and Galex. As already described Chandra and XMM Surveys will be important for discovering the obscured AGN population, if it exists. SWIRE has entered into cooperation with the Galex team so that the SWIRE fields will be included in the Galex Deep Survey.

The SWIRE Optical-Near Infrared goal is to obtain moderate-depth optical multi-band ($`g^{}25.7`$, $`r^{}25`$, $`i^{}24`$; Vega magnitudes, 5$`\sigma `$ detection for a $`2^{\prime \prime }`$ galaxy) data for the entire Survey area. At these limits we expect to detect approximately $`2/3`$ of SWIRE sources detected by both MIPS and IRAC. The ELAIS N1, N2 fields have already been imaged to somewhat shallower limits ($`r^{}24`$) as part of the INT Wide Field Survey and efforts continue to push deeper in the optical and into the near infrared as part of the UK SWIRE Program (ISLES and UKIDSS projects respectively). An extensive program for observations of ELAIS S1 is being undertaken at ESO. Optical and near-infrared imaging of the Lockman, Lonsdale and CDFS fields are being undertaken at KPNO and CTIO with the Mosaic cameras and FLAMINGOS infrared imager.

Two major SWIRE radio surveys are planned. The median 20cm flux density predicted for SWIRE Starburst galaxies is $`43\mu `$Jy — to faint to survey the entire area to this depth. We have therefore planned a deep pencil-beam VLA Survey (F. Owen, PI) and an extended shallow VLA survey (J. Condon, PI):

* SWIRE Lockman Deep Survey — 3$`\mu `$Jy rms @ 20cm; $`\alpha =10^h\mathrm{\hspace{0.17em}46}^m`$ $`\delta =+59^{}\mathrm{\hspace{0.17em}01}^{}`$; 30 VLA primary beam. The Deep VLA Survey is nearly completed and data analysis is just beginning.
* Cosmic Windows VLA Survey — 50$`\mu `$Jy rms @ 20cm in the combined fields of SWIRE, Galex and XMM-LSS which are accessible to the VLA. This Survey is being proposed for the next VLA large survey program.

## 6. Benediction

The SWIRE Legacy Survey is a community Survey; the large dataset which is being accumulated reflects the synergies which between the Legacy program and other community surveys. With a couple of million galaxies and several tens of thousands of AGN, many with redshift estimates and SEDs from x-ray to radio, the SWIRE database will be released to the community through IPAC’s InfraRed Science Archive. We hope that SWIRE will provide a rich datamine for the entire community and will provide answers to many of the questions posed here. If you have projects that you would like to do with SWIRE data, pleas visit the SWIRE WebPages and/or contact one of the team members.

It is a pleasure to thank the local organizers, especially Areg Mickaelian and Ed Khachikian, and the scientific organizing committee for my second stimulating visit to beautiful Armenia. This research was supported by the US NASA.


text:
# Cluster emission and phase transition behaviours in nuclear disassembly

## I Introduction

Phase transition and critical phenomenon is an extensively debatable subject in the natural sciences. Recently, the same concept was introduced into the astronomical objects  and the microscopic systems, such as in atomic cluster  and nuclei , of these the nuclei, as a microscopic finite-size system, are attracting more and more nuclear experimentalists to search for the liquid–gas phase transition and investigate their behaviour. To date, various experimental evidences have cumulated which seem to be related to the nuclear phase transition. For instance, violent heavy-ion collisions break the nuclei into several intermediate mass fragments, which can be viewed as a critical phenomenon as observed in fluid, atomic and other systems. It prompts a possible signature on the liquid–gas phase transition in the nuclear system. The sudden opening of the nuclear multifragmentation and vaporization  channels can be interpreted as the signature of the boundaries of phase mixture . In addition, the plateau of the nuclear caloric curve in a certain excitation energy range gives a possible indication of a first-order phase transition  as predicted in the framework of statistical equilibrium models . On the other hand, the extraction of critical exponents in the charge or mass distribution of the multifragmentation system  can be explained as an evidence of phase transition. More recently, the negative microcanonical heat capacity was experimentally observed in nuclear fragmentation  which relates to the liquid–gas phase transition , and in atomic cluster  which relates to solid to liquid phase transition , respectively. Moreover, some evidence of spinodal decomposition in nuclear multifragmentation was recently obtained experimentally , which shows the presence of liquid–gas phase coexistence region and gives a strong argument in favour of the existence of first-order liquid–gas phase transition in finite nuclear systems.

Meanwhile, several theoretical models have been developed to treat such a phase transition in the nuclear disassembly, e.g. percolation model, lattice gas model, statistical multifragmentation model and molecular dynamics model etc (e.g. see the recent review article of Richert and Wagner  and references therein). In this paper, we are interested in the lattice gas model ( LGM ) and classical molecular dynamics ( CMD ) model . The former is a simple short range interaction model , but it can successfully be applied to nuclear systems with isospin symmetry and asymmetry. LGM is carried assuming that the system is in a ’freeze-out’ density $`\rho _f`$ with thermal equilibrium at temperature T. Previous calculations  with LGM showed that there exists a phase transition for the finite nuclear systems by studying the effective power law parameter ( $`\tau `$ ) of cluster mass or charge distribution, their second moments ( $`S_2`$ ) and the specific heat. More recently, we proposed two novel criteria, namely multiplicity information entropy ( $`H`$ ) and nuclear Zipf’s law to diagnose the onset of liquid–gas phase transition in the framework of the isospin-dependent LGM (I-LGM) and isospin-dependent classical molecular dynamics (I-CMD) model .

In this paper, we show that the emission rate of clusters is a useful tool to diagnose the nuclear liquid–gas phase transition, while the largest fluctuation of cluster multiplicities is simultaneously revealed at the point of the phase transition by investigating the features of light particles (LP), charged particles (CP), intermediate mass fragments (IMF) and the largest fragment (MAX) of disassembling source in I-LGM and I-CMD frameworks.

The paper is organized as follows. Section 2 gives the descriptions of I-LGM and I-CMD. Results and discussions are shown in section 3 where the multiplicities of cluster emissions, their slopes and their fluctuations are investigated. The influence of the ’freeze-out’ density on cluster emission is also presented in the framework of I-LGM and I-CMD and the role of Coulomb interaction is studied by comparing the results of I-LGM and I-CMD at a given ’freeze-out’ density. Finally, the conclusion is given in section 4.

## II Description of models

### II.1  Isospin-dependent lattice gas model

The lattice gas model of Lee and Yang , where the grand canonical partition function of a gas with one type of atoms is mapped into the canonical ensemble of an Ising model for spin 1/2 particles, has successfully described the liquid–gas phase transition for the atomic system. The same model has already been applied to the microscopic nuclear system, e.g. see the papers of Pan and Das Gupta et al . In order to better understand the context of this study, the models are described below.

In the LGM, A ( = N + Z ) nucleons with an occupation number s which is defined as the ‘spin’ s = 1 ( -1 ) for a proton (neutron) or s = 0 for a vacancy, are placed in the L sites of a, three-dimensional cubic lattice. Each cubic lattice has a size 1.0/$`\rho _0`$ = 6.25 $`fm^3`$ and can, at most, be occupied by a single nucleon, where $`\rho _0`$ = 0.16 $`fm^3`$ is the normal nucleon density. Nucleons in the nearest neighbouring sites interact with an energy s i s j. The Hamiltonian of the system is written as

$$H=\underset{i=1}{\overset{A}{}}\frac{P_i^2}{2m}\underset{i<j}{}ϵ_{s_is_j}s_is_j.$$
(1)
The interaction constant $`ϵ_{s_is_j}`$ is related to the binding energy of the nuclei. In order to incorporate the isospin effect in the LGM, the short-range interaction constant $`ϵ_{s_is_j}`$ is chosen to be different between the nearest neighbouring like nucleons and unlike nucleons,

$`ϵ_{nn}`$ $`=`$ $`ϵ_{pp}=0.MeV,`$
$`ϵ_{pn}`$ $`=`$ $`5.33MeV,`$ (2)
which indicates the repulsion between the nearest-neighbouring like nucleons and attraction between the nearest-neighbouring unlike nucleons. This kind of isospin-dependent interaction incorporates, to a certain extent, the Pauli exclusion principle and effectively avoids producing unreasonable clusters, such as di-proton and di-neutron clusters etc. The disassembly of the system is calculated at an assumed ‘freeze-out’ density $`\rho _f`$ = (A/L)$`\rho _0`$ , beyond $`\rho _f`$ nucleons are too far apart to interact.

In this model, N + Z nucleons are put in L sites by Monte Carlo sampling using the canonical Metropolis algorithm . As pointed out in , however, one has to be careful treating the process of Metropolis sampling in order to satisfy the detailed balance principle and therefore warrant the correct equilibrium distribution in the final state. Describing in detail, in this paper, first an initial configuration with N + Z nucleons is established. Second, for each event, a sufficient number of ’spin’-exchange steps are tested, e.g. 20 000 steps to let the system generate states with a probability proportional to the Boltzmann probability distribution with the Metropolis algorithm. In each ’spin’-exchange step, a random trial change on the basis of the previous configuration is made. For instance, we choose a nucleon at random and attempt to exchange it with one of its neighbouring nucleons or vacancies regardless of the sign of its ’spin’ (Kawasaki-like spin-exchange dynamics ), then compute the change $`\mathrm{\Delta }E`$ in the energy of the system due to the trial change. If $`\mathrm{\Delta }`$E is less than or equal to zero, accept the new configuration and repeat the next ’spin’-exchange step. If $`\mathrm{\Delta }E`$ is positive, compute the ’transition probability’ W = $`e^{\mathrm{\Delta }E/T}`$ and compare it with a random number r in the interval . If r $``$ W , accept the new configuration, otherwise retain the previous configuration. Twenty thousand ’spin’-exchange steps are performed to ensure that we get that the equilibrium state (afterwards we will show the ’spin’-xchange step dependence of some observables). Third, once the nucleons have been placed stably on the cubic lattice after 20000 ’spin’-exchange steps for each event, their momenta are generated by Monte Carlo sampling of the Maxwell–Boltzmann distribution. Thus various observables based on phase space can be calculated in a straightforward fashion for each event. One important point of such Monte Carlo Metropolis computations is that the above ’spin’-exchange approach between the nearest neighbors, independently of their ’spin’, is evidenced to be satisfied by the detailed balance condition as noted in . In other words, this sampling method guarantees that the generated microscopic states form an equilibrium canonical ensemble.

Once this is done the LGM immediately gives the cluster distribution by using the rule that two nucleons are part of the same cluster if their relative kinetic energy is insufficient to overcome the attractive bond :

$$P_r^2/2\mu ϵ_{s_is_j}s_is_j<0.$$
(3)
This method has been proved to be similar to the so-called Coniglio–Klein prescription in condensed matter physics  and was shown to be valid in LGM.

### II.2 Isospin-dependent classical molecular dynamics model

Since the LGM is a model of the nearest-neighbouring interaction, a long-range Coulomb force is not amenable to lattice gas type calculation. Pan and Das Gupta  provide a prescription, based on simple physical reasoning, to decide if two nucleons, occupying neighbouring sites form part of the same cluster or not . They first try to map the LGM calculation to a molecular dynamics type prediction, both first carried out without any Coulomb interaction. If the calculations match quite faithfully, then they can study the effects of the Coulomb interaction by adding that to the molecular dynamics calculation. Here we adopt the same prescription to use the molecular dynamics and therefore investigate the Coulomb effect. The results and conclusions are now compared and checked between the LGM and the CMD.Obviously, here we do not perform any ab initio molecular dynamics calculation but only use it for simulating the nuclear disassembly, starting from a thermally equilibrated source which has been produced by the above I-LGM: i.e. the nucleons are initialized at their lattice sites with Metropolis sampling and have their initial momenta with Maxwell– Boltzmann sampling. From this starting point we switch the calculation to CMD evolution under the influence of a chosen force. Note that in this case $`\rho _f`$ is, strictly speaking, not a ‘freeze-out’ density for molecular dynamics calculation but merely defines the starting point for time evolution. However, since classical evolution of a many particle system is entirely deterministic and the initialization does have in it all the information of the asymptotic cluster distribution, we continue to call $`\rho _f`$ the ’freeze-out’ density.

The form of the force in the CMD is also chosen to be isospin dependent in order to compare with the results of I-LGM. The potential for unlike nucleons is expressed as

$`v_{\mathrm{np}}(r)(r/r_0<a)`$ $`=`$ $`C\left[B(r_0/r)^p(r_0/r)^q\right]exp({\displaystyle \frac{1}{(r/r_0)a}}),`$
$`v_{\mathrm{np}}(r)(r/r_0>a)`$ $`=`$ $`0.`$ (4)
where $`r_0`$ = 1.842 fm is the distance between the centres of two adjacent cubes so that $`\rho _0`$ = 1/$`r_{0}^{}{}_{}{}^{3}`$ = 0.16 $`fm^3`$. The parameters of the potentials are p = 2, q = 1, a = 1.3, B = 0.924 and C = 1966 MeV. With these parameters the potential is minimum at $`r_0`$ with the value -5.33 MeV, zero when the nucleons are more than 1.3$`r_0`$ apart and strongly repulsive when r is significantly less than $`r_0`$ . We now turn to the nuclear potential between like nucleons. Although we take $`ϵ_{pp}`$ = $`ϵ_{nn}`$ = 0 in I-LGM, the fact that we do not put two like nucleons in the same cube suggests that there is short-range repulsion between them. We have taken the nuclear force between two like nucleons to be the same expressions as above +5.33 MeV up to r = 1.842 fm and zero afterwards,

$`v_{\mathrm{pp}}(r)(r<r0)`$ $`=`$ $`v_{\mathrm{np}}(r)v_{\mathrm{np}}(r_0),`$
$`v_{\mathrm{pp}}(r)(r>r0)`$ $`=`$ $`0.`$ (5)
Figure 1 shows the above potential $`v_{np}`$ or $`v_{pp}`$. This potential form automatically cuts off at r/$`r_0`$ = a (equation (4)) or r/$`r_0`$ = 1 ( equation (5) ) without discontinuities in any r derivatives, which is a distinct advantage in any molecular dynamics simulation application.

The system evolves with the above potential. The time evolution equations for each nucleon are, as usual, given by

$`\stackrel{}{p_i}/t`$ $`=`$ $`\mathrm{\Sigma }_{ji}_iv(r_{ij}),`$
$`\stackrel{}{r_i}/t`$ $`=`$ $`\stackrel{}{p_i}/m.`$ (6)
Numerically, the particles are propagated in the phase space by a well-known Verlet algorithm , one of the finite-difference methods in molecular dynamics with continuous potentials. At asymptotic times, for instance, the original blob of matter expands to 64 times its volume in the initialization, the clusters are easily recognized: nucleons which stay together after an arbitrarily long time are part of the same cluster. The observables based on cluster distribution in both models are now compared while they are also compared by switching on/off the Coulomb interaction within the molecular dynamics.

## III  Results and discussions

We choose a medium-size nucleus $`{}_{}{}^{129}Xe`$ as an example. The input parameters are temperature T and ’freeze-out’ density $`\rho _f`$ in the model calculations. In this study T is mostly limited to the range of 1.5–7 MeV and the ’freeze-out’ density $`\rho _f`$ changes in a wide range, namely, 0.097$`\rho _0`$, 0.18$`\rho _0`$, 0.38$`\rho _0`$ and 0.60$`\rho _0`$, which corresponds to the total cubic lattices $`11\times 11\times 11`$, $`9\times 9\times 9`$, $`7\times 7\times 7`$ and $`6\times 6\times 6`$, respectively, of which, 0.60$`\rho _0`$ is the maxmium freeze-out density which is allowed for $`{}_{}{}^{129}Xe`$ because a cubic lattice is required in our LGM calculations. One thousand events are simulated for each combination of T and $`\rho _f`$ which ensures good statistics for results.

Before we present the results, we would like to check the role of the ’spin’-exchange step and from that we can know how the system tends to the equilibrium state in the model. Suppose, we wish to determine experimentally the value of a property of a system such as internal energy or cluster emission. In general, such properties depend upon the phase space of A nucleons. Over time, the instantaneous value of the property fluctuates as a result of interactions between the nucleons. However, when the system reaches its equilibrium state, the average of the instantaneous values over huge samples (‘ensemble average’) are viewed as experimental asymptotic values. But when will the system approach its equilibrium in the framework of LGM? To this end, we use the time average of the instantaneous values to investigate this question. In the LGM Monte Carlo simulation, the ’spin’-exchange step is viewed as a time step. We display some step-averaged obsevables that evolve with the step. Figure 2 shows the step-averaged total kinetic energy per nucleon $`E_{kin}`$, total potential energy per nucleon $`E_{pot}`$, the multiplicity information entropy $`H`$ which is defined in section 3.1., the multiplicities of emitted neutrons, protons, CP, IMF, and the mass of the heaviest fragment at $`\rho _f`$ = 0.38$`\rho _0`$ and T= 3 MeV or 6 MeV in I-LGM. From this figure we see that the observables tend to their asymptotic values after 1000 steps, i.e. the system tends to equilibrium after 1000 ’spin’-exchange steps in the lattice gas Monte Carlo simulation. In the following results, we adopt the results after 20 000 ’spin’-exchange steps when the system is in the equilibrium state.

### III.1 I-LGM in different ’freeze-out’ densities

Figure 3 shows that the mean multiplicities of emitted neutrons, protons, CP and IMF and the mean mass for the largest fragment evolve with temperature at different ‘freeze-out’ densities in the I-LGM calculation. Here IMF is defined as 3 $``$ Z $``$ 20. At a fixed ‘freeze-out’ density, average neutron multiplicity ( $`N_n`$ ), proton multiplicity ( $`N_p`$ ), charged particle multiplicity ($`N_{cp}`$ ) and the mean mass for the largest fragment ( $`A_{max}`$ ) display monotonous increase or decrease with temperature as expected. But the multiplicity ( $`N_{imf}`$ )ofIMFshowsariseandfallwith temperature , when the system probably crosses the phase transition boundary. With decreasing ‘freeze-out’ density, $`N_p`$, $`N_n`$, $`N_{cp}`$ and $`A_{max}`$ increase, since larger space separation among nucleons at smaller ‘freeze-out’ density makes the clusters less bound and therefore the sizes of free clusters decrease and then the cluster multiplicities increase. The situation of $`N_{imf}`$ is slightly complicated, i.e. it increases with the decrease of ‘freeze-out’ density in the lower temperature branch contrary to the higher temperature branch.

It seems difficult to discover the possibility of phase transition of nuclei if we only see these mean quantities as shown above ( $`N_{imf}`$ is an exception ). However, when we focus on their slopes to temperature (figure 4), sharp changes are observed at nearly the same temperature at each fixed ‘freeze-out’ density, for instance, 3.5 MeV at 0.097$`\rho _0`$, 4MeV at 0.18$`\rho _0`$, 5MeV at 0.38$`\rho _0`$ and 6 MeV at 0.60$`\rho _0`$. At such a transition point, (1) the multiplicities of emitted clusters increase rapidly and after that the emission rate slows down; and (2) the decrease in the largest fragment size reaches a valley for such a finite system. Physically, the largest fragment is simply related to the order parameter $`\rho _l\rho _g`$ (the difference of density in nuclear ‘liquid’ and ‘gas’ phases). In infinite matter, the infinite cluster exists only on the ‘liquid’ side of the critical point. In finite matter, the largest cluster is present on both sides of the phase transition point. In this calculation, a valley for the slope of $`A_{max}`$ to temperature may correspond to a sudden disappearance of infinite cluster (‘bulk liquid’) near the phase transition temperature. It is not the occasional production of such waves of the slopes; it should reflect the onset of phase transition there. This idea is supported by surveying the other phase transition observables, such as the effective power law parameter $`\rho _f`$rom the mass or charge distribution of fragment and the information entropy $`H`$ of event multiplicity distribution . $`H`$ can be expressed as

$$H=\underset{iEvent}{}p_iln(p_i)$$
(7)
where $`p_i`$ is defined as the event-normalized total multiplicity probability and $`_{iEvent}p_i`$ = 1. Figure 5 depicts these results. The minima of $`\tau `$ and the maxima of $`H`$ appear around respective phase transition temperatures at different values of ‘freeze-out’ density, i.e. about 3.5 MeV at 0.097$`\rho _0`$ , 4–4.25 MeV at 0.18$`\rho _0`$, 5–5.5 MeV at 0.38$`\rho _0`$ and 6–6.5 MeV at 0.60$`\rho _0`$. These temperatures are consistent with those extracted from the choppy position of the above slopes. It indicates that the above slopes (emission rate) can be taken as a probe of a liquid–gas phase transition of nuclei.

Furthermore, the largest fluctuations of cluster multiplicities are found around the phase transition point in the same calculation. Figure 6 illustrates that RMS width ( $`\sigma `$ ) of the multiplicity distributions of neutrons, protons, CP and IMF, and the distribution of the largest fragment masses. These variances generally show peaks at the same phase transition temperatures as those extracted from the above observables for each fixed ‘freeze-out’ density. Note that the fluctuation of $`A_{max}`$ is related to the compressibility of the system. These features are also consistent with one of the phase transition behaviours, i.e. the largest fluctuation at the phase transition point exists . This fluctuation represents an internal feature of the disassembling system, not a numerical fluctuation.

Another way to characterize the fluctuation is the method of the conditional moments introduced by Campi . The normalized second moment $`S_2`$ in each event is defined as

$$S_2=\frac{_{A_iA_{max}}A_{i}^{}{}_{}{}^{2}n_i(A_i)}{_{A_iA_{max}}A_in_i(A_i)}$$
(8)
where $`n_i`$ is the multiplicity of cluster mass $`A_i`$ , and the summation is over all clusters in an event except the heaviest one which corresponds to the bulk liquid in an infinite system. For example, we show the 1000-event Campi scatter plots ( ln($`A_{max}`$) versus ln($`S_2`$) ) as a function of temperature at a fixed density of 0.38$`\rho _0`$ in figure 7. Actually, the plots clearly illustrate the evolution of the disassembling mechanism with temperature. At lower temperatures, only the under-critical (liquid phase) branch with a negative slope of ln($`A_{max}`$) versus ln($`S_2`$) exists while at higher temperature, only the super-critical (gas phase) branch with a positive slope of ln($`A_{max}`$) versus ln($`S_2`$) appears. However, both the branches (the liquid–gas phase coexistence region) meet closely around 5.5 MeV, which indicates the onset of the liquid–gas phase transition . Similar behaviour shows for the calculation at 0.097$`\rho _0`$ ,0.18$`\rho _0`$ and 0.60$`\rho _0`$ at their respective phase transition temperature.

We point out that the ‘freeze-out’ density-dependent phase transition temperature, extracting from $`\tau `$ and $`S_2`$ , was also observed in a previous study . Similar to water, the temperature of liquid–gas phase transition decreases with pressure. In the nuclear case, the decrease in ‘freeze-out’ density is similar to the decrease in the internal pressure inside nuclei, hence it leads to decrease in transition temperature. But this $`\rho _f`$-dependent phenomenon vanishes when excitation energy is used as a variable in the LGM . In other words, the excitation energy has perhaps a good correspondence with the critical temperature because of only one critical point and hence only one critical temperature for a system.

### III.2 Roles of Coulomb force: comparison of I-LGM and I-CMD in fixed ’freeze-out’ density

Considering the absence of long-range Coulomb force in the LGM, we adopt the CMD to investigate the Coulomb effect and check the features of cluster emission and its relation to phase transition behaviour. To this end, first we make a comparison for the results of I-CMD with Coulomb or without Coulomb and those of I-LGM at a certain fixed ‘freeze-out’ density, namely 0.38$`\rho _0`$ . Second, we present the results of all these observables in the frame of I-CMD with Coulomb at four different ‘freeze-out’ densities to check cluster emission and its relation to the phase transition behaviour in next subsection.

Figure 8 shows that $`N_p`$ , $`N_n`$ , $`N_{cp}`$ , $`N_{imf}`$ and $`A_{max}`$ change with temperature in different calculation cases, i.e. I-CMD with Coulomb or without Coulomb and I-LGM (see the meaning of the symbols in the figure.). The multiplicities of clusters and the largest fragment mass are close to each other between I-LGM and I-CMD with Coulomb except for the multiplicities of neutrons and protons, illustrating that the I-LGM is, in general, a good tool to describe the fragmentation if Coulomb interaction can be ignored. When Coulomb interaction is switched on, $`N_p`$ , $`N_{cp}`$ and $`N_{imf}`$ increase due to the repulsive role among protons while $`A_{max}`$ decreases. Meanwhile, $`N_n`$ does not change because of no Coulomb interaction. $`N_{imf}`$ also shows a rise and fall with temperature in the I-CMD cases. Coulomb force makes the turning temperature of $`N_{imf}`$ smaller due to long-range repulsion.

The slopes of multiplicities of emitted clusters and of mean mass of the largest fragment are plotted as a function of temperature in case of I-CMD in figure 9. The definite peaks of slopes are found as in the I-LGM case. The corresponding temperature at the peaks is located about 4 MeV in the I-CMD case with the Coulomb interaction and about 5 MeV in the I-CMD case without the Coulomb interaction. This turning temperature also reflects the onset of phase transition there. If we investigate $`\tau `$ and $`H`$ (figure 10), we find that there are a minima of $`\tau `$ and the maxima of $`H`$ around phase transition temperatures, i.e. about 4–4.25 MeV for I-CMD with the Coulomb, around 5 MeV for I-CMD without the Coulomb and around 5.5 MeV for I-LGM respectively.

Finally, the RMS widths of the multiplicity distributions of the clusters and of the largest fragment mass are checked in the I-CMD case in figure 11. The widths of the multiplicity distributions of neutrons and protons tend to be saturated at higher temperature, while those for CP, IMF and $`A_{max}`$ demonstrate peaks at a certain fixed temperature, i.e. around 4 MeV for the case I-CMD with the Coulomb, 5 MeV for the case of I-CMD without the Coulomb, which is similar to the I-LGM case. These turning temperatures are also consistent with the phase transition temperature as shown in figures 8–10 in the I-CMD cases.

Overall, Coulomb interaction plays a notable role in favour of fragment production and reduces the temperature of the phase transition. When the system is small, the Coulomb interaction is not expected to be important and I-LGM could be a good tool to treat nuclear disassembly. But for large nuclear systems, the neglect of the Coulomb force in LGM is rather serious handicap. In this background, first, we made a suitable selection for the molecular dynamics interaction potential and then got good agreement between the I-LGM and I-CMD without the Coulomb. Thanks to this agreement, we treated the nuclear disassembly more realistically with switching on of the Coulomb interaction in the frame of I-CMD afterwards.

### III.3 I-CMD in different ’freeze-out’ densities

Similar to the I-LGM cases, we also study the cluster emission in I-CMD in a wide range of ‘freeze-out’ density, namely, 0.097$`\rho _0`$ ,0.18$`\rho _0`$ ,0.38$`\rho _0`$ and 0.60$`\rho _0`$ to check phase transition behaviour.

In order to compare it with the I-LGM case, similar figures are plotted in figures 12–15 and compared to figures 3–6. Figure 12 shows that the mean multiplicities of emitted neutrons, protons, CP, IMF and $`A_{max}`$ evolve with temperature at different ‘freeze-out’ densities in the I-CMD calculation with the Coulomb. Their emission rates or slopes and RMS width with temperature are depicted in figures 13 and 15. The effective power-law parameter of fragment distribution $`\tau `$ and multiplicity information entropy $`H`$ is plotted in figure 14. All these figures show the same behaviours as the I-LGM case, i.e. cluster emission rate and their fluctuation can provide us with the temperature of liquid–gas phase transition, namely around 3, 3.5, 4.5, 5 MeV at $`\rho _f`$ = 0.097$`\rho _0`$ ,0.18$`\rho _0`$ ,0.38$`\rho _0`$ and 0.60$`\rho _0`$ respectively.

Similarly, the Campi’s scattering plots indicate the onset of liquid–gas phase transition around the respective transition temperatures. Figure 16 shows an example for I-CMD with the Coulomb. Clearly, phase coexistence takes place around 4.25 MeV. Overall, these phase transition temperatures are consistent with those extracted from the choppy position of the slopes of figure 13. Again, the largest fluctuation simultaneously appears at the point of phase transition.

Overall, the phase transition temperature seems to rely on some ingredients, such as the ‘freeze-out’ density (or pressure), the model and its interaction potential, but the rule that emission rate and fluctuation of cluster multiplicity can be taken as a probe of phase transition has not changed.

## IV Conclusions

In conclusion, the features of the emissions of LP, CP, IMF and MAX are investigated in a wide range of ‘freeze-out’ density for a medium size nucleus $`{}_{}{}^{129}Xe`$ in the frameworks of I-LGM and I-CMD model. $`N_n`$, $`N_p`$, $`N_{cp}`$ and $`A_{max}`$ show monotonous increase or decrease while $`N_{imf}`$ shows rise and fall with temperature. Slopes of these observables versus temperature go through extrema at the same temperature where the largest fluctuation of cluster multiplicity distributions is observed. This temperature is consistent with the phase transition temperature extracted from the extreme values of effective power law parameter $`\tau `$ and information entropy $`H`$ as well as Campi scatter plots. It gives an indication that the cluster emission rate can be taken as a probe of the phase transition of nuclei and furthermore, the largest fluctuation is simultaneously accompanied when the onset of phase transition occurs. In addition, the systematic comparison of I-LGM and I-CMD shows that LGM is a good tool to study nuclear disassembly when the system is not large where the Coulomb interaction can be ignored. But for large nuclear systems, I-CMD should be used to treat the nuclear dissociation and phase transition due to the importance of the Coulmb interaction. In light of this study, we think that the experimental study of cluster emission is rather meaningful, especially in measuring the excitation function of the multiplicities, their slopes and variances, from which some signals of phase transition could be found.

###### Acknowledgements.

I thank Prof S Das Gupta and Dr J Pan for kindly providing the original codes. I also appreciate Prof B Tamain and Prof J B Natowitz for their help. This study was partly supported by the National Natural Science Foundation of China for the Distinguished Young Scholar under Grant no 19725521, the National Natural Science Foundation of China under Grant no 19705012, the Science and Technology Development Foundation of Shanghai under grant no 97QA14038 and the Major State Basic Research Development Program of China under contract no G200077400.