text:
# Estimates of the total gravitation radiation in the head-on black hole collision

## 1 Introduction

Accurate description of the dynamics of relativistic gravitating systems remains one of the challenging problems in general relativity. There is great interest in particular on the description of systems involving very compact objects, as for example black holes. It is expected that these kind of systems will be describable in terms of an appropriate isolated model; which in a relativistic theory of gravity is expressed by an asymptotically flat spacetime. In this description the gravitational radiated energy will be encoded in the behavior of the physical field at future null infinity.

In numerical calculations these ideas have been used in attacking the problem of the time symmetric asymptotically flat initial data corresponding to a two equal mass black hole system. It would be nice to be able to compare these calculations with some analytic solution of the field equation; but in trying to do this one is faced with the fact that there is up to now only one known vacuum asymptotically flat radiating family of solutions, the so called Robinson-Trautman metrics. These spacetimes represent the dynamics of a single black hole which settles down, in the asymptotic future, to the Schwarzschild geometry.

In the process of two black holes falling into each other in a head-on collision, one can distinguish two natural eras: one in which the two black holes are separated and falling, and the other in which they already form a single black hole which is decaying to a stationary state. Since the Robinson-Trautman metrics represent a black hole, one could estimate the total gravitational radiation emitted after the formation of the black hole using these geometries. In doing this one needs to associate the data of the two black hole collision system to the Robinson-Trautman data. We do this by matching, at an arbitrary retarded time, the radiation content of these metrics with the radiation calculated from the quadruple formula with Newtonian dynamics at the moment of collapse.

## 2 Power of the energy radiated from Newtonian dynamics + quadrupole formula

Let us consider the Newtonian problem of two particles of mass $`m_0`$, which start at rest at an initial distance $`d_i`$. One can estimate the power of the energy carried away by gravitational waves when they are at a distance d assuming Newtonian dynamics and the quadrupole approximation. More specifically; the quadrupole radiation formula relates the power of the energy radiated with the time derivatives of the quadrupole:

$$P_Q=\frac{1}{5}\underset{i,j=1}{\overset{3}{}}\left(\frac{d^3Q_{ij}}{dt^3}\right)^2$$

where

$$Q_{ij}=q_{ij}\frac{1}{3}\delta _{ij}q,$$

$$q_{ij}=T^{00}x^ix^jd^3x,$$

and we are using geometric units in which the Newton’s gravitational constant and the speed of light have the unit value. Then, the total power radiated is

$$P_Q(d)=\left(\frac{1}{15}\right)\frac{1}{\left(\frac{d}{M_0}\right)^4}\left(\frac{1}{\left(\frac{d}{M_0}\right)}\frac{1}{\left(\frac{d_i}{M_0}\right)}\right),$$

where $`M_02m_0`$ is the total mechanical mass and d the distance at observation.

## 3 Robinson-Trautman metrics: the quadrupole approximation

In a previous paper it was found the general asymptotic behavior of the Robinson-Trautman metrics in the asymptotic future. The leading behavior of these metrics is governed by the quadrupole structure of the spacetime.

More explicitly, the whole geometry of these metrics is determined by a scalar $`V(u,\theta ,\varphi )`$ , where u, can be considered a retarded time and $`(\theta ,\varphi )`$ the coordinates of a 2-sphere. According to Ref. , for an axisymmetric quadrupole structure, the asymptotic behavior of this scalar, for $`u\mathrm{}`$ is

$$V=1+Ae^{\frac{2u}{M_{\mathrm{}}}}Y_{20}+O(e^{\frac{4u}{M_{\mathrm{}}}}),$$

where A is a constant $`M_{\mathrm{}}`$ is the asymptotic mass of the spacetime in the regime $`u\mathrm{}`$, and $`Y_{20}`$ is a spherical harmonic function. This geometry undoubtedly represents the dynamics of a single black hole; which becomes evident when in the limit for $`u\mathrm{}`$ , i.e., V=1, one recognizes the Schwarzschild metric.

As we have described above, there are two natural eras in the process of the head-on collision: before and after the formation of the remaining black hole. Since the Robinson-Trautman geometries describe the dynamics of a single black hole, it is natural to try to study, with these geometries, the final stage of the process.

The time derivative of the Bondi shear of the Robinson-Trautman sections of future null infinity has the following asymptotic behavior:

$$\dot{\sigma }_0=\sqrt{6}Ae_2^{\frac{2u}{M_{\mathrm{}}}}Y_{20}+O(>),$$

where now $`{}_{2}{}^{}Y_{20}^{}`$ is a spin weight 2 spherical harmonic function and $`O(>)`$ means higher order terms of the time exponential. With this information one can calculate the flux of gravitational Bondi radiation at future null infinity in first order, obtaining:

$$P_B=\frac{6}{4\pi }A^2e^{\frac{4u}{M_{\mathrm{}}}}.$$

To relate the Robinson-Trautman geometry with the problem of the head-on collision of two black holes, we identify the Bondi power at some retarded time with the power calculated from the quadrupole formula at the moment when the two black holes touch; that is when $`d=2M_0`$ . In this way we determine A appearing in the above formulae.

With the Bondi flux one can calculate in first order the total energy radiated, namely:

$`E={\displaystyle _{u_0}^{\mathrm{}}}P_B(u)𝑑u`$ $`=`$ $`{\displaystyle \frac{6}{4\pi }}A^2{\displaystyle \frac{M_{\mathrm{}}}{4}}e^{4u/M_{\mathrm{}}}|_{\mathrm{}}^{u_0}`$
$`=`$ $`{\displaystyle \frac{6}{4\pi }}A^2{\displaystyle \frac{M_{\mathrm{}}}{4}}e^{4u_0/M_{\mathrm{}}}=P_B(u_0){\displaystyle \frac{M_{\mathrm{}}}{4}}.`$
Let $`M_i`$ denote de mass of the spacetime at the *initial* retarded time; then one has the relation

$$M_{\mathrm{}}=M_iE.$$

Therefore, the ratio radiation over initial mass, can be expressed by

$$\frac{E}{M_i}=\frac{M_{\mathrm{}}\frac{P_B(u_0)}{4}}{M_{\mathrm{}}\left(1+\frac{P_B(u_0)}{4}\right)}=\frac{\frac{P_B(u_0)}{4}}{\left(1+\frac{P_B(u_0)}{4}\right)}.$$

Identifying the Bondi with the quadrupole flux one would obtain

$$\frac{E}{M}=\frac{\left(\frac{1}{60}\right)\frac{1}{l_f^4}\left(\frac{1}{l_f}\frac{1}{l}\right)}{1+\left(\frac{1}{60}\right)\frac{1}{l_f^4}\left(\frac{1}{l_f}\frac{1}{l}\right)},$$

with $`l_f\frac{d_f}{M}`$ and $`l\frac{d_i}{M}`$ . In this expression we identify M with $`M_i`$ ; that is we map the Newtonian total mass to the total mass of the spacetime at the initial retarded time.

## 4 Time symmetric 2-black hole initial data

In reference  it was presented the initial data corresponding to a time symmetric 2-black hole system. There the ADM mass $`M_M`$ was calculated, and a magnitude $`D_M`$ (there denoted with capital l) was defined, having the interpretation of some measure of initial distance between the black holes. More precisely $`D_M`$ was defined as the length of the minimal line from the “ring” at minimum radius from one throat to the respective “ring” in the other throat.

If one wants to compare calculations coming from different realizations of the same physical situation, namely the head-on collision of two equal mass black holes, one must provide a relation between the physical quantities involved. The two fundamental physical quantities involved in this situation are the total mass of the system and the separation of the black holes. As we have done in the previous section, the identification of the masses is straight forward, since the natural thing to do is to identify the invariantly defined total masses in each case. There remains however the relation among the notion of distance involved in each case.

The question is: Is there a natural way to identify the Misner magnitude $`D_M`$ with the Newtonian notion of distance appearing in the estimates of the total gravitational energy radiated using the quadrupole formula approach? A simple identification of $`D_M`$ with d would be too naive in the regime of small distances. We seek then for some relation between $`D_M`$ and d such that it captures the most important physical aspects of both systems.

In order to do this let us recall that in reference  the time symmetric initial data for a two equal mass black hole system was represented by the 3-geometry:

$$ds^2=\left(1+\frac{m_1\mathrm{}}{2r_1}+\frac{m_2\mathrm{}}{2r_2}\right)^4ds_F^2,$$

where $`ds_F^2`$ is the line element of a flat 3-geometry, $`r_1`$ and $`r_2`$ are the Euclidean distances from the field point to the deleted points in the flat 3-geometry which represent the location of the first and second black hole respectively, and $`m_1\mathrm{}`$ and $`m_2\mathrm{}`$ are parameters which characterize the spacetime. If in the conformal, scalar appearing in the above line element, one eliminates the third term, the constant time Schwarzschild line element is obtained, with mass $`m_1\mathrm{}`$ expressed in terms of isotropic coordinates.

According to the arguments of reference  the mass of, let us say, the first black hole is

$$m_1=m_1\mathrm{}+\frac{m_1\mathrm{}m_2\mathrm{}}{2r_{12}},$$

where now $`r_{12}`$ is the Euclidean distance between the black holes; while the total mass of the spacetime is

$$M=m_1+m_2\frac{m_1\mathrm{}m_2\mathrm{}}{r_{12}}.$$

This expression agrees with the first relativistic correction to the Newtonian total mass; since in first order one has

$$M_{rel}M_0+E_0=m_0+m_0\frac{m_0^2}{d_i}.$$

Therefore it seems that the Newtonian distance, for large $`d`$, should be identified with the coordinate distance of the flat 3-geometry of the time symmetric problem. So using the line element $`ds^2`$ shown above for the region between the black holes, and identifying the radius of each black hole with double their respective masses, we define the distance between the black holes by

$$D=2m_1+\underset{2m_1}{\overset{d2m_2}{}}\left(1+\frac{m_1\mathrm{}}{2r}+\frac{m_2\mathrm{}}{2(dr)}\right)^2𝑑r+2m_2,$$

which in terms of the total mass M for a two equal mass system can be expressed by

$$D=2M\left(1+\frac{M}{4d}\right)+\underset{M\left(1+\frac{M}{4d}\right)}{\overset{dM\left(1+\frac{M}{4d}\right)}{}}\left(1+\frac{m_{\mathrm{}}}{2r}+\frac{m_{\mathrm{}}}{2(dr)}\right)^2𝑑r;$$

then, defining $`L\frac{D}{M}`$ and $`l\frac{d}{M}`$ one has the relation

