format: text

text:
# Formation of the First Stars by Accretion

## 1 Introduction

The problem of the characteristic mass of the first stars that formed in the universe is still largely unsolved. Recent theoretical and numerical studies are revealing in leaps and bounds the physical state of star forming regions within the first cosmic structures. In the framework of cold dark matter (CDM) cosmology smaller objects tend to collapse earlier. In these models, a moderately rare object of $`3\sigma `$ overdensity collapses and virializes at about $`z30`$. However, only sufficiently massive objects can cool in a Hubble time and become luminous by forming stars. The minimum cooling mass at $`z30`$ is $`M_{\mathrm{min}}10^6M_{\mathrm{}}`$, and is a decreasing function of the collapse redshift (e.g., Haiman, Thoul, & Loeb 1996; Tegmark et al. 1997; Fuller & Couchman 2000). After virialization, objects more massive than $`M_{\mathrm{min}}`$, called “first objects”, can collapse gravitationally owing to molecular hydrogen cooling (e.g., Palla, Salpeter, & Stahler 1983, hereafter PSS), and fragment into high-density clumps. According to numerical studies, the mass scale of these clumps is rather high, of the order of $`10^3M_{\mathrm{}}`$ (Abel, Bryan, & Norman 2000, 2002; Bromm, Coppi, & Larson 1999, 2002; Tsuribe & Inutsuka 2001). However, the formation of fragments of lower mass, down to $`1M_{\mathrm{}}`$, is also possible (e.g., Yoshii & Saio 1986; Uehara et al. 1996; Nakamura & Umemura 1999, 2001, 2002; Uehara & Inutsuka 2000; Omukai 2001). Indeed, the recent discovery of HE0107-5240, a star of 0.8 $`M_{\mathrm{}}`$ in the galactic halo with virtually no metals (\[Fe/H\]=$``$5.3), attests the ability of very metal-poor gas clouds to give birth to low-mass stars (Christlieb et al. 2002).

After fragmentation, the evolution of gravitationally unstable clumps proceeds in a highly non-homologous fashion, with the central parts collapsing first. This runaway phase is induced by cooling provided by H2 line radiation at densities $`n10^{14}\mathrm{cm}^3`$, and by H2 collision-induced emission at higher densities. The resulting gas temperature is nearly constant at several $`10^2`$ K, and the innermost region of $`1M_{\mathrm{}}`$ becomes fully molecular due to the three-body reaction (PSS). At densities $`n10^{16}\mathrm{cm}^3`$, the cloud becomes optically thick to collision-induced absorption and H2 dissociation works as an effective cooling agent. Finally, without further cooling mechanism, at $`n10^{22}\mathrm{cm}^3`$, a small hydrostatic core of mass $`10^3M_{\mathrm{}}`$ is formed (Omukai & Nishi 1998), similar to the conditions found in studies of the collapse of clouds of standard solar composition, i.e., the stellar (or second) core according to Larson (1969).

In the primordial clouds described by the numerical simulations, the outcome of the runaway collapse is a tiny ($`10^3M_{\mathrm{}}`$) protostar, surrounded by a large amount of reservoir gas ($`10^3M_{\mathrm{}}`$). As a result, the protostar can grow by several orders of magnitude in mass by accreting the envelope matter. The mass accretion rate onto the protostar is determined by the radial density distribution at the time of core formation, and is related to the prestellar temperature or, equivalently, to the effective sound speed $`c_\mathrm{s}`$ by (e.g., Stahler, Shu, & Taam 1980):

$$\dot{M}_{\mathrm{acc}}\frac{c_\mathrm{s}^3}{G}.$$
(1)
The evolution of zero-metal stars in the main accretion phase was studied by Stahler, Palla, & Salpeter (1986; hereafter SPS) up to 10.5 M and by Omukai & Palla (2001; hereafter, Paper I) for higher masses under a constant mass accretion rate of $`4.41\times 10^3M_{\mathrm{}}`$ yr-1, corresponding to a gas temperature of 1700 K in eq. (1).

The mass of the forming star is set at the time when protostellar accretion stops. In paper I, we have found that fast accreting protostars enter a phase of rapid expansion at a mass of $`300M_{\mathrm{}}`$ when the luminosity becomes very close to the Eddington limit. This event may determine the onset of a powerful stellar wind driven by radiation pressure that can effectively quench further accretion of circumstellar material. In that event, after a transient internal readjustement, the protostar is expected to settle as an ordinary zero-age main sequence (ZAMS) star.

The upper limit of $`300M_{\mathrm{}}`$ to the protostellar mass has been obtained using a constant value of the mass accretion rate. In realistic collapse calculations, deviations from the isothermal assumption are expected to occur, and even in the ideal case the mass accretion rate is found to vary in time (e.g., Foster & Chevalier 1993; Tomisaka 1996). Therefore, it is important to extend our initial study to a larger set of conditions in order to explore the robustness of the results described above. Along with differences in the mass accretion rate, we will also address the issue of the impact of a small abundance of heavy elements. It is expected that metals may alter the mass limit drastically as a result of the existence of a dust shell and because of steller evolution effects.

In §2, we sketch our numerical approch. Then, in §3 we describe the results of the main accretion phase obtained with a constant and time variable $`\dot{M}_{\mathrm{acc}}`$. The effects of a finite, non-zero metallicity on protostellar evolution are discussed in §4. In §5 we discuss the implications of our work. Finally, in §6, we give the conclusions.

## 2 Numerical Approach

The basic strategy and equations used here are the same as in SPS. In that scheme, the protostellar evolution is treated as a sequence of steady state accretion flows onto a growing core (see Figure 1 of SPS for a definition of the various regions). The core is assumed to be in hydrostatic equilibrium and the ordinary stellar structure equations are applied. If the gas in front of the accretion shock is optically thin, no envelope model is constructed and the boundary condition given by eq. (4a) of SPS is applied at the core surface. In case the preshock gas is optically thick, we integrate the equations inside the radiative precursor from the photosphere to the core surface (see eqs. 7a-d of SPS). Outside the photosphere, we assume a free-falling flow. The main difference with the SPS work is in the use of updated opacity tables: for $`T<7000`$ K, we take the results of Lenzuni, Chernoff, & Salpeter (1991) for a gas composition $`X=0.72,Y=0.28`$, and the OPAL opacity at higher temperatures (Iglesias & Rogers 1996), with a slightly different hydrogen abundance $`X=0.70`$). For partially polluted gas, we have used the results of Alexander & Ferguson (1994) with $`X=0.70,Y=0.30`$ for $`T<7000`$ K, and the OPAL opacity for higher temperatures. Although the primordial He abundance is about $`Y=0.24`$ (e.g., Galli et al. 1995; Izotov & Thuan 1998), we used rather high value of $`Y0.28`$ as stated above because of the lack of opacity tables for primordial gas. This causes no significant error.

The calculations are started at the protostellar mass of $`0.01M_{\mathrm{}}`$. The initial models are arbitrary and constructed using an entropy distribution of the core given by

$$s(M)=s_0+\beta (M/0.01M_{\mathrm{}})^2,$$
(2)
where $`M`$ is the enclosed mass within a spherical surface, $`s_0=13.05`$ and $`\beta =7.488`$ in units of $``$, and the zero-point is arbitrarily set at $`T=2.05\times 10^5`$ K and $`\rho =5.16`$ g cm-3. With the mass accretion rate used by SPS, $`\dot{M}_{\mathrm{acc}}=4.41\times 10^3`$ M yr-1, the initial central temperature is $`10^5`$ K. This value increase (decreases) for higher (lower) accretion rates. The initial conditions do not affect the later evolution, since the protostar adjusts quickly during the “decay-of-transient phase” (for the protostellar mass $`M_{}<0.1M_{\mathrm{}}`$; SPS) to the state appropriate to an accreting object.

## 3 Formation of Primordial Stars:
Effects of Variation of the Mass Accretion Rate

### 3.1 Evolution under different mass accretion rates

First, we study the effects of using different, but constant values of $`\dot{M}_{\mathrm{acc}}`$ on the evolution of primordial protostars. We define $`\dot{M}_{\mathrm{fid}}=4.41\times 10^3M_{\mathrm{}}`$ yr-1 as the fiducial value of the accretion rate used by SPS and in Paper I. To follow the dependence of the evolution on $`\dot{M}_{\mathrm{acc}}`$, we have studied four cases, varying $`\dot{M}_{\mathrm{acc}}`$ from $`1/4\dot{M}_{\mathrm{fid}}`$ to $`2\dot{M}_{\mathrm{fid}}`$, thus covering the range $`1.18.8\times 10^3M_{\mathrm{}}`$ yr-1.

The evolution of the protostellar radius $`R_{}`$ as a function of the accretion rate is shown in Figure Formation of the First Stars by Accretion. As the calculations start from rather arbitrary initial conditions, the protostar adjusts itself to the state appropriate to an accreting object by entropy redistribution. As a result, the fast entropy variations cause the bumps in radius seen in Figure Formation of the First Stars by Accretion at masses $``$0.03 $`M_{\mathrm{}}`$. After this minor transient phase, the evolution proceeds smoothly and the similarity of the curves allows the identification of three distinct phases. The first two are characterized by a common trend where the protostellar radius $`R_{}`$ increases almost linearly with the protostellar mass $`M_{}`$ in the logarithmic scale, reaches a peak between $``$10 and 20 $`M_{\mathrm{}}`$, and then drops substantially. The adiabatic accretion phase (Phase I, $`M_{}`$ a few to $`10M_{\mathrm{}}`$) is marked by a gradual expansion of the radius and lasts until the arrival of the internal luminosity wave at the surface, whose eruption causes the sudden swell of the external layers (e.g., Palla & Stahler 1991). In the Kelvin-Helmholtz contraction phase (Phase II, a few-$`10M_{\mathrm{}}M_{}3060M_{\mathrm{}}`$) the increased gravitational pull of the growing star stops the expansion and forces the star to contract toward the conditions appropriate for a main sequence object. The third phase varies markedly, depending critically on the adopted value of $`\dot{M}_{\mathrm{acc}}`$.

If the rate is comparable to $`\dot{M}_{\mathrm{fid}}`$ radiation pressure causes the abrupt expansion of the surface layers (Phase IIIb) seen in the two upper curves of Figure Formation of the First Stars by Accretion computed for $`\dot{M}_{\mathrm{acc}}=(12)\dot{M}_{\mathrm{fid}}`$. The protostellar radius increases tremendously and the impulsive release of energy likely causes the stripping of the external layers of the protostar by a radiation-driven stellar wind. As shown in Paper I, this event signals the end of the accretion phase and sets the maximum value of the protostellar mass. In the opposite case where $`\dot{M}_{\mathrm{acc}}<\dot{M}_{\mathrm{fid}}`$, contraction proceeds unimpeded by accretion. In the cases shown in Figure Formation of the First Stars by Accretion and computed for $`\dot{M}_{\mathrm{acc}}=1/2`$ and $`1/4\dot{M}_{\mathrm{fid}}`$, the protostar reaches the ZAMS at a mass of $`50M_{\mathrm{}}`$ (Phase IIIa).

The global evolution of the core interior is shown in Figure Formation of the First Stars by Accretion for the low mass accretion rate ($`\dot{M}_{\mathrm{acc}}=1/4\dot{M}_{\mathrm{fid}}`$) and in Figure Formation of the First Stars by Accretion for the case $`\dot{M}_{\mathrm{acc}}=\dot{M}_{\mathrm{fid}}`$. In these figures, the bottom panels illustrate the evolution of the central variables, as well as those of the position of the temperature maximum for the same cases.