$$L=2\left(1+\frac{\mathrm{\hspace{0.33em}1}}{4l}\right)+\underset{\left(1+\frac{\mathrm{\hspace{0.33em}1}}{4l}\right)}{\overset{l\left(1+\frac{\mathrm{\hspace{0.33em}1}}{4l}\right)}{}}\left(1+\frac{1}{4\lambda }+\frac{1}{4(l\lambda )}\right)^2𝑑\lambda .$$

It is observed that with respect of this variable, the black holes do not touch exactly at $`l=2`$ but when $`l^{}=2\left(1+\frac{1}{4l^{}}\right)`$; which corresponds to $`L=l^{}=2.22`$. This picture is in agreement with the numerical calculations since it is found in reference  that for $`L<2.49`$ a common apparent horizon surrounds both holes 111The threshold mentioned in their reference is $`\mu <1.36`$ which corresponds to our $`L<2.49`$..

We relate d to the Misner magnitude $`D_M`$ by setting $`D(d)=D_M`$.

## 5 Upper limit for the radiation from the area theorem

One can estimate an upper limit for the total energy radiated away in the form of gravitational waves when two very separated black holes collide, from the area theorem; which using the relation between total mass and single black hole mass shown above and in terms of the quantities already defined, establishes that

$$\frac{E_{rad}}{M}\left(1\frac{1}{\sqrt{2}}\right)\left(1+\frac{\mathrm{\hspace{0.33em}1}}{4l}\right).$$

## 6 Comparison of estimates with supercomputer results

In ref.  the numerical calculation of the total energy radiated for the time symmetric head-on collision of two equal mass black holes was reported. They improved previous calculations  with the help of latest generation of supercomputers, and further developed analytic and numerical techniques. In the next graph we plot their recent data, our Robinson-Trautman (RT) estimate for the total energy radiated and the upper limit from the area theorem.

It can be observed that the RT estimate keeps bellow the supercomputer data; but this is no surprise since one had expected to only estimate, by this mean, the total energy radiated after the formation of the resulting black hole.

The energy radiated in the first stage, when the black holes fall from $`d_i`$ to $`d`$, before the formation of the resulting black hole, can be estimated from the quadrupole formula and using the Newtonian dynamics; that is

$$E_Q=_{t_{d_i}}^{t_d}P_Q\left(x(t)\right)𝑑t=_{d_i}^dP_Q(x)\frac{1}{\dot{x}}𝑑x;$$

from which it is deduced that

$$\frac{E_Q}{M}=\left(\frac{1}{15\sqrt{2}}\right)\left(\frac{2}{7}\left(\frac{1}{l_f}\frac{1}{l}\right)^{\frac{7}{2}}+\frac{4}{5l}\left(\frac{1}{l_f}\frac{1}{l}\right)^{\frac{5}{2}}+\frac{2}{3l^2}\left(\frac{1}{l_f}\frac{1}{l}\right)^{\frac{3}{2}}\right);$$

with $`l\frac{d_i}{M}`$, $`l_f\frac{d_f}{M}`$ and M is the total mass of the system.

Then adding to the RT estimate the quadrupole Newtonian estimate one obtains the total energy radiated during the entire process in the quadrupole approximation, which is shown in the next graph.

This remarkable agreement between our estimates and the supercomputer calculations allows one to think that the global quantities, as the total energy radiated, for a gravitating system like the one under consideration, is not very dependent on the model realization of the system; since, in particular, the Robinson-Trautman approach is not time symmetric. Therefore the supercomputer results obtained from the time symmetric two black holes head-on collision system might have a wider range of validity than the one assumed up to now.

It is important to note that our estimates can be calculated for values of L corresponding to distance over mass relation above the threshold of the formation of the remaining black hole; which from reference  it is deduced to be around L=2.49. So our estimates are complementary of those calculated in reference  in which, using a completely different approach, they have estimated with high precision the total energy radiated for small L; while their method overestimates the total energy in two orders of magnitude in the vicinity of $`L8`$ (which correspond to their $`\mu _03`$ ).

The total gravitational energy radiated in the large separation regime was represented through semianalytic method in reference  by their equation (6). The approach followed in that reference requires of a series of extrapolation in the applicability of known equations. First of all they use an equation for the total energy radiated calculated in the test particle limit (eq. (1)); that is, one mass much smaller that the other mass ( $`\stackrel{~}{m}<<\stackrel{~}{M}`$ ) falling from infinity, and then they extrapolate the equation to the equal mass case ( $`\stackrel{~}{m}=\stackrel{~}{M}`$ ). This approach needs from several correcting factors. Using an altered quadrupole radiation formula they suggest a correcting factor, appearing in their eq. (4), that takes into account the effects associated to the fact that the infall is not from infinity but from a finite distance. The alteration in the quadrupole formula comes from a substitution of the radial velocity by a “nonunique” 222See reference  of reference . expression. It is observed that the integral in their eq. (4) is calculated up to the upper limit $`2\stackrel{~}{M}`$ , which is consistent with the test particle limit approach, but it overestimates the total gravitational energy radiation in the equal mass case. Therefore they have to introduce another correcting factor smaller than one, there called $`F_h`$, which is intended to take into account the black hole structure of both objects. These points were also studied in references  and . It is important to emphasize that we use, without alterations, the quadrupole radiation formula for the first stage before the collision, up to the point in which the two horizons touch each other ( $`d4\stackrel{~}{M}`$ ); this in contrast to the case of a test particle crossing one horizon ( $`d2\stackrel{~}{M}`$ ). Then, in order to take into account the black hole structure of the remaining object we use the Robinson-Trautman geometry, which radiation is to be added to the quadrupole estimate; this in contrast to the reducing factor $`F_h`$.

The present estimate of the total gravitational energy radiated is based on a couple of simple ideas with direct physical interpretation; which is a contribution to the understanding of the global aspects of the head-on black hole collision.

It would be interesting to reproduce these calculations taken into account post-Newtonian techniques, as those used in reference , to see whether these effects will make any substantial difference to our calculations.

A question that motivated us for this work was how a global quantity, as the total energy radiated, could be dependent on the details of the modeling of a physical situation. On purpose we use a completely different approach, and the result is that this global quantity does not appear to be very dependent on the details of the modeling.

As a final comment let us point out that since the time symmetric initial data obviously incorporates incoming radiation in the system, it is important to be able to know whether this incoming radiation will have an important contribution to the outgoing radiation. Since the Robinson-Trautman geometry can be considered to represent pure outgoing radiating spacetimes, and also the quadrupole estimate only takes into account retarded fields, the agreement of our estimates with the numerical calculations, of the time symmetric problem, suggests that the incoming radiation contribution to the outgoing one is negligible.

## Acknowledgements

It is a pleasure to thank Richard Price for a series of illuminating discussions, for pointing out to us reference  and for a detailed exposition of their work. We are also very grateful to Peter Anninos and collaborators for sending us their recent calculations.


format: text

text:
# Optical feedback and the coupling problem in semiconductor microdisk lasers

(Current address:
Department of Physics, University of Oregon, 1371 E 13th Avenue, Eugene, OR 97403
Published in physica status solidi (a)188, 921 (2001))

## Abstract

The smaller the size of a light-emitting microcavity, the more important it becomes to understand the effects of the cavity boundary on the optical mode profile. Conventional methods of laser physics, such as the paraxial approximation, become inapplicable in many of the more exotic cavity designs to be discussed here. Cavities in the shape of microdisks, pillars and rings can yield low lasing thresholds in a wide variety of gain media: quantum wells, wires and even dots, as well as quantum cascade superlattices and GaN. An overview of the experimental and theoretical status is provided, with special emphasis on the light extraction problem.

Light emission from microcavities is a problem of great fundamental and applied interest. A wide range of possible active media can be used to form microcavities, and consequently the properties of the microscopic dipoles which generate the light can vary significantly. One of our goals is to identify common, material-independent features of microcavity emitters, which are strongly determined by the geometric dimensions and shape of the cavity. Properties of interest are the spectrum of optical modes, their internal intensity distribution and the external field profile of the emitter. These characteristics in turn determine technological figures of merit, such as external quantum efficiency in the case of light emitting diodes (LEDs), or pump threshold and maximum output power in the case of lasers.

A problem in the design of LEDs is their poor external quantum efficiency, owing to the fact that light generated within the diode is not easily extracted. This is because total internal reflection at the interface between the semiconductor ($`n_{in}`$) and the surrounding lower-index medium ($`n_{out}`$) allows only those light rays to escape whose angle of incidence $`\chi `$ with respect to the surface normal satisfies

$$\mathrm{sin}\chi <n_{out}/n_{in}1/n.$$
(1)
A first step toward better output coupling in LEDs is to reduce the vertical dimension (denoted by $`z`$) until a planar cavity is obtained . The formation of Fabry-Perot modes in the $`z`$ direction leads to redistribution of spectral weight into peaks of the density of states $`\rho (k)`$ . As a consequence, the spontaneous emission of the microscopic dipoles into the remaining cavity modes is enhanced according to Fermi’s golden rule .

This can lead to directional LED emission  because the resonant cavity modes are spatially anisotropic. In a microcavity laser, such as a VCSEL, the desired directional emission has to be balanced against the additional goal of obtaining laser oscillation at low threshold power. These two requirements are incompatible if we consider the conventional threshold condition that gain and loss must balance out: reducing the linear cavity dimension also reduces the available gain medium, and this eventually leads to an explosive increase in pump threshold if the mirror reflectivities are kept fixed . On the other hand, increased mirror reflectivity means that the available output power is small.

Thus, it is desirable both for LED and laser applications to have as much freedom as possible in the “mirror design” of a microcavity, to create the analogue of the well-known stable, unstable or confocal configurations of macroscopic laser physics . Clearly, vertical layer structures do not lend themselves to such a level of electromagnetic engineering. However, the lateral shape of a semiconductor microcavity can be chosen arbitrarily.

To calculate the cavity mode structure, recall that the laser oscillation condition implies that the wave equation has a solution with outgoing, but no incoming waves in the far field. This radiation boundary condition can be satisfied only at a discrete set of real numbers $`(k,\gamma )`$, where $`k`$ is the wavenumber and $`\gamma `$ the exponential gain constant  of the active medium. This eigenvalue problem can equivalently be stated as a search for complex wavenumbers $`\stackrel{~}{k}=ki\gamma `$ satisfying the above boundary condition with the index $`n`$ of the cavity at its transparency value; the solutions are known in scattering theory as the quasibound states . This is what we mean by the term “modes” of the leaky cavity.