The overall evolution of the interior luminosity $`L_{}`$ is shown in the top panel for $`\dot{M}_{\mathrm{acc}}=1/4\dot{M}_{\mathrm{fid}}`$ and in the bottom panel of Figure Formation of the First Stars by Accretion for $`\dot{M}_{\mathrm{acc}}=2\dot{M}_{\mathrm{fid}}`$. The contribution to the total luminosity provided by nuclear reactions of deuterium and hydrogen (through the $`pp`$-chain and CN-cycle) at the center of the protostar is also indicated. Notice how the onset of each burnig stage is shifted to higher masses for bigger values of $`\dot{M}_{\mathrm{acc}}`$.

Having sketched the main features of the evolution, we now describe in detail the results of each evolutionary phase identified above.

#### 3.1.1 Phase I: adiabatic accretion and propagation of the luminosity wave

Early in the evolution, the opacity in the interior is due to free-free absorption ($`\kappa _{\mathrm{ff}}\rho T^{3.5}`$) and is very high due to the low temperatures ($`T<10^6`$ K). As shown in Figure Formation of the First Stars by Accretion, the luminosity is extremely low and consequently the cooling time of the internal regions (i.e., the Kelvin-Helmholtz time $`t_{\mathrm{KH}}=GM_{}^2/R_{}L_{}`$) is longer than the evolutionary timescale (i.e., the accretion time $`t_{\mathrm{acc}}=M_{}/\dot{M}_{\mathrm{acc}}`$). As a consequence, the accreted material piles up without further cooling, after passing through the accretion shock and the settling region. During this adiabatic accretion phase, the protostar swells gradually as the mass increases. The protostellar radius obeys the relation (see equation 21 of SPS)

$$R_{}M_{}^{0.27}\dot{M}_{\mathrm{acc}}^{0.41},$$
(3)
that can be derived assuming that the luminosity is dominated by the accretion contribution and that the opacity in the envelope is due to H- bound-free absorption, scaling as $`\kappa T^{14.5}`$.

The increase of the temperature accompanied by the growth of the protostellar mass produces a significant decrease of the opacity. Eventually, the heat contained inside the protostar begins to propagate outward as a luminosity wave (see Figure 10 of SPS). Then, the interior luminosity $`L_{}`$ increases dramatically at a mass of $`4M_{\mathrm{}}`$ for $`\dot{M}_{\mathrm{acc}}=1/4\dot{M}_{\mathrm{fid}}`$, and at $`8M_{\mathrm{}}`$ for $`\dot{M}_{\mathrm{acc}}=\dot{M}_{\mathrm{fid}}`$ (see Figure Formation of the First Stars by Accretion). The arrival of the luminosity wave at the surface produces the sudden swelling of the protostar, marking the end of Phase I. In the course of this expansion, a short-lived surface convective region develops which extends to about 20% in mass (Figs. Formation of the First Stars by Accretion, Formation of the First Stars by Accretion).

In the case of low $`\dot{M}_{\mathrm{acc}}`$ (Figure Formation of the First Stars by Accretion), we can see a temporary, but active phase of deuterium burning around $`3M_{\mathrm{}}`$, when the nuclear energy generation rate exceeds $`0.1L_{}/M_{}`$. This episode is just the consequence of the low value of $`L_{}`$ before the arrival of the luminosity wave (see the top panel of Figure Formation of the First Stars by Accretion). However, the deuterium luminosity $`L_\mathrm{D}`$ is very small and deuterium is hardly consumed with no significant consequences on the evolution.

As a result of the outward increase of the entropy, the temperature is highest off-center, at an internal mass of $`M/M_{}0.1`$. The position of the off-center temperature peak is shown in the top panels of Figs. Formation of the First Stars by Accretion and Formation of the First Stars by Accretion and is maintained throughout the adiabatic accretion phase.

#### 3.1.2 Phase II: Kelvin-Helmholtz contraction

After the relaxation of the interior and the passage of the luminosity wave, the major opacity source in protostars is provided by electron scattering, except in a small central region where free-free absorption is dominant, and in a thin surface layer characterized by bound-free absorption of H and He.

The condition of radiative equilibrium and of an opacity independent of density and temperature throughout most of the interior determines the well known mass–luminosity relation $`L_{}M_{}^3`$, and the Kelvin-Helmholtz time becomes

$$t_{\mathrm{KH}}=\frac{GM_{}^2}{R_{}L_{}}\frac{1}{M_{}R_{}}.$$
(4)
If the radius is sufficiently large, $`t_{\mathrm{KH}}<t_{\mathrm{acc}}`$ and the core shrinks due to fast energy loss, while in the opposite case the core swells adiabatically. The radius of the protostar actually adjusts towards equality of the two timescales $`t_{\mathrm{KH}}t_{\mathrm{acc}}`$, leading to the relation

$$R_{}\dot{M}_{\mathrm{acc}}/M_{}^2.$$
(5)
Small deviations from this relation are caused by minor contributions to the opacity from processes other than electron scattering.

Figure Formation of the First Stars by Accretion shows the evolution of the interior luminosity to mass ratio, $`L_{}/M_{}`$, as a function of the protostellar mass $`M_{}`$. During KH contraction, $`L_{}/M_{}`$ increases rapidly and the mass-luminosity relation follows approximately a unique relation independent of the value of the mass accretion rate, as expected in the ideal case of constant opacity. Also, the photospheric luminosity, which is the sum of the interior ($`L_{}`$) and accretion ($`L_{\mathrm{acc}}`$) luminosity, is independent of $`\dot{M}_{\mathrm{acc}}`$, since $`L_{\mathrm{acc}}\dot{M}_{\mathrm{acc}}/R_{}`$ and $`R_{}\dot{M}_{\mathrm{acc}}`$.

At the beginning of contraction, the maximum temperature exceeds $`10^6`$K, and deuterium can start burning. Because of the off-center temperature peak, burning propagates both inward and outward, scouring out all of the available deuterium (see Fig.s Formation of the First Stars by Accretion and Formation of the First Stars by Accretion). Under the high accretion rates considered here, however, D-burning plays a minor role in the evolution and contributes at most $``$1/4 of the interior luminosity (see Figure Formation of the First Stars by Accretion). The off-center temperature peak is only a transient property of the core and disappears at $`20M_{\mathrm{}}`$ for $`\dot{M}_{\mathrm{acc}}=1/4\dot{M}_{\mathrm{fid}}`$, and $`40M_{\mathrm{}}`$ for $`\dot{M}_{\mathrm{acc}}=\dot{M}_{\mathrm{fid}}`$, as seen in Figs. Formation of the First Stars by Accretion and Formation of the First Stars by Accretion.

After deuterium is completely destroyed and the protostar has contracted further, the $`pp`$-chain becomes operative in the center. However, even in this case, the energy generated by these reactions is not enough to contrast contraction effectively, owing to the saturation of the energy generation rate at high temperatures (see Figure Formation of the First Stars by Accretion). Therefore, despite nuclear heating, protostars can continue the phase of gravitational contraction.

#### 3.1.3 Phase IIIa: Settling onto the ZAMS with continuing accretion ($`\dot{M}_{\mathrm{acc}}1/2\dot{M}_{\mathrm{fid}}`$)

For accretion rates of $`\dot{M}_{\mathrm{acc}}=1/2`$ and 1/4$`\dot{M}_{\mathrm{fid}}`$, during contraction the central temperature becomes high enough to synthesize carbon via the CN-cycle and H-burning becomes a significant source of luminosity. Figure Formation of the First Stars by Accretion shows that at this time gravitational contraction is halted and the interior luminosity to mass ratio $`L_{}/M_{}`$ stops increasing. The protostar relaxes quickly to the ZAMS and the models converge toward a single mass-luminosity-radius relation, as shown in Fig.s Formation of the First Stars by Accretion and Formation of the First Stars by Accretion. With H-burning working as a thermostat, the central temperature remains almost constant at about $`10^8`$ K (see the bottom panel of Figure Formation of the First Stars by Accretion).

At this point, the interior luminosity is essentially supplied by nuclear reactions (see the top panel of Figure Formation of the First Stars by Accretion). Since part of this energy source is used to heat the interior, the specific entropy in the convective core increases as the stellar mass grows. This results in the decrease of the central density, while the temperature remains constant. The mass fraction of the convective core increases (see the top panel of Figure Formation of the First Stars by Accretion) and the CN abundance reaches very quickly a level of $`10^9`$, thereby gradually increasing (see the bottom panel of Figure Formation of the First Stars by Accretion).

As shown in Figure Formation of the First Stars by Accretion, throughout the evolution the luminosity remains below the Eddington value by a rather large margin. The latter has been computed using electron scattering as the main source of opacity in the accretion flow. As a result, accretion can continue without experiencing the effects of the retardation force exerted onto the ionized radiative precursor. Therefore, the upper limit of the protostellar mass formed in these conditions is set by the amount of circumstellar material left during accretion and not by the feedback effect of radiation pressure.

#### 3.1.4 Phase IIIb: Rapid expansion due to radiation pressure ($`\dot{M}_{\mathrm{acc}}\dot{M}_{\mathrm{fid}}`$)

For high accretion rates, the protostar undergoes a phase of violent expansion as soon as the interior luminosity becomes comparable to the Eddington limit. This process is caused by the strong radiation force exerted both on the stellar surface and the radiative precursor. The result is that fast accreting protostars never converge to an ordinary ZAMS structure. Let us see more in detail how this condition is achieved.

During the contraction phase, the luminosity-to-mass ratio increases until the conditions at the center are appropriate to start nuclear burning (see Fig. Formation of the First Stars by Accretion). This happens at progressively large masses as the accretion rate increases. For sufficiently high $`\dot{M}_{\mathrm{acc}}`$, the total luminosity, which is the sum of the interior and accretion luminosity, reaches the Eddington limit before the ignition of hydrogen. This is the case for $`\dot{M}_{\mathrm{acc}}=2\dot{M}_{\mathrm{fid}}`$, while for the lower value of $`\dot{M}_{\mathrm{acc}}=\dot{M}_{\mathrm{fid}}`$ the situation is more complicated. In the accreting envelope, gravitational acceleration is approximately balanced by the deceleration due to the radiation force. Then, the ram pressure onto the protostellar surface decreases suddenly, and the protostar begins to expand without violating the Eddington limit by lowering the accretion luminosity $`L_{\mathrm{acc}}R_{}^1`$. During expansion, the temperature drops, while the radiation force and the opacity increase at the surface, the latter peaking sharply owing to the contribution of ionization of trace amount of atoms. This is mostly a surface effect since the main opacity source is still electron scattering both inside the protostar and in the accreting envelope. By this mechanism, the expansion is accelerated and possibly results in the stripping of the surface layers and in the reversal of the material in the accreting envelope. Thus, the main phase of accretion is over and the protostar can rapidly relax onto a ZAMS structure. According to the results shown in Fig. Formation of the First Stars by Accretion, this critical phase is achieved at a mass of $`90M_{\mathrm{}}`$ for $`\dot{M}_{\mathrm{acc}}=2\dot{M}_{\mathrm{fid}}`$, and $`300M_{\mathrm{}}`$ for $`\dot{M}_{\mathrm{acc}}=\dot{M}_{\mathrm{fid}}`$. Therefore, our models show that the upper limit of the mass shifts toward smaller values for progressively higher values of the mass accretion rate.

#### 3.1.5 A general consideration

We have seen that, depending on the actual value of the accretion rate, the final fate of primordial protostars can be quite different and leads to strongly different values of the maximum limit of their mass. This suggests to look deeper into the physical origin of the bifurcation displayed in the mass-radius relation of Figure Formation of the First Stars by Accretion.