A possible lateral geometry for a VCSEL is that of a circular cylinder, cf. Fig. 1. It allows comparatively simple modeling because each of its layers can be considered as a short section of a cylindrical waveguide for which the modes are known analytically . In particular, there are two types of modes in a cylinder, guided modes whose propagation constant $`\beta _z`$ in the z direction is large, and leaky modes with intensity localized close to the cylinder side wall. When piecing together the fields in the individual layers of the VCSEL stack , it is not surprising that we encouter both type of modes in pixelized VCSELs . However, in a finite-sized, three-dimensionally confined dielectric cavity all modes are leaky, and in fact the large-$`\beta _z`$ modes can turn out to have lower Q than the ones with small $`\beta _z`$, when leakage through the planar mirrors outweighs the lateral radiation losses. This was confirmed in the threshold behavior of a circular AlGaAs multiquantum well (MQW) cavity  of vertical thickness $`L_z\lambda `$ and radius $`R=15\mu `$m.

At large $`\beta _z`$, the sidewalls allow waveguiding  and hence there is predominantly vertical emission, with the planar mirrors controlling the cavity Q factor. On the other hand, small $`\beta _z`$ leads to total internal reflection at the planar mirrors according to Eq. (1), because the angle of incidence satisfies $`\mathrm{tan}\chi =\sqrt{n^2k^2\beta _z^2}/\beta _z`$, where $`k`$ is the free-space wavenumber, and the numerator describes the in-plane component of the internal wave vector. The lateral leakage is also controlled by total internal reflection, provided the light circulates close to the side wall close to grazing incidence; this is just the ray-optics equivalent of the ring-shaped intensity patterns which are commonly known as whispering-gallery (WG) modes, due to an acoustic analogue described first by Lord Rayleigh .

These ring-shaped modes thus avoid refractive escape according to Eq. (1) at all surfaces of the optical cavity. Hereafter, we define $`\chi `$ to be the angle of incidence with respect to the side wall normal to focus on the lateral emission properties. Neglecting absorption in the cavity, the WG mode Q-factor is limited only because internal reflection at the side walls is in fact frustrated by the finite surface curvature. This wavelength-dependent correction to the ray picture can be calculated straightforwardly within the WKB approximation , and leads to infinite Q in the $`\lambda 0`$ limit. The WKB method is here applied to the radial part of the wave equation governing the cavity field (after introducing cylinder coordinates); although it is a short-wavelength approximation, it has been found to provide quantitative results even in microcavities of dimensions comparable to the wavelength .

In contrast to a VCSEL mode, WG modes of a circular cavity emit predominantly sideways, and without any prefered azimuthal (in-plane) direction. For such edge emitting devices, Bragg mirrors on top and bottom of the cavity may be omitted , provided the index-difference between the guiding layer and surrounding cladding is large enough. The strongest index contrast clearly results if we manage to suspend the guiding layer in air. In this way, we get from the circular VCSEL to the microdisk laser : like a thumbtack, the semiconductor disk is supported in the center by a pedestal which is thin enough to have only small or no interaction with the WG modes. Both optically and electrically pumped microdisks have been realized . For better heat-sinking at room-temperature, one can also form the WG disk on pillars of other low-index material, e.g. AlxOy .

Rotational symmetry alone is not sufficient to permit exact analytic solution of the microdisk problem, however. The field of a given cavity mode depends on the cylindrical polar angle $`\varphi `$ only through $`A\mathrm{exp}(im\varphi )+B\mathrm{exp}(im\varphi )`$, where the integer $`m`$ is the analogue of a quantum mechanical angular momentum index. However, the radial ($`r`$) and $`z`$ dependence of the field remain coupled by the combination of boundary conditions at the vertical and lateral interfaces . This makes it necessary either to perform numerical calculations in the $`rz`$ plane , or to neglect this coupling by making a semivectorial approximation  which distinguishes between TE and TM modes. It is in fact confirmed numerically that the electric field of a given mode is either perpendicular to (TM) or in the disk plane (TE), up to small corrections.

If the disk thickness $`L_z`$ is small enough to satisfy $`\sqrt{n^21}kL_z<\pi `$, then only the fundamental slab-waveguide standing waves in the vertical direction are supported for both, TE and TM . The corresponding range of thicknesses for semiconductor disks ($`n>3`$) is therefore approximately $`L_z<\lambda /5`$. The TE modes have higher Q than TM modes because the electric field is better confined in the vertical direction for TE; reported values for TE modes reach as high as $`Q=12000`$ in a $`2\mu `$m diameter disk with InAs self-assembled quantum dots (“boxes”) as emitters . More typical Q factors are in the range from $`Q=100\mathrm{}1000`$; circular disks with such quality have been used to obtain lasing not only from conventional MQW structures  but also from InAs quantum wires in InP  where optical pump powers of $`2`$mW were found at room temperature. In GaN, the introduction of a microdisk cavity led to a pump threshold reduction by an order of magnitude . When only a single quantum dot acts as the emitter in the disk, lasing is not achieved but highly controlled emission of single photons under time-periodiv pumping has been demonstrated; the spontaneous-emission enhancement in the WG microcavity here serves to reduce the timing jitter .

The fact that a microdisk laser emits from its edge could be thought to cause extreme vertical spreading of the output radiation, if we recall that a slab waveguide of thickness $`d`$ has a divergence angle of order $`\pi d/\lambda `$: with disks as thin as $`0.2\mu `$m emitting at $`\lambda =1.55\mu `$m, the spreading should be prohibitive. However, this is not the case . The reason is that the far-field of a microdisk laser cannot be obtained from the near-field on its surface by a simple Fourier transformation, as would be the case for the planar output mirror of a VCSEL . Rotational symmetry implies conservation of the $`\varphi `$-“angular momentum” $`m`$ between the interior and exterior fields, and this impresses a centrifugal potential onto the free-space propagation of the emitted light in the $`rz`$ plane . This in turn favors radial over vertical propagation, reducing the weight of large $`\beta _z`$ in the emission below that expected for a Gaussian beam. A rough estimate yields  a spreading angle (FWHM) $`\mathrm{\Delta }\theta 2/\sqrt{m}`$.

The largest $`m`$ admitted by a microdisk cavity of radius $`R`$ and effective slab index $`n`$ can be estimated from the condition Eq. (1). One can show that semiclassically, a mode with azimuthal index corresponds to rays with angle of incidence satisfying $`m=nkR\mathrm{sin}\chi `$, and since the maximum of $`\chi `$ is $`\pi /2`$, we find that high-Q modes must have $`mnkR`$. As a consequence, $`\mathrm{\Delta }\theta `$ becomes smaller for larger disk radii, provided they lase on WG modes close to grazing incidence at the side walls.

Whereas vertical focusing of the emission around the disk plane is thus surprisingly efficient, there is no easy prescription for inducing focused azimutal emission within this plane. As was recognized early on in the pioneering experiments of Slusher et al., the rotational symmetry of the disk has to be destroyed . In the infinite space of possible shape deformations, it is challenging to find simple design rules that create preferential emission while at the same time preserving the desired high Q-factor. One type of shape perturbation that can never be ruled out is side wall roughness. It turns out, however, that both wavelength and Q factor of WG modes are quite insensitive to small random perturbations of the circular geometry . Additional perturbative effects can be caused by the pedestal .

When the geometry of the disk is distorted to such an extent that perturbative treatments  break down, numerical modeling becomes much more complex because $`m`$ is no longer a “good quantum number” by which the modes can be labeled. Motivated by the success of the WKB method mentioned above, it was suggested early on that short-wavelength approximations can provide valuable insights into the mode structure and emission patterns of microlasers . The main observation is that the internal ray dynamics of a non-circular disk is in general chaotic, so that methods from the fields of classical and quantum chaos theory become applicable. Three chaotic WG modes are shown in Fig. 2, for a wide range of different wavelengths. At the deformation considered here, chaotic ray motion is the dominant mechanism for escape from the cavity, as opposed to diffractive ($`\lambda `$-dependent) effects. This is evidenced by the fact that $`\gamma `$ in Fig. 2 depends only weakly on wavenumber $`k`$. Likewise, one can identify wavelength-independent features in the internal and external intensity profiles; the emission directionality at large deformations is one of these characteristics, which to lowest order can be explained without recourse to the wave nature of light.

The first experimental tests of these predictions were performed on lasing dye microdroplets whose shape is generically oval . When rotational symmetry is absent, $`\chi `$ in Eq. (1) is not conserved between reflections, so that diffusion of photons away from the WG region tends to occur; the modal decay rate then generally increases with deformation because even WG rays can reach Eq. (1) after sufficiently many loops in the cavity.

However, the Q-spoiling due to large deformation is less severe than expected from ray considerations alone, as was pointed out in Ref. . An interference effect called “dynamical localization” actually suppresses decay rates in a manner analogous to Anderson localization in disordered solids , thus making high-Q asymmetric cavities more feasible: the advantages of asymmetric shapes for directional emission can then outweigh the loss in cavity Q.

Experimental realizations of non-circular microdisks  confirm another important prediction of the chaos analysis : oval shapes favor WG lasing only as long as the boundary remains everywhere convex, i.e., has no points of vanishing curvature. The lateral emission from oval disks is found to emanate from the regions of highest curvature in semiconductor disks. The light furthermore escapes tangentially to the surface, as expected for rays at the critical angle of incidence $`\chi _c`$ just satisfying Eq. (1).

A striking exception to this rule is expected in similar structures made of polymers  or other low-index materials: the chaotic ray dynamics in such oval cavities  can counter-intuitively prevent WG rays from reaching $`\chi _c`$ at the high-curvature points. This “dynamical eclipsing” occurs because ray chaos is not random, but deterministic; its underlying structure can in particular rule out certain $`\chi `$ for WG trajectories approaching the high-curvature points. By changing the material (refractive index) of the cavity to bring $`\chi _c`$ into this “forbidden” range, the emission pattern for the same shape as in Fig. 2 then changes significantly.

In contrast to smooth oval deformations, one can also introduce notches or projections in the side walls, but this tends to increase vertical spreading , thus reducing the expected focusing of the laser. Recalling the dependence of vertical spreading $`\mathrm{\Delta }\theta `$ on angular momentum, we should aim for smoothly deformed microdisks emitting near the critical angles of incidence, yet as anisotropic as possible. At large deformations, cavity modes with this property can have internal intensity distributions that are qualitatively different from WG modes.

This can be illustrated by a series of experiments on oval lasers, incrementally ramping up the deviation from circular shape. Such measurements have been performed in quantum-cascade microcylinders in the shape of a flattened quadrupole . The devices were intended for high power and hence had larger area than the examples above: keeping the minor axis at $`b=50\mu `$m diameter and increasing the major axis to $`a=80\mu `$m, it was found that beyond a threshold deformation $`a70\mu `$m, an exponential increase in the maximum achievable output power set in. The lasing mode was identified as belonging to a closed “bowtie” ray pattern which does not exist below the threshold deformation.

The effect of chaos in this experiment is to spoil all cavity modes overlapping the gain region, except for the bowtie modes. This reduces the density $`\rho `$ of available modes, while at the same time the internal intensity profile of the remaining modes becomes more focused in the electrically pumped center of the structure. The four vertices of the bowtie pattern meet the side walls near critical incidence, resulting in directional emission. The combined power of the four emission lobes of up to 10 mW was $`1000`$ times larger than that of the circular device in the test series.

A special feature of the quantum-cascade system that makes it amenable to this near-critical mode is that (unlike the previous examples) it strongly favors emission in TM polarization. This is important because in a high-index cavity the critical angle for total internal reflection is close to the Brewster angle; therefore, near-critical TE modes suffer from much higher lateral emission than their TM counterparts (the reverse holds for the vertical losses).

A reduced $`\rho `$ permits reduced lasing thresholds, and the lesson learned from the bowtie laser is that chaos can reduce the mode density even in relatively large-area cavities. This is because the average of $`\rho `$ is not determined by the cavity volume, but by the available phase space volume for long-lived ray trajectories. The phase space occupied by a ray trajectory is determined by the real-space area it explores, combined with the range of $`\chi `$ it spans. We have recently shown  that in a WG-cavity of mean radius $`R`$ (cf. Fig. 2) and index $`n`$, the average number of long-lived (non-overlapping) modes per unit wavenumber interval is

$$\rho (k)dk=(1/4)n^2kR^2\left[1(2/\pi )\left(\mathrm{arcsin}(1/n)+(1/n)\sqrt{1(1/n)^2}\right)\right]dk.$$
(2)
The first term is determined by the cavity area, but the corrections for small $`n`$ show how the openness of the cavity reduces the density of modes by expanding the classical escape window Eq. (1) which depends on $`\chi `$.

In conclusion, phase space reduction and hence lower threshold can be achieved by shrinking the real area of the cavity, e.g. with holes in the disk  or ring structures , or alternatively by engineering the side wall deformation to obtain stable modes such as the bowtie while simultaneously destablilizing most others. Concepts from nonlinear dynamics provide design guidelines for the latter approach. The rich structure provided by the ray dynamics in chaotic cavities makes it possible to design shapes for laterally directional emission, which makes it possible to collect the emitted light more efficiently than in rotationally symmetric structures.


text:
# Probing hadronic formation times with antiprotons in 𝑝+𝐴 reactions at AGS energies Supported by the European Graduate School ’Complex Systems of Hadrons and Nuclei’, Copenhagen-Giessen

## 1 Introduction

Since the first observation of antiproton production in proton-nucleus  and nucleus-nucleus collisions  the production mechanism has been quite a matter of debate. Especially in nucleus-nucleus collisions at subthreshold energies traditional cascade calculations, that employ free $`NN`$ production and $`\overline{p}N`$ annihilation cross sections, essentially fail in describing the high cross sections seen from 1.5 - 2.1 A$``$GeV . Thus multiparticle nucleon interactions  have been suggested as a possible solution to this problem. On the other hand it has been pointed out that the quasi-particle properties of the nucleons and antinucleons might be important for the $`\overline{p}`$ production process which become more significant with increasing nuclear density . Relativistic transport calculations have indicated that all the low energy data from proton-nucleus and nucleus-nucleus collisions are compatible with attractive $`\overline{p}`$ self energies (at normal nuclear matter density $`\rho _0`$) in the order of -100 to -150 MeV . However, it had been stressed at that time that the high antiproton yield might also be attributed to mesonic production channels  since $`p\overline{p}`$ annihilation leads to multi-pion final states with an average abundancy of 5 pions .

With new data coming up on antibaryon production from proton-nucleus and nucleus-nucleus collisions at the AGS  and SPS  the $`\overline{p},\overline{\mathrm{\Lambda }}`$ enhancement factors seen experimentally were no longer that dramatic as at SIS energies, however, traditional cascade calculations employing free production and annihilation cross sections again were not able to reproduce the measured abundancies and spectra  especially for $`\mathrm{\Xi },\overline{\mathrm{\Xi }}`$ and $`\mathrm{\Omega },\overline{\mathrm{\Omega }}`$. Here additional collective mechanisms in the entrance channel have been suggested such as color rope formation  or hot plasma droplet formation . In another language this has been addressed also as string fusion , a precursor phenomenon for the formation of a quark-gluon plasma (QGP).

The intimate connection of antibaryon abundancies with the possible observation of a new state of the strongly interacting hadronic matter, i.e. the quark-gluon plasma, has been often discussed since the early suggestion in Ref.  that especially the enhanced yield of strange antibaryons – approximately in chemical equilibrium with the other hadronic states – should be a reliable indicator for a new state of matter. In fact, the data on strange baryon and antibaryon production from the NA49 and WA97 Collaborations show an approximate chemical equilibrium  with an enhancement of the $`\mathrm{\Omega }^{},\mathrm{\Omega }^+`$ yield in central $`Pb+Pb`$ collisions (per participant) relative to $`pBe`$ collisions at the same invariant energy per nucleon by a factor $``$ 15. As pointed out in Ref.  the data on multi-strange antibaryons at the SPS seem compatible with a canonical ensemble in chemical equilibrium.

In Ref.  Rapp and Shuryak have proposed multi-meson production channels for $`p\overline{p}`$ pairs to describe the antiproton abundancies in central $`Pb+Pb`$ collisions at the SPS. Later on, Greiner and Leupold  have applied the same concept for the $`\overline{\mathrm{\Lambda }}`$ production by a couple of mesons including a $`K^+`$ or $`K^0`$ (for the $`\overline{s}`$ quark). In fact, their estimates have been supported by the microscopic multi-particle calculations in Ref.  employing detailed balance relations for 2 hadron $``$ n hadron transitions and vice versa. Thus, in case of relativistic nucleus-nucleus collisions at AGS and SPS energies the low effective antibaryon absorption seen experimentally can naturally be explained by the backward reaction channels, i.e. multi-meson fusion channels.

The latter mechanism, however, does not apply to proton-nucleus reactions since there is no longer an approximately thermal ’meson bath’ feeding the baryon-antibaryon production channel as pointed out in Ref. . Nevertheless, the $`\overline{p}`$ yield measured by the E910 Collaboration in $`p+A`$ reactions at 12.3 GeV/c and 17.5 GeV/c was found to be surprisingly high . An analysis within the Glauber model showed that either the annihilation cross section $`\sigma _{ann}`$ should be small compared to the cross section in free space or the $`\overline{p}`$ formation time $`\tau _F`$ should be a couple of fm/c when employing the free annihilation cross section . The authors of the latter work conclude that the ’naive’ first-collision model with subsequent annihilation of antiprotons is incompatible with their data supporting the earlier finding from Ref. , where previous data from $`p+A`$ and $`A+A`$ reactions had been analysed within a similar model.

In this work we will address $`\overline{p}`$ production in proton-nucleus collisions at 12.3 GeV/c and 17.5 GeV/c within the covariant HSD111Hadron-String-Dynamics transport approach  which is based on the concept of string formation and decay at invariant collision energies $`\sqrt{s}`$ 2.6 GeV for baryon-baryon collisions and $`\sqrt{s}`$ 2.3 GeV for meson-baryon collisions, while the low-energy reactions are described by known vacuum cross sections or resonance formation and decay cross sections. A single sensitive parameter in this approach is the string duration or hadron formation time $`\tau _F`$ 0.7-0.8 fm/c which has been fitted to the proton rapidity distribution in nucleus-nucleus collisions at SPS energies  and kept fixed furtheron. For more details we refer the reader to the review .

## 2 Qualitative considerations

The physical picture within this transport approach for $`\overline{p}`$ production in $`p+A`$ reactions is as follows: The impinging proton might create a $`N\overline{N}`$ pair in a first collision at position $`x_1`$ (cf. Fig. 1) and the prehadronic quark-antiquark configuration travels essentially in forward direction by a time $`\gamma \tau _F`$ – where $`\gamma `$ denotes the Lorentz dilatation factor – until the hadronic $`\overline{p}`$ emerges at position $`x_2`$. The hadronic state then propagates further in the same direction and reacts with the remaining nucleons elastically or inelastically. In the latter phase the dominant channel will be annihilation with a cross section

$$\sigma _{ann}(\sqrt{s})=\frac{50[mb]}{v_{rel}},$$
(1)
where the relativistic relative velocity is given by

$$v_{rel}=\frac{\sqrt{(sM_1^2M_2^2)^24M_1^2M_2^2}}{2E_1E_2},$$
(2)
which is evaluated in the lab. frame. As shown in Fig. 3 of Ref.  the expression (1) well describes the experimental data  on $`\overline{p}`$ annihilation for relative momenta up to 2–3 GeV/c.

In rapidity space the produced antibaryons will be localized around midrapidity $`y_m`$; fast antiprotons hadronize later than slow $`\overline{p}`$’s due to the larger $`\gamma `$ factor and in view of (1) also experience a smaller annihilation cross section. Thus antiprotons from first-chance collisions in heavy nuclei should result in rapidity distributions that are highly asymmetric with respect to midrapidity and a larger suppression should be observed for $`yy_m`$ than for $`yy_m`$, if the target is located at $`y=`$ 0. In fact, a simple look at the experimental $`\overline{p}`$ rapidity distribution from Fig. 2 of Ref.  does not support the first-collision picture.

On the other hand, in Ref.  it was found that a dominant fraction of antiprotons emerge from secondary pion-nucleon reactions in $`p+A`$ reactions at lower energies. Such secondary reaction channels might also be a reason for the observations of Ref.  since the $`\overline{p}`$ production cross section in $`\pi N`$ reactions close to threshold is large compared to the production cross section in $`pp`$ collisions (cf. Fig. 6 of Ref. ). The dymical effects of such secondary production channels are easy to map out within the string/hadron dynamical picture: In a first $`pN`$ collision a ’meson’ ($`q\overline{q}`$ pair) $`M`$ is produced at position $`x_1`$ which hadronizes after time $`\gamma _M\tau _F`$ at position $`x_2`$. Then it propagates on average for a distance $`\lambda _M(\sigma _{tot}\rho _0)^1`$ 2 fm (in the lab. frame of reference) before reacting with another nucleon. In this reaction (at position $`x_3`$) it may also produce a pre-$`N\overline{N}`$ pair which hadronizes after time $`\gamma \tau _F`$ in position $`x_4`$. Then the antinucleon propagates in the same direction and may be rescattered, annihilated or emerge to free space. The antibaryons from secondary – or even ternary – processes have a much shorter way to escape the nucleus without being annihilated than those from first-chance collisions. Furthermore, their initial rapidity distribution will be shifted closer to target rapidity because the primary mesons are produced symmetrically around $`y_m`$ and only the energetic mesons will give a sizeable contribution to the $`\overline{p}`$ production provided that they hadronize in the medium.