The total luminosity of a protostar, $`L_{\mathrm{tot}}`$, consists of the interior, $`L_{}`$, and accretion luminosity, $`L_{\mathrm{acc}}`$:

$$L_{\mathrm{tot}}=L_{}+L_{\mathrm{acc}}L_{}+GM_{}\dot{M}_{\mathrm{acc}}/R_{}.$$
(6)
Imagine that the protostar reaches the ZAMS star while continuing accretion, as in the cases of $`\dot{M}_{\mathrm{acc}}1/2\dot{M}_{\mathrm{fid}}`$. Since a unique mass-luminosity-radius relation holds for ZAMS stars, the total luminosity would be

$$L_{\mathrm{tot}}L_{\mathrm{ZAMS}}+GM_{}\dot{M}_{\mathrm{acc}}/R_{\mathrm{ZAMS}},$$
(7)
where we have set $`R_{}=R_{\mathrm{ZAMS}}(M_{})`$ and $`L_{}=L_{\mathrm{ZAMS}}(M_{})`$. If the accretion rate is larger than the critical value

$$\dot{M}_{\mathrm{crit}}=\frac{4\pi cR_{\mathrm{ZAMS}}}{\kappa _{\mathrm{es}}}\left(1\frac{L_{\mathrm{ZAMS}}}{L_{\mathrm{Edd}}}\right)4\times 10^3M_{\mathrm{}}/\mathrm{yr},$$
(8)
where $`\kappa _{\mathrm{es}}`$ is the electron scattering opacity, the total luminosity exceeds the Eddington limit $`L_{\mathrm{Edd}}`$, implying that the star cannot adjust to the ZAMS structure at these accretion rates. Although the critical accretion rate is a function of the mass of the central star, the dependence is weak considering the $`M_{}L_{}R_{}`$ relations.

Note that the critical accretion rate $`\dot{M}_{\mathrm{crit}}`$ is the value where the total luminosity $`L_{\mathrm{tot}}`$ reaches $`L_{\mathrm{Edd}}`$ at the same time as the protostar reaches the ZAMS. Thus, we can derive $`\dot{M}_{\mathrm{crit}}`$ by using the ZAMS quantities, although protostars start expanding before reaching the ZAMS for $`\dot{M}_{\mathrm{acc}}>\dot{M}_{\mathrm{crit}}`$. For $`\dot{M}_{\mathrm{acc}}<\dot{M}_{\mathrm{crit}}`$, $`L_{\mathrm{tot}}`$ remains below $`L_{\mathrm{Edd}}`$ up to the point where the protostar reaches the ZAMS, while for $`\dot{M}_{\mathrm{acc}}>\dot{M}_{\mathrm{crit}}`$, $`L_{\mathrm{tot}}`$ equals $`L_{\mathrm{Edd}}`$ in the course of the Kelvin-Helmholtz contraction and the protostar starts expanding.

The weak dependence of $`\dot{M}_{\mathrm{crit}}`$ on the protostellar mass $`M_{}`$ can be seen more clearly if we define the Eddington radius as the radius of the core where the total luminosity (defined in eq. 7) equals the Eddington limit $`L_{\mathrm{Edd}}`$:

$$R_{\mathrm{Edd}}=\frac{\kappa _{\mathrm{es}}\dot{M}_{\mathrm{acc}}}{4\pi c(1\mathrm{\Gamma })},$$
(9)
where $`\mathrm{\Gamma }L_{}/L_{\mathrm{Edd}}`$.

In the two panels of Figure Formation of the First Stars by Accretion, we show the Eddington radius and the core radius for the four values of the accretion rate. In all cases, when the core mass is small and $`\mathrm{\Gamma }1`$, the Eddington radius is constant. Later on, the evolution differs. Consider first the case for low $`\dot{M}_{\mathrm{acc}}`$. During the contraction phase, the luminosity-to-mass ratio increases, the core shrinks, but the Eddington radius starts a slow rise. However, $`R_{\mathrm{Edd}}`$ never climbs up to $`R_{}`$, since the core meets the ZAMS conditions before and its radius also increases as a result of the heating provided by nuclear burning. Thereafter, $`R_{\mathrm{Edd}}`$ remains below the ZAMS radius and accretion can proceed unimpeded by radiation forces. Since the $`R_{\mathrm{Edd}}/R_{}`$ ratio remains approximately constant, the critical accretion rate is almost independent of the protostellar mass.

If, on the other hand, the protostar reached the ZAMS structure with an accretion flow $`\dot{M}_{\mathrm{acc}}>\dot{M}_{\mathrm{crit}}`$, the total luminosity would exceed the Eddington limit. Unlike the previous case, at some point of the evolution the Eddington and core radii would become equal (see the bottom panel of Fig. Formation of the First Stars by Accretion). The core now undergoes the runaway expansion, keeping the luminosity slightly below the Eddington limit. The increased radius reduces the accretion luminosity, thus allowing more gas to flow for a short period.

#### 3.1.6 Properties of the radiative precursor

Except for the short-lived phase of the peak radius around $`10M_{\mathrm{}}`$, primordial protostars are embedded within an optically thick envelope. The spatial extent of the precursor, which is the optically thick part of the accretion flow, corresponds approximately to that of the ionized region. However, the ionized region around the accreting protostar differs from an ordinary HII region in that it is opaque to optical photons, and both matter and radiation are in the thermodynamic equilibrium.

In Figure Formation of the First Stars by Accretion, we show the evolution of the photospheric radius for two cases with $`\dot{M}_{\mathrm{acc}}=1/4\dot{M}_{\mathrm{fid}}`$ and $`\dot{M}_{\mathrm{fid}}`$. During the adiabatic phase, most of the luminosity is generated at the accretion shock and the ratio of the photospheric to core radius remains nearly constant, with a value of $``$1.4 (see also eq. 19 of SPS). After the luminosity wave reaches the surface, the interior luminosity exceeds the accretion contribution. As a result of core contraction and of the increased heat deposit in the precursor, the ratio of photospehric to core radius increases dramatically. The trend is reversed during the final stages of the evolution.

The evolutionary tracks of the photosphere in the HR diagram are displayed in Figure Formation of the First Stars by Accretion for the four accretion histories. For comparison, we also show the evolution of the protostellar surface and locus of the ZAMS for metal-free stars. The distinction between the surface temperature and the ZAMS effective temperature after reaching the MS phase is a consequence of the difference in boundary conditions at the stellar surface. Because of the optically thick envelope, the hot surface of the protostar remains invisible in the optical. The photospheric temperature is locked at $``$6000 K due to strong sensitivity of the H- bound-free opacity on temperature. However, for low $`\dot{M}_{\mathrm{acc}}`$, the tracks make an abrupt shift toward higher temperatures for $`M_{}100M_{\mathrm{}}`$. As the protostellar mass increases, the density in the accreting envelope decreases, while the temperature and ionization degree both increase. Electron scattering replaces the H- b-f absorption as the major source of opacity at the photosphere. Since this process is independent of temperature, the effective temperature rises as the density in the accreting envelope decreases. Correspondingly, the photospheric radius decreases sharply (Figure Formation of the First Stars by Accretion). The tracks would then join the corresponding ZAMS models at masses exceeding $`1000M_{\mathrm{}}`$, but we cannot determine exactly where this would occur since such models are not available in the literature.

### 3.2 Protostellar Evolution with a mass dependent accretion rate

The results described so far have been obtained under the rather limiting assumptions that the accretion rate is constant throughout the evolution and of large magnitude ($`10^3M_{\mathrm{}}\mathrm{yr}^1`$), at least two orders of magnitude more than what is normally considered in the collpase phase of present-day molecular cores. That $`\dot{M}_{\mathrm{acc}}`$ is so large depends on the thermal and chemical properties of the primordial gas, which is known to be a very poor radiator owing to the lack of molecules and dust grains. Since the gas temperature during collapse is mainly determined by the excitation of H2 molecules, it is natural to expect that collapse models find typical values well in excess of several hundred degree, thus justifying the large values of $`\dot{M}_{\mathrm{acc}}`$ through eq. (1). On the other hand, the use of a single value of the accretion rate is less motivated since, as we have noted in §1, numerical simulations have shown a rather strong dependence on time of the amount of material that can actually collapse, even in the case of strict isothermality.

Thus, it is natural to ask how robust is the sequence of events that we have described in the previous sections in view of a more general behavior of the gas dynamics, and hence of the time history of the accretion rate. Since the key factor that determines the ultimate fate of massive protostars is whether $`\dot{M}_{\mathrm{acc}}`$ is above or below the critical value $`\dot{M}_{\mathrm{crit}}4\times 10^3M_{\mathrm{}}`$ yr-1, we should consider how and when this situation is met under more realistic conditions.

A useful guide is provided by the recent 3D hydrodynamical simulations of the formation of the first stars presented by Abel, Bryan & Norman (2002; hereafter ABN). Starting from cosmological initial conditions, these authors have followed the evolution of a collapsing primordial clump up to the stage when optical depth effects become very large, i.e., when the central number density reaches $`10^{12}\mathrm{cm}^3`$. Due to numerical limitations and their optically thin approximation for H2 line radiation, ABN could not follow correctly the subsequent history of protostar formation and accretion phases. However, they provide the run of the accretion time as a function of enclosed gas mass which can be converted into an effective accretion rate that depends on the growing protostellar mass. Using the results shown in Fig. 5 of ABN, we obtain the following fit