As mentioned above we here adopt the HSD transport approach that has e.g. been tested in Ref.  for pion production and proton stopping in $`p+A`$ reactions at AGS energies using $`\tau _F`$ = 0.8 fm/c. Since the description of the related data is very satisfactory (cf. Fig. 8 in Ref. ) we continue with aspects of meson formation and secondary interactions, that are relevant for antiproton production and propagation in $`p+A`$ reactions.

In order to achieve a quantitative idea about the time scales for meson production and propagation we show in Fig. 2 for $`p+Au`$ at 12.3 GeV/c the time dependence of the ’formed’ meson number (solid line) in comparison to the number of string-like $`q\overline{q}`$ excitations (dashed line) that convert to mesons on a time scale of $`\tau _F`$ in their individual rest frame. The reference frame for the transport calculation is the nucleon-nucleon cms. (at $`y_{lab}1.63`$) with a Lorentz $`\gamma `$-factor of $`2.67`$ for the projectile proton and target nucleus, respectively, such that the dimensions of the target in space - as seen by the hadrons - are squeezed by the factor $`\gamma `$. The sharp increase of the dashed line at t= 0 denotes the beginning of the $`p+Au`$ reaction with the excitation of strings that may lead to the formation of $`N\overline{N}`$ pairs, too. The latter also appear as hadronic states on a similar time scale as the mesons – approximately 2 fm/c later – and may be annihilated or rescattered on the residual target nucleons. We recall that only the formed hadronic states are allowed to interact with the nucleons of the target again and to produce secondary mesons or baryon-antibaryon pairs. Thus antiprotons or mesons, that move fast with respect to the target at rest, may only form in the vacuum and not interact in the nucleus again. In fact, this is approximately the case for a $`Be`$-target at this energy, while the dominant fraction of hadrons are formed inside the $`Au`$-nucleus and reinteract again (see below). For completeness we note, that the increase of the meson number for $`t7`$ fm/c is essentially due to the decays of mesonic resonances ($`\rho ,\omega ,\varphi ,K^{}`$) that increase the net meson number.

The distribution in the invariant collision energy $`\sqrt{s}`$ is shown in Fig. 3 for $`p+Au`$ at 12.3 GeV/c for baryon-baryon ($`BB`$) collisions (solid line) in comparison to the reactions of formed mesons with nucleons ($`mB`$, dashed line). The high energy peak of the $`BB`$-distribution (around $`\sqrt{s}`$ 5 GeV) is due to primary $`pN`$ collisions; it is smeared out substantially due to Fermi motion. The low energy peak of the $`BB`$-distribution reflects multiple scattering of baryons in the target that, however, do not contribute to $`p\overline{p}`$ production due to the threshold of $``$ 3.75 GeV (solid arrow). The $`\sqrt{s}`$ distribution for interactions of formed mesons with baryons (dashed line) extends to about 4.5 GeV, which is well above the $`\overline{p}p`$ production threshold of $``$ 2.8 GeV (dashed arrow). The distribution in the collision number above threshold is substantially larger for $`BB`$-reactions than for meson-nucleon collisions, however, the distributions have to be ’folded’ with the respective $`p\overline{p}`$ production cross sections. We recall that in $`\pi N`$ reactions close to threshold this cross section is large compared to the production cross section in $`pp`$ collisions (cf. Fig. 6 of Ref. ) such that the $`mB`$ channels might well compete with the $`BB`$ channels.

## 3 Comparison to data

It should be noted that for the present analysis the perturbative production scheme of Ref.  has been employed due to the low multiplicities for $`\overline{p}`$ production ($`2\times 10^4`$ for $`p+Au`$ at 12.3 GeV/c). All cross sections are the same as in the latter work; the only modification introduced is a new parametrisation of the individual rapidity $`(y)`$ and transverse momentum $`(p_T)`$ dependence for antiprotons in the elementary $`pN`$ and $`\pi N`$ production channels, since the kinematical regime is substantially different from the latter study and the final states can no longer be described accurately by 4- or 3-body phase-space, respectively. To this aim high statistics event generation has been performed within the LUND string model  – as implemented in the HSD approach – for $`\overline{p}`$ production and the resulting two-dimensional spectra have been fitted by a Gaussian in rapidity and exponential in $`p_T^2(y)`$. As demonstrated in Refs.  the LUND model works rather well for this purpose in comparison to the available data , however, uncertainties in the order of 50% for the elementary cross sections cannot be ruled so far since there are no data available for invariant energies close to threshold.

With the input cross sections defined we can step on with a quantification of the physics addressed above. In Fig. 4 we show the $`\overline{p}`$ rapidity spectra for $`p+Be`$ at 12.3 GeV/c in comparison to the experimental data from . The dashed line ($`BBprim.`$, upper part) gives the $`\overline{p}`$ spectrum from $`BB`$ production channels prior to final state interactions (FSI). The contribution from secondary $`mB`$ collisions is not visible on a linear scale for the $`Be`$ target. When including $`\overline{p}`$ rescattering as well as annihilation the dash-dotted line ($`BBfinal`$) is obtained for $`BB`$ production channels. It is clearly seen that $`\overline{p}`$ annihilation is more pronounced for low rapidities – with the target at $`y=0`$ – than for $`yy_m`$ 1.63. Thus for $`p+Be`$ the $`pN`$ production channel clearly dominates and the comparison with the experimental rapidity distribution (lower part of Fig. 4) essentially tests the adequacy of the elementary differential $`\overline{p}`$ production cross section in $`pN`$ collisions which appears to be underestimated slightly (by about 30%) in the present approach using the input from Ref.  and assuming that the data point for $`y`$ = 1.1 is solid. One might rescale the elementary cross section accordingly, however, the focus of our present work is not on fitting data. On the other hand, antiprotons from a feeddown of $`\overline{\mathrm{\Lambda }}`$ and $`\overline{\mathrm{\Sigma }}^0`$ might also be contained in the data while they are excluded in the calculation222In Ref.  such a feeddown contribution was estimated to be less than 5%.. Thus, provided that the reactions are described properly in the HSD transport approach, the cross sections for all targets should be low by about 30%, too, in case the $`\overline{p}`$ yield is dominated by the $`BB`$ channels (see below).

The numerical results for $`p+Cu`$ and $`p+Au`$ at 12.3 GeV/c are shown in Figs. 5 and 6 adopting the same assignment of the lines as in Fig. 4. It is clearly seen that for the $`Cu`$-target the contribution from $`mB`$ collisions – as displayed by the short dashed line ($`mBprim.`$) – becomes visible and that the annihilation of $`\overline{p}`$’s from the secondary interactions is rather small as expected from the schematic Fig. 1. On the other hand, in case of the $`Au`$-target the antiprotons from $`NN`$ channels are annihilated to a large extent especially for $`yy_m`$ and the surviving $`\overline{p}`$’s in the experimental rapidity window $`1.0y2.0`$  practically stem from the secondary production channels. Nevertheless, the experimental data are underestimated on average again by about 30% as in case of the $`Be`$-target which might be due to improper elementary cross sections and/or still missing production channels in the transport approach.

The relative antiproton absorption according to Figs. 4–6 strongly depends on rapidity due to a larger formation time $`\tau _F`$ and lower annihilation cross section with increasing rapidity $`y`$. In case of $`p+Be`$ there is practically no $`\overline{p}`$ annihilation for $`y`$ 2 since the fast antinucleons only form in the vacuum. For the $`p+Cu`$ system this threshold is shifted to $`y`$ 2.6 while for the $`p+Au`$ system the annihilation becomes negligible only for $`y`$ 3. The sensitivity to the formation time $`\tau _F`$ is demonstrated in the lower part of the Figs. 5 and 6, where also the final antiproton rapidity spectra are displayed for $`\tau _F`$ = 0.4 fm/c (dashed line) and 1.6 fm/c (dotted line), respectively. In case of the light $`Be`$-target no sensitivity to $`\tau _F`$ was found within the numerical statistics of the transport calculations since there is practically no contribution from the $`mB`$ channels. As can be extracted from Figs. 5 and 6 a shorter formation time $`\tau _F`$ leads to an enhancement of the secondary production channel by $`mB`$ collisions, whereas a larger formation time ($``$ 1.6 fm/c) essentially suppresses the contribution from secondary meson induced reactions in line with Fig. 1. In fact, the experimental spectra for $`Cu`$ and $`Au`$ appear to be better described by $`\tau _F`$ = 0.4 fm/c for the elementary cross sections employed.

The average $`\overline{p}`$ multiplicity in the rapidity interval $`1y2`$ at 12.3 GeV/c is compared in Fig. 7 as a function of target mass $`A`$ with the data from  for different $`\tau _F`$. As mentioned before, for the $`Be`$-target the calculations for $`\tau _F`$ = 0.4, 0.8 and 1.6 fm/c give approximately the same result, which in comparison to the experimental data is down by about 30% just scratching the lower level of the experimental error bars. For the $`Cu`$-target both calculations for $`\tau _F`$ = 0.8 fm/c and $`\tau _F`$ = 1.6 fm/c fall down too low. In case of the $`Au`$-target only the calculation for $`\tau _F`$ = 1.6 fm/c is out of range in comparison to the data and also incompatible with the shape of the $`dN/dy`$ spectra such that a formation time in the order of 1.6 fm/c should be excluded.

We, furthermore, compare our calculations to the system $`p+Au`$ at 17.5 GeV/c where the E910 Collaboration has taken data, too . In this case (Fig. 8) the most important production channel is due to primary $`pN`$ inelastic collisions, however, the surviving antiprotons dominantly appear above midrapidity $`y_m`$ 1.81 due to the strong absorption. Within the rapidity window $`1y2`$ of the data the contribution from secondary $`mB`$ channels is almost compatible (thin solid line, $`mBfinal`$, compared to the dot-dashed line, $`BBfinal`$). The sum of the contributions is reasonably well in line with the experimental data from  (lower part) for $`\tau _F`$ =0.8 fm/c (solid line). A slightly better description is achieved with $`\tau _F`$ = 0.4 fm/c for $`y`$ 1-1.5 in line with the data at 12.3 GeV/c, however, now the $`\overline{p}`$ rapidity distribution is slightly overestimated for $`y`$ 2. In view of the fact, that the elementary antiproton production cross section is underestimated by about 30% in our calculations, this provides a lower limit for the formation time $`\tau _F`$.

## 4 Summary

In this work we have analyzed antiproton production in $`p+A`$ reactions at AGS energies within the HSD transport approach to examine the apparent low $`\overline{p}`$ annihilation seen experimentally . Whereas the rather high antibaryon yield in nucleus-nucleus collisions at AGS and SPS energies can be attributed to multi-meson fusion channels employing detailed balance  the experimental data in proton-nucleus collisions at AGS energies appear to be compatible with primary proton-nucleon and secondary meson-baryon production channels as also found in Ref.  for lower energies.