$$\dot{M}_{\mathrm{ABN}}=\{\begin{array}{cc}7.76\times 10^3(M_{}/M_{\mathrm{}})^{0.24}M_{\mathrm{}}\mathrm{yr}^1\hfill & 0.6M_{\mathrm{}}<M_{}<13M_{\mathrm{}}\\ 1.51\times 10^3(M_{}/M_{\mathrm{}})^{0.41}M_{\mathrm{}}\mathrm{yr}^1\hfill & 13M_{\mathrm{}}<M_{}<60M_{\mathrm{}}\\ 11.2(M_{}/M_{\mathrm{}})^{1.76}M_{\mathrm{}}\mathrm{yr}^1\hfill & 60M_{\mathrm{}}<M_{}.\end{array}$$
(10)
In order to avoid computational difficulties, we have used $`\dot{M}_{\mathrm{acc}}=8.8\times 10^3M_{\mathrm{}}`$ yr-1 for lowest mass models ($`M_{}<0.6M_{\mathrm{}}`$). Thus, the accretion rate starts at very high values, it remains well above $``$10$`{}_{}{}^{3}M_{\mathrm{}}^{}`$ yr-1 up to about 60 $`M_{\mathrm{}}`$, and then falls rapidly for very massive stars. Most importantly, $`\dot{M}_{\mathrm{ABN}}`$ drops below the critical value for $`M_{}>95M_{\mathrm{}}`$.

The evolution of the protostar with the accretion rate given by equation (10) is displayed in Figure Formation of the First Stars by Accretion. For comparison, we also show the results at fixed $`\dot{M}_{\mathrm{acc}}`$ discussed before. In Figure Formation of the First Stars by Accretion, the evolution of the photospheric radius is compared to that of the core. During the initial adiabatic phase, the protostar expands more slowly than in the constant $`\dot{M}_{\mathrm{acc}}`$ case in response to the decrease of $`\dot{M}_{\mathrm{ABN}}`$ for $`M_{}<13M_{\mathrm{}}`$. The thermal relaxation and the appearance of the luminosity wave with the expansion of the radius occur at approximately the same mass as for $`\dot{M}_{\mathrm{acc}}=\dot{M}_{\mathrm{fid}}`$, since the instantaneous value of $`\dot{M}_{\mathrm{ABN}}`$ is nearly the same as the fiducial one ($`0.97`$ $`\dot{M}_{\mathrm{fid}}`$ at 13 $`M_{\mathrm{}}`$).

During Phase II, the protostar contracts more gradually than in the fiducial case since $`\dot{M}_{\mathrm{ABN}}`$ reaches a minimum at $`13M_{\mathrm{}}`$ and increases thereafter. The core shrinks down to its minimum value of $`20R_{\mathrm{}}`$ at about 50 $`M_{\mathrm{}}`$ and then rebounds when radiation pressure effects begin to appear. This is the most delicate part of the evolution that depends critically on the behavior of $`\dot{M}_{\mathrm{ABN}}`$. As Fig. Formation of the First Stars by Accretion shows, the phase of violent expansion is avoided in this case and the protostar resumes further contraction. At about $`90M_{\mathrm{}}`$, hydrogen is ignited at the center and a convective core begins to develop. Then, the accretion rate falls below the critical value at $`95M_{\mathrm{}}`$. Shortly after, the protostar relaxes to a ZAMS star, without suffering from surface expansion or a blow out of the accretion envelope. Therefore, accretion is expected to continue unimpeded to very high masses.

With continuing accretion, the final mass of the forming star is set by how much material can be accreted during the stellar lifetime (e.g., Larson & Starrfield 1971). Because of the sharp decrease of the accretion rate, the accretion time $`t_{\mathrm{acc}}=M_{}/\dot{M}_{\mathrm{ABN}}3\times 10^4(M_{}/100M_{\mathrm{}})^{2.76}\mathrm{yr}`$ (for $`M_{}>60M_{\mathrm{}}`$) exceeds the stellar lifetime ($`3\times 10^6`$ yr) at $`M_{}500M_{\mathrm{}}`$. Therefore, even if accretion goes on throughout the entire stellar lifetime, the stellar mass would not exceed $`500M_{\mathrm{}}`$ by a large margin. A more detailed analysis has concluded that if the accretion continues, the mass attained during the entire stellar lifetime is about $`600M_{\mathrm{}}`$ (ABN). In reality, when the accretion rate drops below a certain value, the star is hot enough to produce a copious amount of ionizing photons that create an ordinary HII region, whose expansion can effectively reverse the infalling gas. Omukai & Inutsuka (2002) have shown that the mass limit set by this process is very similar to the limit imposed by the main-sequence lifetime ($`460M_{\mathrm{}}`$; see their eq. 27).

Returning to the evolution of the photospheric radius, in the high mass regime ($`M_{}>60M_{\mathrm{}}`$), the rapidly decreasing accretion rate causes the abrupt turn around shown in Fig. Formation of the First Stars by Accretion, after reaching the peak value of $`1000R_{\mathrm{}}`$ at about $`100M_{\mathrm{}}`$. At a mass of $`500M_{\mathrm{}}`$, the radiative precursor disappears and the accretion envelope becomes optically thin: At this point the hot stellar surface becomes directly visible from outside.

## 4 The formation of second-generation stars:
Effects of metal enrichment

Second-generation stars are formed out of material that has been contaminated by metals released from Population III SNe. Heavy elements affect the chemistry and dynamics of the gas in several ways. For example, detailed models have shown that a metal abundance of $`10^2Z_{\mathrm{}}`$ can substantially change the cooling properties of gas, thereby lowering the minimum mass of cosmological objects (e.g., Fall & Rees 1985). A metal abundance as low as $`10^4Z_{\mathrm{}}`$ can influence the fragmentation mass scale of star-forming clouds (Omukai 2000; Bromm et al. 2001a; Schneider et al. 2002). Protostellar evolution is also altered through important changes in the opacity and in the initial abundance of the CN elements. Therefore, it is important to consider how even a minute enrichment of heavy metals modifies the picture of star formation that we have been describing so far.

Since the preceding discussion has demonstrated that the qualitative evolution of accreting protostars is the same for fixed and variable mass accretion rates, we have run several models at constant $`\dot{M}_{\mathrm{acc}}`$ for nonzero metallicities. Figure Formation of the First Stars by Accretion shows the resulting mass-radius relation obtained with $`\dot{M}_{\mathrm{acc}}=1/2\dot{M}_{\mathrm{fid}}=2.2\times 10^3M_{\mathrm{}}`$ yr-1 for $`Z=10^6=5\times 10^5Z_{\mathrm{}}`$ and $`10^4=5\times 10^3Z_{\mathrm{}}`$, along with the standard zero-metal case. The evolution is basically the same, except for the temperature at which the CN-cycle starts operating. At higher $`Z`$, significant burning begins at lower temperatures. Therefore, the contraction phase is less pronounced, yielding a larger core radius on the main sequence. This explains the shift of the minimum radius in Figure Formation of the First Stars by Accretion at higher $`Z`$. On the other hand, the protostellar luminosity at a given mass is the same in all models. The main effect of metals consists of an increase of the critical accretion rate that varies from $`\dot{M}_{\mathrm{crit}}=4\times 10^3`$ M yr-1 for $`Z=0`$ to $`9\times 10^3M_{\mathrm{}}`$ yr-1 for $`Z=10^4=5\times 10^3Z_{\mathrm{}}`$. The change of a factor $``$2 in $`\dot{M}_{\mathrm{crit}}`$ reflects the similar variation of the protostellar radius (through equation 8). Of course, a higher value of the critical accretion rate implies that protostars can keep accreting without entering the phase of violent expansion.

The effect of metals is much more dramatic if we consider other values of $`\dot{M}_{\mathrm{acc}}`$. Figure Formation of the First Stars by Accretion displays the mass-radius relation computed with $`\dot{M}_{\mathrm{acc}}=\dot{M}_{\mathrm{fid}}=4.4\times 10^3M_{\mathrm{}}`$ yr-1 for the three metal compositions. The accretion rate is now slightly higher than the critical value for the metal-free case, but smaller than that at $`Z=10^6=5\times 10^5Z_{\mathrm{}}`$ ($`6\times 10^3M_{\mathrm{}}`$ yr-1). As a result, the violent expansion that occurs at $`Z=0`$ is no longer found, even in the lowest metallicity case.

Thus, we conclude that the presence of metals, by increasing the critical accretion rate, favors the formation of very massive objects. This result is not intuitive since metals, locked up in grains, are usually believed to reduce the upper mass limit because of the enhanced radiation force onto the dust shell surrounding the protostar (e.g., Kahn 1974; Wolfire & Cassinelli 1987). To understand this behavior, let us remind that the mechanism that should prevent infall is the violent expansion of the surface layers caused by the radiation force on the ionized radiative precursor. In the case of high $`Z`$, radiation pressure is imparted onto the dust shell which is located at a much greater distance from the precursor.

We can derive a simple estimate of the critical metallicity that separates the two regimes. The radiation force on the dust shell becomes bigger than that on the ionized precursor if the opacity in the shell $`\kappa _{\mathrm{dust}}=\kappa _{\mathrm{dust}}^{(0)}(Z/Z_{\mathrm{}})`$, where $`\kappa _{\mathrm{dust}}^{(0)}`$ 30 cm2 g-1 is the present value of dust opacity (e.g., Beech & Mitalas 1994), exceeds that due to electron scattering in the ionized precursor, $`\kappa _{\mathrm{prec}}=`$ 0.4 cm2 g-1. Here, we have assumed that the dust properties (e.g., composition, size distribution, etc.) do not change with metallicity. In this case, the critical value of metallicity is given by

$$Z>Z_{\mathrm{dr}}=0.01Z_{\mathrm{}}\left(\frac{\kappa _{\mathrm{dust}}^{(0)}}{30\mathrm{cm}^2\mathrm{g}^1}\right)^1.$$
(11)
Note that Wolfire & Cassinelli (1987) argued that the dust must be depleted by a factor of $`4`$ from the average interstellar value for the formation of massive stars. If so, the opacity will be lower by the same factor, and this will affect the critical metallicity $`Z_{\mathrm{dr}}`$ as well. For $`Z>Z_{\mathrm{dr}}`$, the infall is reversed by blowing out the outer envelope before the luminosity reaches the Eddington limit due to electron scattering. In this regime, metals act against prolonged accretion. For $`Z<Z_{\mathrm{dr}}`$, the radiation force onto the outer envelope is always below gravity, since the protostellar luminosity never exceeds the Eddington limit. When $`L_{}L_{\mathrm{Edd}}`$, the radiation force on the radiative precursor and on the surface layers causes the violent expansion of the protostar. In this regime, metals work in favor of continuing accretion by increasing the critical accretion rate.

## 5 Discussion

Much of the recent attention on primordial stars has focussed on the properties of massive objects. Several authors have argued that the first generation of cosmic structures formed mainly massive and very massive stars, and indirect evidence for a top-heavy IMF includes chemical abundances and supernova rates in our Galaxy, as well as star formation rates in high redshift galaxies (e.g., Larson 1998). Also, the apparent lack or dearth of stars in the galactic halo with metallicity less than $`Z/Z_{\mathrm{}}10^4`$, coupled with the theoretical prediction of a critical metallicity for efficient fragmentation at essentially the same value of $`Z`$ (Bromm et al. 2001a) have led to the suggestion that the IMF at early times differed substantially from its present form. Finally, the initial chemical enrichment of our Galaxy and the abundance pattern in the intracluster medium have been considered as compelling evidence for a generation of stars with mass $``$100 M (Schneider et al. 2002). In particular, these authors have postulated that in order to explain the apparent average minimum metallicity of $`Z/Z_{\mathrm{}}10^4`$, the first stars should have had masses in the rather narrow interval between $`140M_{\mathrm{}}`$ and $`260M_{\mathrm{}}`$, ending their lives as pair-unstable supernovae and polluting the pristine gas with the right amount of heavy elements. Higher and lower mass stars would in fact produce massive black holes, without releasing metals (e.g., Heger & Woosley 2002).

While there are good theoretical reasons to believe that massive stars were common in the early universe, the question of their actual formation has not been studied in detail. As we know from studies of the local universe, forming massive stars is a rare event that requires special conditions and it is not clear that inefficient fragmentation, lack of metals and fast mass accretion are per se sufficient ingredients to justify the expectations. The discovery of an extremely metal-poor halo star of $``$solar mass with an iron abundance \[Fe/H\]=$``$5.3, twenty times lower than the previous limit of $`Z/Z_{\mathrm{}}10^4`$ (Christlieb et al. 2002), clearly indicates that low-mass star formation in such extreme environments is more pervasive than otherwise believed.

The results of our calculations have shown that the final outcome of the accretion process is a sensitive function of the mass accretion rate and of its time variation. Notwithstanding the role of other important physical processes that we have neglected (presence of circumstellar disks, turbulence and magnetic fields), even in the idealized case in which star formation is mainly determined by $`\dot{M}_{\mathrm{acc}}`$ leads to uncertain predictions on the characteristic mass scale of the first stars. The existence of a critical mass accretion rate $`\dot{M}_{\mathrm{crit}}4\times 10^3M_{\mathrm{}}`$ yr-1 separates solutions in which very massive objects with mass $`100M_{\mathrm{}}`$ can form, provided there is sufficient material in the parent clumps, from others where the upper mass limit decreases as $`\dot{M}_{\mathrm{acc}}`$ increases. In the case of the time dependent accretion rate from the numerical simulations of Abel et al. (2002), protostellar growth can continue unimpeded by radiation forces onto the radiative precursor. However, slight deviations of the actual accretion rate from the empirical form of eq. (10) may result in a large difference in the final mass of the formed stars. Obviously, more refined models of the collapse phase should be carried out to verify the robustness or uniqueness of the behavior of the accretion rate in the high-mass regime ($`M_{}60M_{\mathrm{}}`$) which is critical for prolonged accretion. If the decrease of $`\dot{M}_{\mathrm{acc}}`$ is slower than that given by eq. (10), accretion can be stopped at around a few tens to 100 $`M_{\mathrm{}}`$, thus not allowing the formation of the objects which are considered responsible for the initial enrichment of the primordial gas.

Another important result of our study is the appreciation of the sensitivity of the evolution on the metal abundance. In the $`Z=0`$ case, the value of the actual accretion rate from eq. (1) is similar to the critical one. For the metal enriched case, $`\dot{M}_{\mathrm{acc}}`$ is expected to be smaller, as a result of the reduced gas temperatures of the prestellar clumps. On the other hand, $`\dot{M}_{\mathrm{crit}}`$ increases at higher $`Z`$, as shown in §4. Hence, in a slightly metal-enriched gas, $`\dot{M}_{\mathrm{acc}}`$ would be lower than $`\dot{M}_{\mathrm{crit}}`$, leading to continued accretion and favoring the formation of very massive stars.

Considering the mass scale for cloud fragmentation, recent studies have found a transition from massive ($`10^3M_{\mathrm{}}`$) to low-mass ($`1M_{\mathrm{}}`$) clumps at a metallicity of $`Z_{\mathrm{frag}}10^4Z_{\mathrm{}}`$ (Bromm et al. 2001a; Schneider et al. 2002). However, as discussed in the previous section, there appears to be another critical metallicity. If $`ZZ_{\mathrm{dr}}0.01Z_{\mathrm{}}`$, radiation pressure onto the dust shell surrounding a protostar becomes important and limits the amount of matter that can be accreted. Therefore, one can argue that there are three regimes that determine the maximum value of the stellar mass: for $`Z<Z_{\mathrm{frag}}`$, the limit is determined by the the stellar lifetime and/or by the expansion of the ionized region; for $`Z_{\mathrm{frag}}<Z<Z_{\mathrm{dr}}`$ by the fragmentation scale, and by the reversal of the gas inflow due to radiation pressure onto the dust shell for $`Z_{\mathrm{dr}}<Z`$. In the first case, massive stars form predominantly, while in the other two cases their production is possible, although not typical.

In this work, we have found that if $`\dot{M}_{\mathrm{acc}}>\dot{M}_{\mathrm{crit}}`$, radiation pressure on the ionized outer layers of the protostar causes a violent expansion of the radius. Since we have treated protostellar evolution as a series of hydrostatic configurations, the final outcome of this dynamic phase could not be fully investigated. We have suggested that the rapid swelling leads to the ejection of the external layers, to the onset of powerful mass loss through a wind, and eventually to the complete halt of accretion. Note that stellar pulsations due to the $`ϵ`$-mechanism triggered by the onset of the CN-cycle can lead to moderate mass loss in metal-rich stars, but in the case of metal-free stars the loss rate is strongly reduced (Baraffe, Heger, & Woosley 2001).

There are, however, other possibilities that need to be considered. For example, the stellar wind may not be strong enough to blow away all the envelope material. With no accretion, the central star would quickly contract to the main sequence, decreasing substantially its mass loss activity. Then, the remaining circumstellar matter would be able to collapse and recurrent accretion might be possible.

Deviations from spherical symmetry is also crucial. The stellar wind might blow in a bipolar jet-like configuration, without affecting accretion from the disk. Even if the accretion rate through the disk is higher than the critical rate, the protostar would adjust itself in a state marginally below $`\dot{M}_{\mathrm{crit}}`$ by shedding the remainder material as a jet-like wind. Such energetic jet could disrupt the entire parent clump and quench the supply of matter to the disk (e.g., Nakano, Hasegawa, & Norman 1995).

When the accretion proceeds in a disk-like fashion, the star loses energy faster than in the case of spherical accretion since the stellar radiation can escape in the optically thin polar directions. The reduced entropy level in the star results in a smaller protostellar radius (Palla & Stahler 1993; Hartmann, Cassen, & Kenyon 1997). Thus, at the same protostellar mass, the accretion luminosity is higher in the disk accretion than in the spherical case. This suggests that accretion in the polar direction is halted earlier on. On the other hand, if the parent clump is not disrupted by sweeping-up of the matter in the polar direction, accretion can proceed through the disk since only the matter at large angles from the disk is subject to the radiation force (Nakano 1989). To decide whether the disk-like configuration is more suitable or not for continuing accretion, 2D or 3D hydrodynamical studies must be carried out, similar to those now available for the present-day case (e.g., Yorke & Sonnhalter 2002). Recently Tan & McKee (2002) studied rotational effects on halting the protostellar accretion onto metal-free stars using a simplified treatment. They point out the difficulty of forming very massive objects (M¿¿100 Msun) because of the flow reversal in the polar directions. However, we believe that this conclusion is uncertain since accretion can still proceed through the circumstellar disk. Thus, the problem of making very massive stars is still open.

## 6 Conclusions

We have studied the evolution of extremely low metallicity ($`Z/Z_{\mathrm{}}=05\times 10^3`$) protostars accreting at very fast rates $`\dot{M}_{\mathrm{acc}}10^310^2M_{\mathrm{}}`$ yr-1, as is expected to occur in the first star forming sites of cosmological structures. The main results can be summarized as follows:

* The earliest stages of the protostellar growth are qualitatively the same independent of the details of the evolution of the mass accretion rate, or of the presence of trace abundances of heavy elements. These phases correspond to the adiabatic accretion, the propagation of the luminosity wave, and the Kelvin-Helmholtz (KH) contraction. However, the mass range covered in each phase does depend on $`\dot{M}_{\mathrm{acc}}`$, and in the parameter space explored here extends up to $`60M_{\mathrm{}}`$.
* For more massive protostars, the evolution follows two fundamentally different branches, depending on the instantaneous value of the mass accretion rate. If $`\dot{M}_{\mathrm{acc}}<\dot{M}_{\mathrm{crit}}4\times 10^3M_{\mathrm{}}`$ yr-1, protostars relax to ordinary H-burning objects while still accreting. In the opposite case, $`\dot{M}_{\mathrm{acc}}>\dot{M}_{\mathrm{crit}}`$, the protostellar luminosity reaches the Eddington limit before the onset of nuclear burning at the center. This phase is soon followed by a rapid expansion of the radius, possibly leading to reversal of the accretion flow.
* The two regimes yield upper mass limits that differ considerably: for $`\dot{M}_{\mathrm{acc}}\dot{M}_{\mathrm{crit}}`$, the maximum mass is $`M_{\mathrm{max}}300M_{\mathrm{}}`$ and decreases with $`\dot{M}_{\mathrm{acc}}`$; in the opposite case, $`M_{\mathrm{max}}300M_{\mathrm{}}`$, allowing the formation of very massive objects.
* Under a more realistic case of time dependent accretion, such as that proposed by Abel et al. (2002) where an initially high $`\dot{M}_{\mathrm{acc}}`$ begins to decrease only at high protostellar masses ($`60M_{\mathrm{}}`$), the evolution resembles the case $`\dot{M}_{\mathrm{acc}}<\dot{M}_{\mathrm{crit}}`$. Therefore, the star is expected to grow in mass until the accretion rate drops substantially and the protostar joins the main sequence. At this time, the expansion of an ionized region would likely stop further accretion ($`460M_{\mathrm{}}`$; Omukai & Inutsuka 2002). This limit is similar to that imposed by the amount of mass that can be accreted by the star during the entire main-sequence lifetime ($`600M_{\mathrm{}}`$; Abel et al. 2002).
* In all cases considered, deuterium burning and the $`pp`$-chain do not produce enough energy to counteract the effect of Kelvin-Helmoltz contraction. It is only with the beginning of the CN-cycle that nuclear burning has significant consequences for the protostellar evolution.
* A slight amount of metals ($`Z/Z_{\mathrm{}}10^410^2`$) favors continuing accretion by increasing the value of the critical accretion rate below which accretion can continue unimpeded. Thus, a non-zero metallicity favors the formation of very massive protostars. Above a threshold metallicity of $`Z0.01Z_{\mathrm{}}`$, the situation is reversed because of the feedback of radiation pressure on dust grains, as in the standard solar metallicity case.
* Almost throughout the protostellar phase, the accreting envelope is optically thick. Although the stellar surface is as hot as $`10^5`$ K, the photospheric temperature remains limited to $``$6000 K. Only after the protostellar mass exceeds $`100M_{\mathrm{}}`$ and the density in the envelope decreases substantially ($`5\times 10^{14}\mathrm{g}\mathrm{cm}^3`$ at the photosphere), the photospheric temperature begins to increase. In the case of the Abel et al. evolution, the photosphere coincides with the stellar surface for protostellar masses above $`500M_{\mathrm{}}`$. Beyond this value, the hot stellar surface becomes optically visible.
* Because of the high effective temperatures ($`10^5`$K), population III stars have been considered an efficient source of H and He ionizing photons, particularly in the case of star clusters consisting only of primordial ZAMS stars (Tumlinson & Shull 2000; Bromm et al. 2001b; Schaerer 2002). However, during the main accretion phase, the ionizing photons are degraded to optical photons in the precursor and cannot create a standard HII region. If the accretion phase is prolonged, the total emissivity of ionizing photons might be reduced by a non-negligible fraction with important consequences on the history of the re-ionization of the universe.

The work of K.O. is supported in part by Research Fellowships of the Japan Society for the Promotion of Science for Young Scientists, grant 6819. The research of F.P.has been partly supported by grant COFIN 2002-MIUR to INAF–Osservatorio Astrofisico di Arcetri.
Figure Caption


text:
# Nuclear absorption of Charmoniums in pA and AA collisions

## I Introduction

Recently, in Quark Matter 2002, NA50 collaboration presented their analysis of the nuclear absorption of $`J/\psi `$ in high statistics 450 GeV pA collisions . They estimated the $`J/\psi `$ nucleon absorption cross section ($`\sigma _{abs}^{J/\psi N}`$) in the framework of Glauber model. High statistics 450 GeV pA data yield $`\sigma _{abs}^{J/\psi N}=4.4\pm 1.0`$ mb . They also estimate a common $`\sigma _{abs}^{J/\psi N}`$ from latest pA and NA38 200 GeV/c S+U data , $`\sigma _{abs}^{J/\psi N}`$ =4.4 $`\pm `$ 0.5 mb. The extracted absorption cross section is much smaller than the earlier value of 6.4 $`\pm `$ 0.8 mb extracted from fit to earlier NA50 data  or 7.1 $`\pm `$ 3.0 mb obtained from a fit to NA38 S+U data . Within error, the S+U cross sections are compatible with pA cross sections. In Quark matter 2002, NA50 collaboration also presented their preliminary analysis of 2000 Pb+Pb run . $`E_T`$ dependence of the $`J/\psi `$ over Drell-Yan ratio, in the rapidity range of 2.9-4.5 was presented. Compared to 1998 run , 2000 data are flatter, suppression being more at low $`E_T`$ and less at high $`E_T`$. The difference may be attributed to the difference in method of analysis. In 1998, NA50 collaboration presented combined result of standard as well minimum bias analysis. In the 2000 run, only the standard analysis was performed. Centrality dependence of $`J/\psi `$ absorption is still anomalous in the sense that Glauber model of nuclear absorption fails to explain the data. The 1998 data of NA50 collaboration  were analyzed in a variety of models, with and without assumption of QGP . We have shown that a QCD based model, where $`J/\psi `$’s are absorbed in nuclear medium could explain the data . Parameters of the model were obtain from a fit to then existing pp/pA/AA total $`J/\psi `$ cross section data .

Preliminary analysis of high statistics pA data  with the implication of less absorption of $`J/\psi `$ in nuclear medium, along with the changed suppression pattern observed in the preliminary analysis of NA50 2000 Pb+Pb run , necessitated the re-examination of the QCD based nuclear absorption model . In the present brief report, we have analyzed the latest NA50 data on pA/SU total $`J/\psi `$ cross section  to obtain the parameters of the QCD based nuclear absorption model. The model was then used to explain the centrality dependence of $`J/\psi `$ over Drell-Yan ratio in 200 GeV/c S+U collisions and in 158 GeV/c Pb+Pb collisions (preliminary data). As will be shown in the following, once the model parameters are fixed from $`J/\psi `$ total cross sectional data, the QCD based nuclear absorption model give a consistent description of the centrality dependence of $`J/\psi `$ suppression in S+U collisions as well as in Pb+Pb collisions. The plan of the paper is as follows: in section 2, we will briefly describe the model and obtained the model parameters from the new high statistics pA/AA NA50 data. In section 3, the NA50 data on $`E_T`$ dependence of $`J/\psi `$ over Drell-Yan ratio in S+U and in Pb+Pb collisions are analyzed. Summary and conclusions are drawn in section 4.

## II $`J/\psi `$ production and suppression

In the QCD based nuclear absorption model , $`J/\psi `$ production is assumed to be a two step process, (a) formation of a $`c\overline{c}`$ pair, which is accurately calculable in QCD and (b) formation of a $`J/\psi `$ meson from the $`c\overline{c}`$ pair, which is non-perturbative but can be conveniently parameterized. The $`J/\psi `$ cross section in $`AB`$ collisions, at center of mass energy $`\sqrt{s}`$ is written as,

$`\sigma ^{J/\psi }(s)`$ $`=K{\displaystyle \underset{a,b}{}}{\displaystyle 𝑑q^2\left(\frac{\widehat{\sigma }_{abcc}}{Q^2}\right)𝑑x_F\varphi _{a/A}(x_a,Q^2)}`$ (2)
$`\varphi _{b/B}(x_b,Q^2){\displaystyle \frac{x_ax_b}{x_a+x_b}}\times F_{c\overline{c}J/\psi }(q^2),`$
where $`_{a,b}`$ runs over all parton flavors, and $`Q^2=q^2+4m_c^2`$. The $`K`$ factor takes into account the higher order corrections. The incoming parton momentum fractions are fixed by kinematics and are $`x_a=(\sqrt{x_F^2+4Q^2/s}+x_F)/2`$ and $`x_b=(\sqrt{x_F^2+4Q^2/s}x_F)/2`$. $`\widehat{\sigma }_{abc\overline{c}}`$ are the sub process cross section and are given in . $`F_{c\overline{c}J/\psi }(q^2)`$ is the transition probability that a $`c\overline{c}`$ pair with relative momentum square $`q^2`$ evolve into a physical $`J/\psi `$ meson. It is parameterized as,

$`F_{c\overline{c}J/\psi }(q^2)=`$ $`N_{J/\psi }\theta (q^2)\theta (4m_{}^{}{}_{}{}^{2}4m_c^2q^2)`$ (4)
$`(1{\displaystyle \frac{q^2}{4m_{}^{}{}_{}{}^{2}4m_c^2}})^{\alpha _F}.`$
In a nucleon-nucleus/nucleus-nucleus collision, the produced $`c\overline{c}`$ pairs interact with nuclear medium before they exit. It is argued  that the interaction of a $`c\overline{c}`$ pair with nuclear environment increases the square of the relative momentum between the $`c\overline{c}`$ pair. As a result, some of the $`c\overline{c}`$ pairs can gain enough relative square momentum to cross the threshold to become an open charm meson. Consequently, the cross section for $`J/\psi `$ production is reduced in comparison with nucleon-nucleon cross section. If the $`J/\psi `$ meson travel a distance $`L`$, $`q^2`$ in the transition probability is replaced to $`q^2q^2+\epsilon ^2L`$, $`\epsilon ^2`$ being the relative square momentum gain per unit length. Parameters of the model ($`\alpha _F`$,$`KN_{J/\psi }`$ and $`\epsilon ^2`$) can be fixed from experimental data on total $`J/\psi `$ cross section in pA/AA collisions. Earlier NA50 data  are well described with $`KN_{J/\psi }`$=0.458, $`\epsilon ^2=0.225GeV^2/fm`$ and $`\alpha _F=1.0`$ . As mentioned in the beginning, NA50 collaboration presented high statistics pA data on $`J/\psi `$ cross section. They have measured $`B_{\mu \mu }\sigma (J/\psi )/\sigma (DY)`$. In Fig.1, experimental data are shown as a function of nuclear length. The Drell-Yan cross sections donot have any A or alternately any L-dependence. The observed $`L`$ dependence is then due to $`J/\psi `$’s only. We fit the data with two parameters, $`N_{norm}=KN_{J/psi}/\sigma _{NN}^{DY}(nb)`$ and square momentum gain factor $`\epsilon ^2`$ ($`\alpha _F`$ being kept fixed at 1). In Fig.1, the fit obtained to the data is shown. The two sets of data at 200 GeV/c and 450 GeV/c could be fitted with a common square momentum gain factor, $`\epsilon ^2`$= 0.1875 $`GeV^2/fm`$, a value 20% lower than the value obtained earlier . Lowering of $`\epsilon ^2`$ indicate less absorption of $`J/\psi `$’s in nuclear medium, in agreement with the Glauber model calculations. While the square momentum gain factor do not show energy dependence, the evident energy dependence of the cross section ratios shows up in the other parameter of the model $`N_{norm}`$. We obtain $`N_{norm}`$ = 10.18 at 200 GeV/c and $`N_{norm}`$ = 4.43 at 450 GeV/c. The energy dependence of $`J/\psi `$ cross section being taken care of in the model (Eq.2), the energy dependence of $`N_{norm}`$ is due to the Drell-Yan cross sections only. In the mass range, $`2.9>M>4.5`$ GeV, the Craigie parameterization , of the Drell-Yan cross section, $`\sigma (DY)e^{14.9M/\sqrt{s}}`$, gives for the ratio $`\sigma (DY)_{450GeV}/\sigma (DY)_{200GeV}`$ = 2.1 -3.1, consistent with the presently obtained ratio of 2.29.

## III $`E_T`$ dependence of $`J/\psi `$/DY ratio

In the present section, we analyzed the $`E_T`$ dependence of $`J/\psi `$ over Drell-Yan ratio in 200 GeV/c S+U collisions and in 158 GeV/c Pb+Pb collisions. Details of the calculation can be found in ref..

In Fig.2, NA38 experimental ratio of $`B_{\mu \mu }\sigma (J/\psi )/\sigma (DY)`$ for 200 GeV/c S+U collisions  is compared with the present model calculation (solid line). The data are well described in the model. For comparison purpose, we have also shown the fit obtained to the data in the Glauber model. The dashed line is the Glauber model calculations with $`\sigma _{abs}^{J/\psi N}`$= 7.1 mb. However, quality of fit gets poorer if we use lower value (4.4 mb) for $`\sigma _{abs}^{J/\psi N}`$ as extracted from the most recent data (the dotted line). It is evident that in the Glauber model of nuclear absorption, the transverse energy dependence of $`J/\psi `$’s prefer larger absorption of $`J/\psi `$ in S+U collisions than indicated in pA/AA total $`J/\psi `$ cross section data. The QCD based nuclear absorption model, on the other hand, give consistent description of the data.

Transverse energy dependence of the ratio, $`B_{\mu \mu }\sigma (J/\psi )/\sigma (DY)`$ in 158 GeV/c Pb+Pb collisions could not be explained in the Glauber model of nuclear absorption. The preliminary NA50 data obtained in 2000, are flatter compared to 1996/1998 data , still the usual Glauber model of nuclear absorption fails to describe the data. In Fig.3 the experimental data are shown along with the Glauber model fit to it obtained with $`\sigma _{abs}^{J/\psi }`$ =4.4 mb. The Glauber model calculation (the dotted line) agree only for the peripheral collisions. For more central collisions, the model predict less absorption than in experiment. Also the Glauber model shows saturation at more central collisions, in contrast to experiment. The 2nd drop in the ratio is mainly due to transverse energy fluctuations at a fixed impact parameter . In the present model, we take into account $`E_T`$ fluctuations at a fixed impact parameter $`b`$, by the replacement: $`L(b,s)L(b,s)E_T/<E_T>(b)`$. In Fig.3, the dashed is obtained in the QCD based nuclear absorption model, without incorporating the $`E_T`$ fluctuation effect. It explains the data up to 100 GeV, the knee of the $`E_T`$ distribution, beyond which only, $`E_T`$ fluctuations are important. If $`E_T`$ fluctuation is included, the data are explained throughout the $`E_T`$ range (the solid line). Consistent description of $`E_T`$ dependence of $`J/\psi `$ over Drell-Yan ratio in 200 GeV/c S+U and 158 GeV/c Pb+Pb collisions clearly demonstrate that the QCD based nuclear absorption model, together with transverse energy fluctuation could explain consistently the $`J/\psi `$ absorption data at SPS energy.

## IV summary and Conclusions

To summarize, we have analyzed the latest NA50 data  on $`J/\psi `$ suppression obtained in 2000 Pb+Pb run. The data were analyzed in the QCD based nuclear absorption model. The $`J/\psi `$ production is assumed to be two step process, (i) production of $`c\overline{c}`$ pair, purturbatively calculable and (ii) formation of $`J/\psi `$ from the $`c\overline{c}`$ pair. In a nuclear medium, due to random collisions with medium, relative square momentum of the $`c\overline{c}`$ increases and some of the pair may gain enough square momentum to cross the threshold for open charm meson. We obtain the parameter, square momentum gain factor per unit length, $`\epsilon ^2`$, from a fit to the latest NA50 data , and obtain $`\epsilon ^2`$=0.1875 $`GeV^2/fm`$. The value is lower than the value $`\epsilon ^2`$= 0.225 $`GeV^2/fm`$, obtained in our earlier work . There, we had fitted the earlier version of the NA50 data . Reduction in $`\epsilon ^2`$ indicate that the $`J/\psi `$’s are less absorbed in nuclear medium. This is in accordance with the $`J/\psi `$-nucleon absorption cross section, extracted from the data . $`\sigma _{abs}^{J/\psi N}`$ extracted from the earlier version of the NA50 data  was large, $`\sigma _{abs}^{J/\psi N}=6.4\pm 0.8`$ mb, compared to the value extracted from the latest data set , $`\sigma _{abs}^{J/\psi N}=4.4\pm 1.0`$ mb. With the new model parameter, the model could explain the centrality dependence of $`J/\psi `$ over Drell-Yan ratio in 200 GeV/c S+U collisions as well as the latest NA50 data on the centrality dependence of $`J/\psi `$ over Drell-Yan ratio in 158 GeV/c Pb+Pb collisions.


text:
# Classicality and anticlassicality measures of pure and mixed quantum states

## 1 Introduction

Since Glauber’s paper , quantum states for which the so-called $`P`$-distribution  is nonpositive or more singular than delta function are called “nonclassical”. For the past decades, many authors proposed different quantitative measures of “nonclassicality”. It seems impossible to reduce all variety of quantum states in infinite-dimensional Hilbert spaces to some unique parameter. Therefore, different existing approaches should be considered sooner as complementary rather than competitive.

Historically, the first approach was based on the analysis of deviations from the Poissonian photon statistics inherent to the Klauder–Glauber–Sudarshan  coherent states, which are known to be the only “classical” pure states . The examples are Mandel’s “$`Q`$-parameter”  and its various generalizations . The second direction is to evaluate the volume of that part of phase space where some quasiprobability distribution assumes negative values  (see also ). There are also other approaches .

Here we follow the direction opened by Hillery . It consists in evaluating some kinds of distances in the Hilbert space between the state concerned and a family of states which are assumed to be “classical”. There exist many different definitions of the distance, besides the “trace distance” used in . In particular, distances along geodesics on curved manifold supplied with Riemannian metrics for coherent, squeezed, displaced, and other states were studied in . Other examples are the Monge distance , “polarized” and “classical-like” distances . The most simple from the point of view of calculations is the Hilbert–Schmidt distance, used in ,

$`d_{HS}^2(\widehat{\rho },\widehat{\rho }_c)`$ $`=`$ $`\mathrm{Tr}(\widehat{\rho }\widehat{\rho }_c)^2`$ (1)
$`=`$ $`\mathrm{Tr}(\widehat{\rho }^2)+\mathrm{Tr}(\widehat{\rho }_c^2)2\mathrm{Tr}(\widehat{\rho }\widehat{\rho }_c).`$
Here $`\widehat{\rho }`$ is the statistical operator of the quantum state concerned, and $`\widehat{\rho }_c`$ is related to the reference “classical” state.

For pure quantum states, $`\widehat{\rho }=|\psi \psi |`$, the reference states are usually identified with coherent states $`\widehat{\rho }_c=|\alpha \alpha |`$, as soon as the latter are assumed to be “the most classical” ones . Then the calculation of distance is reduced to calculating the scalar product $`\alpha |\psi `$. This scalar product is the principal ingredient of many other distances, such as the Fubiny–Study distance  $`d_{FS}^2=1|\psi _1|\psi _2|^2`$ (we use here slightly different notation), or Wootters’ distance  $`d_W=\mathrm{cos}^1(|\psi _2|\psi _1|)`$ (see  for more references).

However, coherent states (whose $`P`$-distributions are delta-functions) represent only a small subset “at the border” of the set of all “classical” states. Therefore, it seems reasonable to enlarge the family of reference states $`\widehat{\rho }_c`$, incorporating  all displaced thermal states, in complete correspondence with original paper .

Unfortunately, even the Hilbert–Schmidt distance becomes complicated in such a case, because the presence of three terms given in the second line of Eq. (1) does not permit to find the minimum of this expression analytically (except for the most simple cases), when one deals with mixed states. A possibility to simplify calculations by adjusting the “purity” $`\mathrm{Tr}(\widehat{\rho }_c^2)`$ of the reference states to the purity $`\mathrm{Tr}(\widehat{\rho }^2)`$ of the state concerned was discussed in .

The aim of our article is to consider another possibility, which permits us to find analytical expressions for most known families of nonclassical states. Namely, in the next section we introduce a new “classicality measure”, which is proportional to (but not identical with) the last term in (1), $`\mathrm{Tr}(\widehat{\rho }\widehat{\rho }_c)`$. We analyze this measure for the Fock and generic Gaussian mixed states. In Sec. 3 we consider the “anticlassicality measure”, based on using “the most nonclassical” Fock states as the reference ones. A brief discussion of results and perspectives is given in the concluding Sec. 4

## 2 Classicality measure

As a matter of fact, practically all information about the closeness of the states $`\widehat{\rho }`$ and $`\widehat{\rho }_c`$ is contained in the term $`\mathrm{Tr}(\widehat{\rho }\widehat{\rho }_c)`$, whereas the presence of the term $`\mathrm{Tr}(\widehat{\rho }_c^2)`$ in (1) only complicates the search of minimum of $`d_{HS}^2`$, adding no significant information. Therefore, it seems reasonable to consider only the term $`\mathrm{Tr}(\widehat{\rho }\widehat{\rho }_c)`$, searching for its maximum with respect to the enlarged family of “classical” states, consisting of all displaced thermal states

$$\widehat{\rho }_c=\widehat{D}(\alpha )\widehat{\rho }_{th}\widehat{D}^1(\alpha ).$$
(2)
Here $`\widehat{D}(\alpha )=\mathrm{exp}\left(\alpha \widehat{a}^{}\alpha ^{}\widehat{a}\right)`$ is the well-known displacement operator  written in terms of the bosonic ahhihilation and creation operators, $`[\widehat{a},\widehat{a}^{}]=1`$, and $`\widehat{\rho }_{th}`$ is the thermal state,

$$\widehat{\rho }_{th}=(1\eta )\underset{k=0}{\overset{\mathrm{}}{}}\eta ^k|kk|,0\eta <1.$$
(3)
If we confine ourselves to the special case of $`\eta =0`$, when $`\rho _c=|\alpha \alpha |`$, and consider only pure quantum states $`\widehat{\rho }=|\psi \psi |`$, then the quantity $`\mathrm{max}_\alpha |\alpha |\psi |^2`$ seems a good measure of “classicality”, because it equals $`1`$ for coherent states and is less than $`1`$ for all other pure states. However, dealing with the generic case of mixed states, one meets certain difficulties: on the one hand, the measure of classicality for any thermal state must be the same as for coherent states, i.e., it must be equal to $`1`$ according to the normalization chosen, but on the other hand, $`\mathrm{Tr}(\widehat{\rho }_c^2)<1`$. Therefore, one has to find a suitable generalization of the scalar product between pure states to the case of mixed states, which would result in the unit value for identical states.

One possible solution, satisfying many additional requirements, was given by Uhlmann , who showed that the generalization of the quantity $`|\psi _1|\psi _2|`$ can be taken in the form $`\text{Tr}\sqrt{\widehat{\rho }_1^{1/2}\widehat{\rho }_2\widehat{\rho }_1^{1/2}}`$. However, the calculation of this trace is rather complicated problem, although certain progress has been achieved recently for Gaussian states (in the case of infinite-dimensional Hilbert space) .

We prefere to follow a more straightforward way: simply to replace the operators $`\widehat{\rho }`$ and $`\widehat{\rho }_c`$ in $`\mathrm{Tr}(\widehat{\rho }\widehat{\rho }_c)`$ by the renormalized operators $`\widehat{\rho }^{}`$ and $`\widehat{\rho }_c^{}`$, where

$$\widehat{\rho }^{}\widehat{\rho }/\sqrt{\text{Tr}\left(\widehat{\rho }^2\right)}.$$
(4)
Therefore, we define the “classicality measure” as

$$𝒞=\underset{\rho _c}{\mathrm{max}}\text{Tr}\left(\widehat{\rho }^{}\widehat{\rho }_c^{}\right).$$
(5)
For pure states, $`\widehat{\rho }=|\psi \psi |`$, we have

$$𝒞=\underset{\rho _c}{\mathrm{max}}\psi |\widehat{\rho }_c^{}\psi .$$
(6)
In modern literature, the “generalized scalar products”, such as $`\text{Tr}\sqrt{\widehat{\rho }_1^{1/2}\widehat{\rho }_2\widehat{\rho }_1^{1/2}}`$ or $`\mathrm{Tr}(\widehat{\rho }_1\widehat{\rho }_2)`$, are frequently called  “fidelities” (“Bures–Uhlmann” or “Hilbert–Schmidt”, respectively). In particular, a possibility of using the “modified fidelity” $`\mathrm{Tr}(\widehat{\rho }_1^{}\widehat{\rho }_2^{})`$ was mentioned in . Note that if only one of the states is pure, then the “modified fidelity” does not coincide with $`\psi _1|\widehat{\rho }_2|\psi _1`$, but differs from it by the factor $`[\mathrm{Tr}(\widehat{\rho }_2^2)]^{1/2}`$. Sometimes  this is considered as a flaw. But in our problem the factor discussed has real physical meaning, because it permits us to distinguish mixed states from their pure “partners” possessing the same probability distribution functions of quanta (see examples below). Another difference between the cited papers and our approach is that we consider “fidelity” not to a fixed state, but to the whole family of “classical” states.

We shall use the notation $`f(\eta ,\alpha )`$ for the function $`\mathrm{Tr}(\widehat{\rho }^{}\widehat{\rho }_c^{})`$, whose maximum gives the value of classicality. Sometimes, it can be difficult to find a maximum with respect to both variables. Then one can calculate partial maxima with respect to one variable, fixing zero value for another. We shall denote such “reduced classicalities” as $`\stackrel{~}{𝒞}_\alpha `$ or $`\stackrel{~}{𝒞}_\eta `$, where the subscript indicates the parameter, with respect to which the maximum was calculated.

### 2.1 Classicality of the Fock states

For the Fock state $`|n`$, the quantity $`\mathrm{Tr}(\widehat{\rho }\widehat{\rho }_c)`$ is reduced to the probability $`n|\widehat{\rho }_c|n`$ of discovering $`n`$ quanta in the displaced thermal state. This probability was calculated by many authors, e.g., in . Using the results of  and taking into account the expression for purity of the thermal state (it does not depend on the shift parameter $`\alpha `$),

$$\text{Tr}\left(\widehat{\rho }_c^2\right)=\frac{1\eta }{1+\eta },$$
(7)
we arrive at the function

$`f_n(\eta ,\alpha )`$ $`=`$ $`\sqrt{1\eta ^2}\eta ^n\mathrm{exp}\left(|\alpha |^2\sqrt{1\eta ^2}\right)`$ (8)
$`\times L_n\left(|\alpha |^2(1\eta )^2/\eta \right),`$
where $`L_n(z)`$ is the Laguerre polynomial. For $`\eta =0`$ we have the well-known Poissonian distribution of the coherent state,

$$f_n(0,\alpha )=\mathrm{exp}\left(|\alpha |^2\right)|\alpha |^{2n}/n!,$$

which has maximum at $`|\alpha |^2=n`$. Therefore, confining ourselves to pure coherent reference states, we obtain for the reduced classicality of the Fock state the formula (cf. ) $`\stackrel{~}{𝒞}_\alpha ^{(n)}=e^nn^n/n!`$, which results in $`\stackrel{~}{𝒞}_\alpha ^{(n)}(2\pi n)^{1/2}`$ for $`n1`$.

On the other hand, the function $`f_n(\eta ,0)`$ has maximum at $`\eta _n=\sqrt{n/(n+1)}`$, so that

$$\stackrel{~}{𝒞}_\eta ^{(n)}=\sqrt{n^n/(1+n)^{n+1}}.$$
(9)
One can verify that $`\stackrel{~}{𝒞}_\eta ^{(n)}>\stackrel{~}{𝒞}_\alpha ^{(n)}`$ for any $`n1`$. In particular, $`\stackrel{~}{𝒞}_\eta ^{(1)}=1/2`$ and $`\stackrel{~}{𝒞}_\alpha ^{(1)}=1/e`$, whereas for $`n1`$ we have $`\stackrel{~}{𝒞}_\eta ^{(n)}(en)^{1/2}`$. Moreover, the function $`f_n(\eta ,\alpha )`$ has negative derivative with respect to $`|\alpha |^2`$ at the point $`|\alpha |^2=0`$ for any $`\eta _n`$. Therefore, the maxima of $`f_n(\eta ,\alpha )`$ at the points $`(\eta _n,0)`$ are global, whereas the maxima at the points $`(0,n)`$ are only local: see Fig. 1. This means that the classicality of the Fock state is given by (9), and that the pure state $`|n`$ is more close not to the pure (displaced) coherent state with $`\overline{n}_{coh}|\alpha |^2=n`$, but to the mixed (undisplaced) thermal state with

$$\overline{n}_{th}\eta /(1\eta )=\sqrt{n}\left(\sqrt{n}+\sqrt{n+1}\right)$$

(i.e., $`\overline{n}_{th}2n`$ for $`n1`$) and

$$\text{Tr}\left(\widehat{\rho }_c^2\right)=\left(\sqrt{n}+\sqrt{n+1}\right)^2.$$

### 2.2 Classicality of Gaussian states

Generic Gaussian states are characterized  by two variances of the quadrature components, $`\sigma _q`$ and $`\sigma _p`$, their covariance, $`\sigma _{pq}`$, and two displacement parameters in the phase plane $`qp`$. It is clear, however, that for the fixed (co)variances, the quantity $`\mathrm{Tr}(\widehat{\rho }\widehat{\rho }_c)`$ is maximal for coinciding displacement parameters of $`\widehat{\rho }`$ and $`\widehat{\rho }_c`$, which means that it is sufficient to consider the states with zero displacements, i.e., to find the maximum of the function of single variable $`F(\eta )f(\eta ,0)`$. Since the statistical operator of the thermal state is diagonal in the Fock basis, the function $`F(\eta )`$ can be easily expressed in terms of the generating function

$$G(z)\underset{n=0}{\overset{\mathrm{}}{}}p_nz^n$$

for the photon distribution $`p_nn|\widehat{\rho }|n`$:

$$F(\eta )=\mu ^{1/2}\sqrt{1\eta ^2}G(\eta ),\mu \mathrm{Tr}(\widehat{\rho }^2).$$
(10)
The function $`G(z)`$ for generic Gaussian states was calculated by means of different methods, e.g., in . Using the form given in , we arrive at the function

$$F(\eta )=2\left[\frac{\mu (1\eta ^2)}{a2b\eta +c\eta ^2}\right]^{1/2},$$
(11)
where

$$a=1+\mu ^2+2T\mu ^2,b=1\mu ^2,c=1+\mu ^22T\mu ^2,$$

$$\mu ^2=\left[4\left(\sigma _p\sigma _q\sigma _{pq}^2\right)\right]^1,T=\sigma _p+\sigma _q.$$

One can easily verify that the maximum of function (11) is attained for $`\eta _{}=(1\mu )/(1+\mu )`$, which means that the thermal state which is most close to the given Gaussian state has the same purity: $`\mu _c=\mu `$. The classicality of the Gaussian states is given by a simple formula

$$𝒞^{(G)}=\sqrt{\frac{2}{1+\mu T}}.$$
(12)
Note that $`\mu T1`$ for Gaussian states, and the equality holds for thermal states. The parameter $`T`$ is related to the mean number of quanta $`\overline{n}`$ in the unshifted Gaussian state as $`T=1+2\overline{n}`$. Therefore, the classicality of pure unshifted Gaussian states (which are nothing but squeezed vacuum states) equals $`𝒞^{(sqv)}=(\overline{n}+1)^{1/2}`$. If $`\mu <1`$, then the minimal mean number of quanta in the Gaussian state equals $`\overline{n}_{min}=(1\mu )/(2\mu )`$. However, mixed Gaussian states remain unsqueezed (that means, no one quadrature component can have the variance less than the vacuum state value $`1/2`$ in dimensionless units), if $`\overline{n}<\overline{n}_c=(1\mu ^2)/(2\mu ^2)`$. Such states possess positive $`P`$-distributions , and it was proposed in  to extend the family of “classical” states, including all mixed unsqueezed Gaussian states. We suppose to study such an “extended classicality” in another paper.

There exist pure states possessing the same photon distribution function as thermal ones; usually they are called coherent phase states :

$$|\epsilon =\sqrt{1|\epsilon |^2}\underset{n=0}{\overset{\mathrm{}}{}}\epsilon ^n|n.$$
(13)
Although the state (13) is most close (at least with respect to the reduced $`\stackrel{~}{𝒞}_\eta `$-measure) to the thermal state with $`\eta =|\epsilon |^2`$, its $`\eta `$-classicality is less than $`1`$ (precisely because we use renormalized operators in the definition (5), the function $`F(\eta )`$ in Eq. (10) contains the factor $`\mu ^{1/2}`$, but this factor equals unity now, as soon as $`\widehat{\rho }_\epsilon ^{}=\widehat{\rho }_\epsilon `$):

$$\stackrel{~}{𝒞}_\eta ^{(\epsilon )}=\sqrt{\frac{1|\epsilon |^2}{1+|\epsilon |^2}}=\left(1+2\overline{n}\right)^{1/2},$$

in accordance with our feeling that pure states, in general, are less classical than their mixed partners. For highly excited states with the same mean value $`\overline{n}1`$ we obtain the asymptotical relations

$$𝒞^{(sqv)}\sqrt{2}\stackrel{~}{𝒞}_\eta ^{(\epsilon )}\sqrt{e}𝒞^{(Fock)}.$$

## 3 Anticlassicality measures

As soon as the Fock states are usually considered as “the most quantum” states, we can use them as reference states for the definition of the degree of “anticlassicality”:

$$𝒜=\underset{n}{\mathrm{max}}n|\widehat{\rho }n.$$
(14)
Perhaps, it could be worth including the displaced Fock states, $`|\stackrel{~}{n}\widehat{D}(\alpha )|n`$, into the set of reference states, but here we shall not do it. Actually, we can define two measures: one given by (14), where integer $`n`$ runs over all integers, including $`n=0`$, and another (denoted as $`𝒜_1`$), where the vacuum state $`|0`$ (which is distinguished from all other states) is excluded from the set of reference states. To see the difference, let us consider an example of coherent states $`|\alpha `$. One can easily find that the dependence of $`𝒜`$ on the mean photon number $`\overline{n}=|\alpha |^2`$ has different analytical forms in the intervals $`k\overline{n}k+1`$, $`k=0,1,\mathrm{}`$:

$$𝒜^{(\alpha )}(\overline{n})=\mathrm{exp}\left(\overline{n}\right)\overline{n}^k/k!,k\overline{n}k+1.$$
(15)
The function given by (15) is continuous, but it has jumps of derivatives at the points $`\overline{n}_k=k`$: the right derivatives at these points are equal to zero, while the left derivatives have finite negative values. However, $`𝒜^{(\alpha )}(\overline{n})`$ goes to $`1`$ as $`\overline{n}0`$, and this may seem strange.

The function $`𝒜_1^{(\alpha )}(\overline{n})`$ coincides with $`𝒜^{(\alpha )}(\overline{n})`$ for $`\overline{n}1`$. But for $`\overline{n}<1`$ its behaviour is quite different:

$$𝒜_1^{(\alpha )}(\overline{n})=\overline{n}\mathrm{exp}\left(\overline{n}\right),0\overline{n}1.$$

The dependence $`𝒜_1^{(\alpha )}(\overline{n})`$ is shown in Fig. 2. Asymptotically, for $`\overline{n}1`$, we have (using the Stirling formula) $`𝒜_1^{(\alpha )}(2\pi \overline{n})^{1/2}`$.

As another example we consider the squeezed vacuum state. It is known that the photon distribution function $`p_n`$ of this state is different from zero only for even values $`n=2m`$, and it can be expressed in terms of the mean number of photons as

$$p_{2m}^{(sqv)}=(1+\overline{n})^{1/2}\frac{(2m)!}{(2^mm!)^2}\left(\frac{\overline{n}}{1+\overline{n}}\right)^m.$$
(16)
The sequence $`p_{2m}`$ decreases monotonously. Therefore we obtain a monotonous dependence

$$𝒜^{(sqv)}(\overline{n})=(1+\overline{n})^{1/2},$$

whereas the function

$$𝒜_1^{(sqv)}(\overline{n})=\frac{\overline{n}}{2(1+\overline{n})^{3/2}}$$

has maximum $`𝒜_{1max}^{(sqv)}=1/(3\sqrt{3})`$ for $`\overline{n}=2`$ (and the most close Fock state is $`|2`$).

For the coherent phase state (13) we obtain

$$𝒜^{(\epsilon )}(\overline{n})=(1+\overline{n})^1,𝒜_1^{(\epsilon )}(\overline{n})=\overline{n}(1+\overline{n})^2,$$

with $`𝒜_{1max}^{(\epsilon )}=1/4`$ at $`\overline{n}=1`$. According to definition (14), the same results must hold for thermal states, since they have the same photon distribution function $`n|\widehat{\rho }|n`$. However, we can distinguish pure and mixed states having the same diagonal elements in the Fock basis, playing again with the purity parameter $`\mu =\text{Tr}(\widehat{\rho }^2)`$. But in contrast to the case of classicality parameter, now we should, sooner, multiply the quantity $`n|\widehat{\rho }|n`$ by $`\mu `$ (or some power of $`\mu `$, if one wishes), rather than divide it by $`\mu `$ (or $`\sqrt{\mu }`$), because the intuition tells us that mixed states are more anticlassical than pure ones. Adopting the definitions

$$\stackrel{~}{𝒜}=\mu \underset{n}{\mathrm{max}}n|\widehat{\rho }n,\stackrel{~}{𝒜}_1=\mu \underset{n1}{\mathrm{max}}n|\widehat{\rho }n,$$
(17)
we obtain for the thermal states the expressions

$$\stackrel{~}{𝒜}^{(th)}(\overline{n})=\left[(1+\overline{n})(1+2\overline{n})\right]^1,$$

$$\stackrel{~}{𝒜}_1^{(th)}(\overline{n})=\frac{\overline{n}}{(1+\overline{n})^2(1+2\overline{n})}.$$

These three cases are illustrated in Fig. 3. The asymptotical dependences at $`\overline{n}1`$ are as follows:

$$𝒜_1^{(sqv)}\frac{1}{2\sqrt{\overline{n}}},𝒜_1^{(\epsilon )}\frac{1}{2\overline{n}},𝒜_1^{(th)}\frac{1}{2\overline{n}^2}.$$

For generic (mixed and shifted) Gaussian states, the dependences $`𝒜_1(\overline{n})`$ and $`\stackrel{~}{𝒜}_1(\overline{n})`$ may be more complicated, especially in the regime of strong irregular oscillations of the photon distribution function $`p_n`$ . We reserve this case for another study.

## 4 Conclusion

We have introduced new measures of “classicality” and “anticlassicality” of quantum states, which characterize the “fidelities” of the chosen state to the whole families of states considered as “classical” or “maximally quantum”. Due to simplicity, these measures can be calculated analytically for many important sets of quantum states, thus providing additional information on their properties. In particular, we have shown that the most close to arbitrary Gaussian states are thermal states with the same degree of quantum mixing (purity). At the same time, there are (mixed) thermal states which are more close to the Fock states than (pure) coherent states. As soon as the new measures distinguish coherent quantum superpositions (“cat states”) from quantum mixtures, they could be used in studies of decoherence processes. Also, it would be interesting to establish relations between the new and existing measures of classicality, such as Mandel’s parameter, “nonclassical depth” , and so on, or to try to “order” different families of quantum states in accordance with their degrees of (anti)classicality. We hope to report on the results of studies in these directions somewhere.

## Acknowledgments

VVD and MBR acknowledge support of the Brazilian agencies CNPq and FAPESP, respectively.