Detailed calculations employing different hadronic formation times rule out long formation times of the order of 1.6 fm/c (in the hadron rest frame) since such times strongly suppress the meson reinteraction rate. The comparison with the data from the E910 Collaboration instead favors formation times of 0.4-0.8 fm/c when adopting the free antibaryon annihilation cross section. The latter cross section might be somewhat reduced in the nuclear medium, but it is not likely to be down by factors $``$ 5 as quoted in the analysis of Ref. . We have pointed out, furthermore, that especially the high momentum/rapidity part of the $`\overline{p}`$ spectrum, which dominantly arises from primary $`pN`$ collisions, shows a sensitivity to the $`\overline{p}`$ formation time since ’fast’ antiprotons may only form in free space and no longer get annihilated in the medium. Unfortunately, this rapidity regime is out of the E910 acceptance such that further experimental studies will be necessary to provide a final answer to the issue of antiproton formation times.

The authors like to acknowledge valuable discussions with C. Greiner.


text:
# Numerical simulations of the cosmic star formation history

## 1 Introduction

A detailed understanding of galaxy formation remains one of the primary goals of modern cosmology. Whereas gravitation is basically responsible for the large scale structure of the universe, many physical phenomena, covering a broad dynamical range, are involved on galactic scales. A reliable treatment of baryonic physics is essential in order to obtain a consistent picture of galaxy formation and evolution, where star formation and feedback to the interstellar medium (ISM) play a key role.

The star formation history of the universe, i.e. the global star formation rate (SFR) density as a function of redshift, is a crucial test for galaxy formation scenarios. Many measurements of this quantity have been undertaken during the past few years, and the number and accuracy of available observations are still rapidly increasing.

Photometric redshift estimates based on Lyman-limit systems (Steidel et al. 1996a ; 1996b ) stimulated the study of the rest frame UV continuum flux around 2000 Å as a preferred wavelength range for objects in the far universe. This emission traces the presence of massive (and therefore young) stars and can be related directly to the actual star formation rate. It is relatively easy to detect at high $`z`$ as it is redshifted into the optical band. However, older stellar populations and AGN can make a significant contribution as well, leading to an overestimation of the total SFR. On the other hand, dust enshrouding star formation regions absorb very efficiently the UV light, re-emitting it in the far IR. Unfortunately, there is still some uncertainty to account for dust extinction and its evolution with redshift (Nagamine et al. (2001)).

Recombination lines from H II regions can be in principle a more reliable estimator of the instantaneous SFR, since the ionizing radiation ($`\lambda 912`$Å) comes from more massive (younger) stars than the softer UV continuum. These optical lines are less dramatically affected by dust extinction and can be mapped with higher resolution in the local universe, but they require infrared spectroscopy at moderate redshifts, which has not been possible until very recently due to the faintness of the sources involved (Pettini et al. (1998)). For Hα, this happens at $`z0.5`$, and for forbidden lines, such as \[O II\]3727Å, at $`z1.5`$. There is also a considerable degree of uncertainty in the relation between recombination lines and star formation activity.

Finally, the far infrared (FIR) spectrum of a galaxy can be roughly separated into two components, namely the thermal emission of dust (heated by star formation bursts) around $`\lambda 60\mu `$m and an infrared cirrus powered by the global radiation field, which dominates for $`\lambda 100\mu `$m. Neglecting cirrus and AGN contamination, the FIR luminosity would also be an excellent tracer, but instrumentation is not as developed in the IR and radio bands as it is in the optical, and the conversion factor between FIR luminosity and SFR is for the moment rather model dependent.

Although there is a large scatter among different indicators at any given redshift, basically all studies find a significant increase (by about one decade) in the cosmic SFR density from the present day to $`z=1`$. However, it is still a controversial issue (Hopkins et al. (2001)) whether it reaches a broad maximum there (Madau et al. (1996); Gispert et al. (2000)) or it declines gradually towards higher redshifts (Sawicki et al. (1997); Pascarelle et al. (1998); Steidel et al. (1999)).

From a theoretical point of view, star formation is a highly non-linear process, which precludes any simple analytical treatment. Most efforts towards modeling the star formation history of the universe resort to semi-analytical methods or numerical simulations to tackle the feedback mechanisms that self-regulate the SFR within the hierarchical clustering scenario of structure formation.

Semi-analytical models (Kauffmann et al. (1993); Cole et al. (1994)) construct a large sample of Monte Carlo realizations of halo merging histories, using the Press-Schechter formalism. Hydrodynamics, star formation and feedback are implemented through simple recipes, whose parameters are fixed in order to match the observed properties of real galaxies, specially the luminosity function and the Tully-Fisher relation (Somerville & Primack (1999)). Their high computational efficiency allows a fast exploration of the parameter space in search of a viable model, but the number of free parameters asks for supplementary modeling based on more fundamental physics.

Numerical simulations are intended to solve directly the set of differential equations of gravitation and gas hydrodynamics, and hence require much fewer model assumptions than the semi-analytic approach. However, despite the staggering advances in computer technology and numerical algorithms, resolution in large-scale hydrodynamical simulations is severely restricted by computational and data storage requirements. Although gas can be treated self-consistently, some phenomenological recipes are still required in order to model star formation and feedback.

The cosmic SFR density is a direct outcome of semi-analytical methods (Somerville et al. (2001)) and hydrodynamical simulations (Nagamine et al. (2001)). Both of them predict a significant star formation activity at high $`z`$, favoring the hypothesis of a gradual evolution of the comoving SFR as a function of time.

The present analysis is based on a series of numerical simulations, within different cosmological models. In previous papers, we have studied basic observable properties of simulated galaxy samples, such as the luminosity function and color dependence on environment (Elizondo et al. 1999a ), and the Tully-Fisher relation (Elizondo et al. 1999b ). In this paper, we will focus on the evolution of the SFR on a cosmological scale, comparing our results with recent observations.

The following section is devoted to the observational estimation of the cosmic SFR. The standard procedure is briefly summarized and a compilation of recent results is included with an analytical prescription for the transformation between cosmological models. Sect. 3 describes the code and the different numerical experiments on which the present paper is based. Our results are presented in Sect. 4, where we compare with observations and evaluate the separate contribution of cosmology, feedback effects and environment. Our main conclusions are summarized in Sect. 5.

## 2 Observations

The first step to obtain an estimate of the SFR density from observations is the selection of a complete sample of galaxies at a given redshift and wavelength. Some authors (Steidel et al. (1999)) claim that the HDF covers a very small area ($`5`$ arcmin2), and hence it is not statistically representative of the whole universe. This cosmic variance is difficult to quantify, but the strong clustering of high redshift galaxies (and the fact that most estimates beyond $`z=2`$ are based on the HDF) seems to indicate that it could be an important effect.

Another crucial issue is the correction for incompleteness in the construction of the luminosity function, which can be done following the $`V_{\mathrm{max}}`$ formalism (Schmidt (1968)) or resorting to more elaborate Monte Carlo algorithms (e.g. Steidel et al. (1999)). The comoving luminosity density is easily obtained by integration over all magnitudes (extrapolating the faint end as a Schechter function), and then it can be converted to a SFR density with the help of a population synthesis model.

A compilation of recent results is given in Table 1. The first column indicates the bibliographic reference, as well as the survey from which the galaxy sample was extracted. Next columns correspond to the chosen wavelength and redshifts of each analysis, and finally the comoving luminosity and SFR densities.

Conversion factors between the last two quantities differ significantly from author to author, and in some cases only the luminosity density was provided. If there is a SFR estimate in the original paper, then it is quoted in Table 1 as corrected by dust extinction; else, we follow Kennicutt’s (1998) prescription for an exponential burst and a Salpeter IMF between 0.1 and 100 M:

$$\dot{\rho }_{}[\mathrm{M}_{}\mathrm{yr}^1]=\{\begin{array}{cc}1.4\times 10^{28}\mathrm{L}_{\mathrm{UV}}[\mathrm{erg}\mathrm{s}^1\mathrm{Hz}^1]\hfill & \\ 1.4\times 10^{41}\mathrm{L}_{3727}[\mathrm{erg}\mathrm{s}^1\mathrm{Hz}^1]\hfill & \end{array}$$
(1)
For dust extinction, we have applied a correction of A(1500-2000Å)=1.2 mag and A(2880Å,O II)=0.625 mag, which correspond to factors of 3.02 and 1.78 respectively. All modifications to the original data are shown as a footnote on the bibliographic reference in the Table 1.

Most observational values given in Table 1 were calculated assuming a standard cold dark matter (SCDM) cosmology with dimensionless Hubble parameter $`h=0.5`$ and mean matter density $`\mathrm{\Omega }_\mathrm{m}=1`$. Since the conversion factor from luminosity to SFR and the correction for dust extinction are completely independent of the cosmological scenario, the computation of the comoving SFR density in a $`\mathrm{\Lambda }`$CDM ($`h=0.7`$) model only involves the transformation of the luminosity densities.

In a flat universe, the volume enclosed by a solid angle $`\mathrm{\Delta }\omega `$ between $`z\mathrm{\Delta }z`$ and $`z+\mathrm{\Delta }z`$ is

$$V(z,\mathrm{\Delta }z)=\frac{4\pi }{3}\mathrm{\Delta }\omega \left(d_\mathrm{m}^3(z+\mathrm{\Delta }z)d_\mathrm{m}^3(z\mathrm{\Delta }z)\right),$$
(2)
where $`d_\mathrm{m}(z)`$ is the comoving distance in each cosmological model, while luminosity scales with luminosity distance as $`d_\mathrm{L}(z)=(1+z)d_\mathrm{m}(z)`$. Then, assuming that all galaxies are located at the centers of the interval in $`z`$, the luminosity density (and hence the total SFR density) will be proportional to

$$\dot{\rho }_{}(z)\frac{L(z)}{V(z,\mathrm{\Delta }z)}\frac{d_\mathrm{m}^2(z)}{d_\mathrm{m}^3(z+\mathrm{\Delta }z)d_\mathrm{m}^3(z\mathrm{\Delta }z)}$$
(3)
To first order in $`\mathrm{\Delta }z`$, the conversion between SFR densities in $`\mathrm{\Lambda }`$CDM and variants of the SCDM ($`h=0.5`$) cosmologies can be approximated as

$$\frac{\dot{\rho }_{}^{\mathrm{\Lambda }\mathrm{CDM}}}{\dot{\rho }_{}^{\mathrm{SCDM}}}\frac{h^{\mathrm{\Lambda }\mathrm{CDM}}}{h^{\mathrm{SCDM}}}\sqrt{\mathrm{\Omega }_\mathrm{m}+\mathrm{\Omega }_\mathrm{\Lambda }(1+z)^3}$$
(4)
Fig. 1 shows this factor, as well as the exact expression for different values of $`\mathrm{\Delta }z`$, which requires a numerical evaluation of the luminosity distance in the $`\mathrm{\Lambda }`$CDM model. Although the results are very similar, we have chosen to use the exact expression (3) to compute the SFRs in Table 1 in order to minimize errors.

At low redshifts, the data are consistent for all the independent tracers of the star formation activity once the effect of dust obscuration is taken into account. The steep increase in the cosmic SFR density from the present day to $`z1`$ is firmly established. However, there is not yet a full agreement about whether there is a characteristic epoch of star formation in the universe around $`z1.5`$ or on the contrary that the SFR declines very slowly towards higher $`z`$, as more recent observations seem to point out. The results of ongoing deep surveys will be extremely helpful in order to reach a definitive conclusion from the observational point of view.

## 3 Simulations

### 3.1 Cosmological Models

Up to now, the $`\mathrm{\Lambda }`$CDM model has proven to be very successful in describing most of the observational data at both low and high redshifts. Moreover, recent observations indicate the existence of a positive cosmological constant (Perlmutter et al. (1999); Melchiorri et al. (2000)). For our study we have chosen a spatially flat universe with a cosmological constant $`\mathrm{\Lambda }/3H_0^2\mathrm{\Omega }_\mathrm{\Lambda }=0.65`$ and a present-day Hubble constant of $`H_0=100h=70`$ km/s/Mpc. In order to study the influence of cosmology on the global SFR, we have investigated three other scenarios: $`\mathrm{\Lambda }`$CDM with lower $`\mathrm{\Omega }_\mathrm{\Lambda }=0.6`$, the standard CDM model (dominated by dark matter density), and the BSI model with a Broken Scale Invariant perturbation spectrum (as predicted by double inflation, Gottlöber et al. (1991)). The latter model is also dark matter-dominated, but it describes a primordial power spectrum in which fluctuations on small scales are significantly suppressed. It is almost the same spectrum as in the $`\tau `$CDM model, in which the primordial perturbation spectrum has been changed due to the hypothetical decay of massive $`\tau `$ neutrinos (Efstathiou et al. (1992)). The model parameters are summarized in Table 2. All spectra are COBE-normalized. The normalization is given by the linear dark matter mass variance $`\sigma _8`$ in a $`8h^1`$ Mpc sphere extrapolated to present time.

### 3.2 Numerical code

The code used to performed the numerical experiments consists of a combination of a Particle-Mesh N-body gravity solver and an Eulerian Piecewise Parabolic Method to account for gasdynamics. A detailed modeling of the most relevant non-adiabatic baryonic processes (i.e. multiphase ISM, cooling, photoionization, star formation and supernovae feedbacks) acting at a sub-grid level is also incorporated. The code is described in detail elsewhere (Yepes et al. (1997)). Here we briefly review its main features regarding star formation and feedback.

The ISM is modeled as a two phase medium: hot gas ($`T_\mathrm{h}>10^4`$ K) and a cold phase ($`T_\mathrm{c}=10^4`$K), which pretends to mimic the molecular clouds from which stars are formed.

Cold clouds can be formed from the hot gas by radiative and Compton cooling or by thermal instability. Radiative cooling rates, $`\mathrm{\Lambda }_\mathrm{r}(T,Z)`$, are taken from Sutherland & Dopita (1993) for two metallicity cases: either primordial composition or solar, depending on the star formation history of each computational element.

Stars are created, from the cold phase, as collisionless particles at each timestep. Their mass is proportional to the density of cold gas in each cell:

$$\dot{\rho }_{}=\dot{\rho }_\mathrm{c}=(1\beta )\frac{\rho _\mathrm{c}}{t_{}}$$
(5)
where $`\beta `$ is the fraction of massive stars ($`M>10M_{}`$) that explode as supernovae. For a Salpeter IMF, $`\beta =0.12`$. For the characteristic time of star formation $`t_{}`$ we assume a constant value of $`10^8`$ yrs.

Several feedback mechanisms are included: photoionisation by UV light, gas heating and evaporation by supernova explosions. Massive stars are assumed to create a homogeneously distributed background of ionizing radiation, and gas should be dense enough to screen it in order to be able to cool into the cold phase and form stars (e.g. Giroux & Shapiro (1986); Mücket & Kates (1997)). A simple prescription, without resorting to a detailed computation of photoionization heating is to impose an overdensity threshold $`𝒟`$ for star formation to take place. On the other hand, supernova explosions can evaporate the cold gas phase and increase the temperature of the hot component. In our model, the supernova feedback parameter $`A`$ controls the amount of energy (relative to the supernovae mass) that goes into cold cloud evaporation:

$$\left(\dot{\rho }_\mathrm{h}\right)_{\mathrm{SN}}=\left(\dot{\rho }_\mathrm{c}\right)_{\mathrm{SN}}=\frac{\beta \rho _\mathrm{c}}{t_{}}A$$
(6)
$$\left(\dot{E}_\mathrm{h}\right)_{\mathrm{SN}}=\frac{\beta \rho _\mathrm{c}}{t_{}}\left(ϵ_{\mathrm{SN}}+Aϵ_\mathrm{c}\right)$$
(7)
where the subscripts $`h`$ and $`c`$ denote hot and cold gas, $`ϵ_{\mathrm{SN}}4.45\times 10^{49}`$ erg M$`{}_{}{}^{1}{}_{}{}^{}`$ is the supernovae specific energy (Salpeter IMF) and $`ϵ_\mathrm{c}`$ is the specific energy of the cold gas phase ($`T_\mathrm{c}=10^4`$ K).

### 3.3 Numerical experiments

Table 2 summarizes our set of experiments. In order to improve statistics, several realizations were run when the number of particles was low enough ($`N=128^3`$), changing the random seeds used to generate the initial conditions. The number of realizations is given by $`n_\mathrm{r}`$ in the second column. The cosmological model and the normalization $`\sigma _8`$ are shown in the next 5 columns, as well as the number of particles, box size, cell size (which determines our spatial resolution) and mass of each dark matter particle (mass resolution). In the last two columns we present the values of the supernova feedback parameter $`A`$ and the overdensity threshold $`𝒟`$ required for star formation.

The comparison of experiments $`\mathrm{\Lambda }`$CDM5, SCDM and BSI allows us to determine the effect of cosmology and feedback from supernovae on the global SFR history. Effects of photoionization are discussed from the results of experiment $`\mathrm{\Lambda }`$CDM4, and the dependence on resolution and volume can be studied from the whole set of $`\mathrm{\Lambda }`$CDM experiments. In $`\mathrm{\Lambda }`$CDM 1 and 2, the volume is high enough to include a cluster of galaxies, so we have also investigated the behavior of the SFR in different environments.

### 3.4 Resolution effects

In pure $`N`$-body simulations the lack of mass and force resolution leads to the well known overmerging problem (see e. g. Klypin et al. (1999); van Kampen (2001)), i. e. substructures of large halos or halos inside of halos are erased.

This effect becomes important in galaxy clusters. To avoid overmerging high force and mass resolution are necessary which may be reached later in adaptive mesh and multimass N-body/gasdynamical codes.

Our code is based on a fixed Eulerian mesh, so resolution is limited to the cell size. High spatial resolution can only be achieved by increasing the number of cells uniformly, or reducing the computational volume. As a trade off between the two, we have decided to restrict us on medium spatial resolution, in favour of large simulation volumes. Resolution effects are discussed in detail in the next section.

The numerical code has been parallelized using OpenMP compiler directives. Therefore, it runs very efficiently in parallel computers with either Shared or Non Uniform Memory Access (NUMA) architectures. The simulations have been performed on a variety of machines.

## 4 Results

### 4.1 Cosmology dependence

The global star formation histories from simulations of the three cosmological scenarios considered are compared in Fig. 2. The initial power spectrum normalization determines the evolution of structures, which represent a crucial factor in determining the rate of conversion of gas to stars.

Cosmological scenarios with CDM alone cannot explain structure formation on both small and very large scales. However, we have included SCDM as a reference model in our comparison. Due to its large power at small scales, SCDM produces many structures at high redshift. Therefore, we get a flat star formation rate even beyond $`z6`$, with a plateau about one decade over the observations (Fig. 2, middle panel). It does not seem very likely that this difference can be explained by dust absorption. In this model, too many baryons are transformed into stars too early with respect to observations, which poses an additional problem for the SCDM model.

On the contrary, the BSI model leads to a clear maximum in the SFR at $`z2`$, but the absolute value is too small and the decrease at higher redshifts is clearly too steep (Fig. 2, right). This model has not enough power at small scales for producing enough stars at high redshift. As a side effect, the cosmic star formation rate density remains too high at low redshifts, since a large part of cold gas remains in high density regions triggering further star formation activity.

A cosmological model with non-zero cosmological constant offers the best agreement with observations (Fig. 2, left). Both the overall trend and the height of the plateau fit reasonably well to the observational data points. At present, the $`\mathrm{\Lambda }`$CDM model fits best all known observational data (e.g. Bahcall et al. (1999)). Our results confirm that this model can also reproduce the observed star formation history of the universe when dust extinction is taken into account, in which case a significant amount of star formation is expected to take place at high redshift.

A further feature (common to all cosmologies) is the smooth redshift dependence of the mean SFR and the small range of scatter. This is in sharp contrast to the evolution of individual baryonic clumps, or the behavior of different regions in a single time-step of the simulation. There we observe huge bursts of star formation with an increase in activity by a factor of 10 and higher, as well as a quiescent regime with a star formation rate lower than the average.

In these simulations, the first stars are created around $`z10`$, which seems to be not too unrealistic, but it may depend on resolution. This was not a primary issue of the present study.

### 4.2 Feedback dependence

Fig. 3 shows the results of experiments $`\mathrm{\Lambda }`$CDM5 and $`\mathrm{\Lambda }`$CDM4, from where the effects of supernova feedback and photoionization were studied. The same trends described in this section were also seen in SCDM and BSI simulations.

Supernovae (left panel of Fig. 3) play a key role in the slope of the SFR density at low $`z`$. Gas heating and evaporation act as a self-regulating mechanism that inhibits subsequent gas infall when many new stars are formed. This is specially important in the most massive objects, which show a much higher star formation activity when feedback is not implemented. Supernova explosions can also blow away the gas reservoirs of dwarf galaxies. The effect on the global star formation rate is more noticeable at high redshift, when there are still very few galaxies massive enough to retain the heated gas within their gravitational potential wells. However, photoionization (right panel) is the dominant effect at very early epochs, preventing the gas from cooling and forming stars in low density regions (Efstathiou (1992)). As we have discussed above, our model takes into account this effect through the overdensity threshold parameter $`𝒟`$ for star formation. Only in regions denser than $`𝒟`$, molecular clouds will be able to form, and the onset of cosmic star formation activity is thus retarded until more and more objects become dense enough to screen the photoionizing background.

We found a SFR density consistent with observational data for $`𝒟100`$. This is in roughly agreement with results from a more detailed theoretical estimate of the same parameter (Mücket & Kates (1997)), although in that case the actual value of $`𝒟`$, as well as the temperature range in which thermal instability is an efficient process, depend on the local gas density and the intensity of the ionizing flux.

### 4.3 Environmental dependence

In the hierarchical scenario of structure formation, the first objects collapsed around the highest overdensity peaks and merged to eventually form the giant ellipticals in the centers of clusters and the bulges of the most massive isolated galaxies. These objects contain the vast majority of stars older than $`z=3`$, and hence their star formation history must be completely different from that of galaxies living in less dense environments.

The experiments $`\mathrm{\Lambda }`$CDM 1 and 2 sample a volume containing clusters, filaments, and groups of galaxies as well as isolated galaxies. Gottlöber et al. (2001) have studied the merging history of dark matter halos as function of halo environment. They found that halos located inside clusters have formed earlier than isolated halos of the same mass. Moreover, they showed that at higher redshifts ($`z14`$), progenitors of cluster and group galaxies have 3–5 times higher merger rates than isolated galaxies. Mergers of galaxies are thought to play a crucial role in the evolution of galaxies. In particular, the inflow of material may serve as a source of fresh gas and therefore increase the star formation rate.

In Fig. 4 we plot the contribution to the cosmic SFR density of progenitors of objects more massive than $`10^{13}`$ M at present (clusters, large groups) and less massive (isolated galaxies or small groups). The most massive objects have their peaks of star formation activity around $`z3`$, showing passive evolution since $`z2`$. This is in excellent agreement with the merging history found by Gottlöber et al. (2001). The progenitors of these objects were responsible for the bulk of star formation at high redshift, and can be very likely associated with the Lyman Break Galaxy population.

Isolated galaxies and galaxies in small groups are younger. They are responsible for most of the current stellar production. In these areas, halos massive enough to retain their gas content have a nearly constant SFR (with sporadic burst episodes).

The simulated volume in our numerical experiments is too small to be a representative cosmological volume. In fact, the larger volume $`\mathrm{\Lambda }`$CDM1 contains only one galaxy cluster of approximately $`3\times 10^{14}`$ M whereas the samller volume $`\mathrm{\Lambda }`$CDM2 contains four massive groups of about $`5\times 10^{13}`$ M. Comparing the solid lines in Fig. 4 we find qualitatively the expected scenario: the peak of star formation in the progenitors of the more massive cluster happened earlier (about $`z=4.5`$) than for the less massive groups ($`z=3`$).

In the absence of stellar feedback, most of the baryonic matter would have cooled into small mini-galaxies (Cole (1991); White & Frenk (1991); Blanchard et al. (1992)). Photoionization and supernova explosions inhibit cooling and star formation in very low-mass objects, so early star formation was biased towards the high overdensity regions that will collapse later to form the cores of the first clusters.

As the collapse proceeds, gas is heated to temperatures that make the cooling times too long and the SFR drastically drops once the cold gas has been exhausted.

At the same time, new galaxies will form and produce stars. Some of them will remain in rather isolated areas until present (dashed lines in Fig. 4) and some others will fall into the potential wells of clusters. It seems very likely that tidal interactions trigger intense star formation bursts in these infalling galaxies (Lavery & Henry (1988)), but ram pressure of the ICM may also deprive them from their cold gas reservoirs (Gunn & Gott (1972); Quilis et al. (2000)). This later mechanism is supported by the observational evidence of the displacement of the neutral hydrogen disk of post-starburst galaxies with respect to the stellar population (Valluri & Jog (1991)).

Cluster galaxies would then stop/decay their star formation activity and become redder, giving rise to the Butcher & Oemler (1978) effect. At the same time, they will experience a morphological transformation, due to mergers with their neighboring cluster galaxies, from spirals into S0 (Dressler et al. (1997)). At present, most of the observed clusters are at $`z0.5`$ (Fasano et al. (2000), and references therein), so deeper surveys are clearly needed in order to check this scenario.

### 4.4 Resolution and volume dependence

Star formation and feedback processes are incorporated in gasdynamical simulations (either SPH or Eulerian) by means of rather simple recipes which pretend to extrapolate the effects of local physics acting on scales of a few pc and masses of $`10^5`$ M below the resolution limit (often at kpc and higher). In our case, part of this extrapolation is hidden in the two parameters: The supernova feedback parameter $`A`$ and the overdensity threshold for star formation $`𝒟`$. Whether this parameterization is a reasonable extrapolation of the underlying subgrid physics depends on the stability of results against changes in spatial resolution of the simulations.

On the other hand, when one tries to estimate cosmological averaged quantities, like the SFR density, it is necessary to check for deviations of these quantities due to small number statistics. In the case that matters here, a volume averaged quantity, the cosmic variance can have a relatively large contribution in the particular determination of the comoving SFR density. Thus, we have simulated different computational volumes with similar resolutions, depending on the computational facilities.

To check for these two effects we compare in Fig. 5 our results for $`\mathrm{\Lambda }`$CDM simulations with different resolutions and volumes.

The most noticeable feature that we derive from this Figure is the steeper slope of the local SFR for larger volume simulations. As we increase the volume, larger structures are being simulated (i.e. clusters) and as explained above, the reduced star formation activity in the cores of clusters in recent epochs can explain this behavior. In simulations with smaller volumes, most of the galaxies that form are either isolated or in small groups, which have a completely different star formation activity. Thus, we can conclude from our simulations that the effects of cosmic variance are quite important in determining the proper behavior of the cosmic SFR density from $`z1`$ to the present.

We do not find significant differences in our estimates of the SFR density from simulations with different volumes and spatial resolutions, but there is a general trend to have more SFR density at early times for simulations with higher resolution. This is a consequence of the overcooling problem (see Balogh et al. (2001) for a recent review), since the number of resolved low mass halos increases in these simulations. A detailed analysis of this question would require a more realistic treatment of photoionization.

Keeping all these considerations in mind, we should be extremely cautious when making quantitative estimates of the cosmic SFR density at any given $`z`$. However, this quantity shows a similar behavior in all $`\mathrm{\Lambda }`$CDM experiments. This is a very promising result indicating that its qualitative evolution (particularly at high redshift) is a robust prediction of our simulations.

No characteristic epoch of cosmic star formation is found for a $`\mathrm{\Lambda }`$CDM universe. The comoving SFR density increases with lookback time until it reaches a nearly constant mean plateau beyond $`z=2`$. Stars probably formed gradually at high redshift, according to most recent observations and taking dust extinction into account. At some point around $`z1`$, the global star formation activity would have declined sharply until the present day, although individual SFRs of active galaxies are kept approximately constant.

The star density $`\rho _{}(z)`$ shows a similar behavior for both CDM and $`\mathrm{\Lambda }`$CDM models. The most remarkable difference is found in the present value of $`\rho _{}(0)`$ (see Table 3), which can be observationally deduced from the local luminosity density, assuming a constant mass-to-light ratio (Madau et al. (1996); Fukugita et al. (1998)):

$$\rho _{}^{\mathrm{obs}}(0)7\times 10^8h\mathrm{M}_{}\mathrm{Mpc}^3$$
(8)
Based on the fraction of stars in bulges, ellipticals and S0s, Renzini (1998) estimated that about $`30\%`$ of the stellar content of the universe should have already been formed by $`z=3`$, regardless of the cosmological model. The amount of stars produced in experiments $`\mathrm{\Lambda }`$CDM1 and 2 (Table 3), both at $`z=3`$ and $`z=0`$, are in fair agreement with these observationally derived quantities. For the other experiments, the fraction of stars older than $`z=3`$ is significantly lower, which is not only due to a low star formation activity at high redshift, but also to a higher SFR in the local universe and the effects of cosmic variance due to the smaller simulated volumes.

## 5 Conclusions

We have performed 66 numerical simulations, grouped in 7 different experiments, in order to obtain a theoretical prediction of the cosmic star formation history. We compared our predictions with a number of published observational estimates. We summarize our results as follows:

1. $`\mathrm{\Lambda }`$CDM shows the best agreement with observations. Standard CDM tends to overpredict the cosmic SFR density, while BSI fails to form enough stars at high $`z`$ due to its reduced power at small scales.
2. Photoionization fixes the onset of star formation in the universe, supressing star formation in low-mass objects at high $`z`$. Supernovae feedbacks are very efficient to self-regulate the process of converting gas into stars, preventing the most massive halos to form too many stars at later epochs.
3. Star formation histories are very different for objects that end up in the cores of clusters and isolated galaxies. Galaxies in cluster cores have much older stellar populations, and their activity is drastically reduced when they lose their cold gas reservoirs.
4. In contrast to the first observational estimates (Madau et al. (1996)), but in agreement with most recent data (Steidel et al. (1999)), the $`\mathrm{\Lambda }`$CDM model predicts almost no drop in the cosmic SFR for $`2<z<5`$. Star formation seems to be a gradual process with no characteristic epoch.
5. The comoving simulated volumes in our experiments are comparable to those covered by observations. Thus, one should expect a non negligible error in the observational measurements due to cosmic variance associated with small volume statistics. This is most important when one wants to determine the steep slope of the SFR density from $`z1`$ to the present time.

Selection effects, dust extinction and the conversion from luminosities to SFRs are the major sources of uncertainty in the observational estimates. The small volume sampled by the HDF may also introduce some statistical bias at high redshift.

Uncertainty in our numerical experiments comes from spatial resolution and the star formation and feedback prescriptions. The use of Adaptive Mesh Refinement Eulerian codes (AMR) will be useful in order to improve resolution in statistically significant cosmological volumes. Photoionization and chemical enrichment have been given a phenomenological treatment in our code. It would be extremely interesting to self-consistently model both processes, computing the intensity of the ionizing background at each time-step, as well as metal production and advection.

###### Acknowledgements.

This work has been partially supported by the SEUID (Spain) under project number PB96-0029, by the Acciones Integradas Hispano-Alemanas HA2000-0026 and by Deutscher Akademischer Austauschdienst DAAD. We thank the Centro Europeo de Paralelismo de Barcelona and the Centro de Computación Científica de la Universidad Autónoma de Madrid for allowing us to use their computational facilities in which most of the simulations reported in this paper were performed.