text:
# The Origin of Jovian Planets in Protostellar Disks: The Role of Dead Zones

## 1 Introduction

The discovery of exosolar Jovian-mass planets has stimulated intense efforts to understand how planets form (e.g. Lin et al., 1999; Wuchterl et al., 2000). The current theoretical models are based on two general types of physical processes – disk instability or core accretion. In the disk instability picture (e.g. Cameron, 1978; Boss, 1997; Mayer et al., 2002), gravitational instability in a sufficiently massive disk fragments it into the clumps out of which planets may form. This process may require as little as $`10^3`$ years to form Jovian planets, but may not be able to explain sub-Jupiter mass planets (Boss, 2000). This is in stark contrast to the core accretion model (e.g. Mizuno, 1980; Pollack et al., 1996) wherein the coagulation of dust grains ultimately leads to the formation of planetesimals that collide with one another to form planetary cores. When a core reaches a critical mass in the range $`1015M_{}`$, gas around the protoplanet quickly accretes onto it (“runaway accretion”) forming a gas giant. This process takes at least $`10^610^7`$ years – the average disk lifetime (e.g., Feigelson & Montmerle, 1999), and may not be able to explain massive Jovian planets.

Our paper focuses on an important aspect of Jovian planet formation that is common to both pictures. Regardless of whether proto-Jovian planets form through core accretion or gravitational instability, accretion onto a planet will continue until it is massive enough for the tidal torque that it exerts upon the surrounding disk to overcome the viscous disk torque which continues to fill-in any gap (e.g. Lin & Papaloizou, 1985). Planets that are sufficiently massive to open gaps in the disks have final masses that depend upon the disk viscosity (Nelson et al., 2000). In absence of viscosity, the tidal torque induced by the protoplanet opens a gap sooner thereby terminating the accretion earlier and leaving the planet with a smaller mass.

Thus the final mass of Jovian planets comes down to the question of the origin of the viscosity of protostellar disks at the time that significant protoplanetary cores are present within them. If disks are turbulent, then turbulent viscosity is likely to be the major factor in determining these masses. The Magneto-Rotational Instability (MRI) is the most promising source of turbulent disk viscosity (Balbus & Hawley, 1991). Even in the presence of an initially very weak magnetic field, the MRI will generate a significantly magnetized disk. This MRI viscosity requires good coupling between the gas and the magnetic field. Poor coupling – which occurs when the gas is insufficiently ionized – damps out the MRI and leads to a so-called “dead zone” where the viscosity is effectively zero (Gammie, 1996). The high column density for protoplanetary disks ensures that dead zones are expected to occur during planet formation (Sano et al., 2000).

Disk viscosity will not entirely vanish, even in the dead zone, however. The gravitational interaction between the fairly massive protoplanetary cores – the forerunners of the Jovian planets – and the disk generates density waves within the disk which can transport angular momentum and provide an effective “viscosity” in the disk. This disk-planet interaction leads to both planetary migration (Goldreich & Tremaine, 1980) as well as evolution of the disk through the nonlinear damping and shocking of the waves (Larson, 1989; Goodman & Rafikov, 2001, – GR01). The deposition of the angular momentum back into the disk as a consequence of this shock-induced wave damping is equivalent to the action of a viscous process, (Spruit, 1987; Larson, 1989) with an equivalent dimensionless $`\alpha `$ parameter of the order $`\alpha _{damp}10^410^3`$ (e.g. GR01). This viscosity parameter is one to two orders of magnitude smaller than is expected for MRI turbulence in the regions of the disk where the magnetic field is well coupled to the fluid.

We can readily compare the masses of planets that are predicted for each of these mechanisms of disk viscosity. For the case of a turbulent viscous disk (e.g. Bryden et al., 1999; Rafikov, 2002), the gap will open when the angular momentum transfer rate by tidal torque $`\dot{H_T}=0.23(M_p/M_{})^2\mathrm{\Sigma }a_p^4\mathrm{\Omega }^2(a_p/h)^3`$ exceeds that driven by viscous torque $`\dot{H_\nu }=3\pi \mathrm{\Sigma }\nu a_p\mathrm{\Omega }`$. Here $`\mathrm{\Omega }`$ is the Keplerian angular velocity and $`\nu `$ is the viscosity. Expressing the viscosity in terms of a standard $`\alpha `$ parameter $`\nu \alpha _{turb}c_sh`$ (Shakura & Sunyaev, 1973), the ratio of a gap opening planetary mass $`M_p`$ to a stellar mass $`M_{}`$ is

$$\frac{M_{p,turb}}{M_{}}\sqrt{40\alpha _{turb}\left(\frac{h}{a_p}\right)^5}.$$
(1)
where the disk height and the radius of the planet’s orbit are $`h`$ and $`a_p`$, respectively.

In the case that disk viscosity is provided solely by damped density waves, then the ratio of the planetary masses has been calculated by Rafikov (2002);

$$\frac{M_{p,damp}}{M_{}}\frac{2}{3}\left(\frac{h}{a_p}\right)^3\mathrm{min}[5.2Q^{5/7},3.8\left(Q\frac{a_p}{h}\right)^{5/13}]$$
(2)
where the Toomre parameter is $`Q=\mathrm{\Omega }c_s/(\pi G\mathrm{\Sigma })`$, $`c_s`$ is the sound speed, and $`\mathrm{\Sigma }`$ is the surface mass density. Note that the gap opening condition for an inviscid disk by Lin & Papaloizou (1985) corresponds to $`Q=1`$ case. For the standard disk model by Chiang et al. (2001), we calculate that $`Q=54`$ and $`h/a_p0.04`$ at 1 AU, so that the planetary mass is very small.

By taking the ratio of these two predicted planetary mass scales, we immediately see that the planetary masses in well coupled zones in which the MRI is active will be much larger than in the dead zone for which only the $`\alpha _{damp}`$ value pertains. As an example, at 1 AU, and for a Toomre Q value of 54 calculated above, the mass of a Jovian planet in a well-coupled zone dictated by MRI turbulence $`M_{p,turb}`$ will exceed the mass that a Jovian protoplanet will attain in a dead zone through density wave damping $`M_{p,damp}`$ provided that $`\alpha _{turb}`$ exceeds a minimum value of;

$$\alpha _{turb}>\alpha _{min}\mathrm{min}[1.0\times 10^3\frac{h}{a_p},7.4\times 10^3\left(\frac{h}{a_p}\right)^{23/13}]$$
(3)
which is less than $`10^4`$ for the aspect ratio at 1 AU ($`h/a_p0.04`$). Numerical simulations of the MRI instability show that $`\alpha _{turb}`$ is much larger than this – by two orders of magnitude or more. Thus, a distinct jump in planetary masses – by a factor that is similar to the ratio of the masses of the Jovian to terrestrial planets – is expected at the disk radius that determines the maximum radial extent of the dead zone.

In this paper, we perform a detailed investigation in order to calculate precisely where the dead zone occurs in well constrained protostellar disk models. This allows us to determine where in the disk that this “jump” in planetary masses is likely to occur. In our second paper (to be submitted), we calculate precisely what these planetary masses are. We proceed by carefully computing the ionization balance within the protostellar disk, since it is this that determines the coupling of the magnetic field to the disk and hence, where the MRI induced turbulence can be sustained. The main sources of ionization in a circumstellar disk are X-rays from the stellar plasma (probably produced by magnetospheric reconnection events) as well as cosmic rays from the outside of the stellar system. Glassgold et al. (2000) showed that X-ray ionization dominates cosmic ray ionization out to $`1000`$ AU in an optically thin disk. We compared these ionization sources for the self-consistent disk model developed by Chiang et al. (2001). We include radioactive elements as a global and constant ionization source. We used the disk model developed by Chiang et al. (2001) primarily because it reproduces the spectral energy distributions (SEDs) of two T-Tauri and three Herbig Ae stars extremely well.

Our major finding is that dead zones typically encompass the region where terrestrial planets are found in our solar system ($`2.03.3`$ AU). The physical definition of a dead zone is reviewed in §2. We introduce the disk model of Chiang et al. (2001) and elucidate the connection between the disk model and the X-ray ionization rate (using Glassgold et al., 1997) in §3. We state our results in §4 and apply them to the well studied case of AA Tau in §5. The important symbols used in this paper are summarized in Table 1.

## 2 Dead Zones in Disk Models

A dead zone is the region in a disk where MRI turbulence is unsustainable due to poor coupling of the field to the disk. The absence of turbulence is what gives this region its negligible viscosity. Gammie (1996) showed that MRI turbulence will be damped at a physical scale of $`\lambda `$ whenever the local growth rate of the MRI instability ($`V_A/\lambda `$) is balanced by the Ohmic diffusion at that scale ($`\eta /\lambda ^2`$). For complete damping at some radius of the disk, one requires that all the turbulence scales less than the local pressure scale height $`h(a)`$ are damped; i.e. $`\lambda h`$. By equating these two time-scales, one finds that the MRI turbulence damps in regimes with a small magnetic Reynolds number:

$$R_{eM}=\frac{V_Ah}{\eta }1.$$
(4)
Here, $`V_AB/\sqrt{4\pi \rho }\alpha _{turb}^{1/2}c_s`$ is the Alfvén speed, $`\rho `$ is the mass density, $`c_s=\sqrt{kT/\mu m_H}`$ is the sound speed ($`k`$ is the Boltzmann constant, $`T`$ is the temperature, $`\mu `$ is the mean molecular weight, and $`m_H`$ is the mass of the hydrogen) and $`\alpha _{turb}`$ is the viscosity parameter that measures the strength of the turbulence (Shakura & Sunyaev, 1973). This formula is one of the most important in the paper because it is the condition that determines the dead zone region in a disk once the disk structure and ionization state are specified. The importance of the latter is readily established by noting that the Reynolds number in this equation is sensitive to the ionization of the disk, through its dependence on the diffusivity of the magnetic field, $`\eta `$, which takes the form;

$$\eta =\frac{234}{x_e}T^{1/2}\mathrm{cm}^2\mathrm{s}^1,$$
(5)
where $`T`$ is the disk temperature that is obtained from disk models (either interior $`T_i`$ or disk surface $`T_{ds}`$ depending on the disk height we are interested in – see §2.3, eq. (3) and (4) in Chiang et al., 2001). The electron fraction is $`x_en_e/n`$ where $`n_e`$ and $`n`$ are the number density of electrons and neutral atoms respectively. The lower the electron fraction, the higher is the diffusivity coefficient for the field and hence the smaller is the magnetic Reynolds. Thus, poorly ionized disks will tend to have lower values of $`R_{eM}`$, resulting in larger dead zones.

One of the recent numerical MRI disk simulations by Fleming et al. (2000) defined the magnetic Reynolds number as

$$R_{eM}^{}=\frac{c_sh}{\eta },$$
(6)
and determined the critical value $`R_{eM,crit}^{}=100`$ assuming the presence of a uniform vertical magnetic field. Another recent analysis by Balbus & Terquem (2001) suggested that the Hall effect should be included into the diffusion equation. The Hall effect is due to the velocity difference between electrons and ions, and causes a small transverse potential difference. The effect also diffuses the magnetic field just as the ohmic diffusion does. It implies that the magnetic Reynolds number might be smaller than the case of Fleming et al. (2000).

Summarizing these results, it appears that the critical value of magnetic Reynolds number $`R_{eM,crit}`$ is somewhat uncertain and that there are a few different definitions for $`R_{eM}`$. Taking these uncertainties into account, we used the critical magnetic Reynolds number of 1 and 100 with the alpha parameter of 1, 0.1, and 0.01. The extreme cases of $`(R_{eM,crit},\alpha _{turb})=`$(1, 1), (1, 0.01), (100, 1), and (100, 0.01) correspond to $`R_{eM}^{}=1,10,100,\mathrm{and}1000`$ respectively in the definition by Fleming et al. (2000). The following is the criterion we used to determine the dead zone:

$$R_{eM}=\frac{V_Ah}{\eta }R_{eM,crit}.$$
(7)
The dead zone is the region in the disk where the magnetic Reynolds number $`R_{eM}`$ is smaller than some critical value $`R_{eM,crit}=1\mathrm{or}100`$. This region is characterized as having no MRI turbulence.

The calculation of the needed electron fraction $`x_e`$ (see eq. (5)) is well established, at least in principle. Assuming steady state ionization balance, Oppenheimer & Dalgarno (1974) derived an equation for the electron fraction in the gas; namely,

$$x_e^3+\frac{\beta _t}{\beta _d}x_Mx_e^2\frac{\zeta }{\beta _dn}x_e\frac{\zeta \beta _t}{\beta _d\beta _rn}x_M=0,$$
(8)
where $`x_M`$ is the metal fraction. In our case, $`\zeta `$ represents the ionization rate of X-rays, cosmic rays or radioactive elements, or total ionization rate of these three sources. Their mathematical forms are found in §3. Three beta terms represent three kinds of recombination rate coefficients – the dissociative recombination rate coefficient for electrons with molecular ions ($`\beta _d=2\times 10^6T^{1/2}\mathrm{cm}^3\mathrm{s}^1`$), the radiative recombination coefficient for electrons with metal ions ($`\beta _r=3\times 10^{11}T^{1/2}\mathrm{cm}^3\mathrm{s}^1`$) and the rate coefficient of charge transfer from molecular ions to metal atoms ($`\beta _t=3\times 10^9\mathrm{cm}^3\mathrm{s}^1`$).

In this paper, we only take account of the two extreme cases – the metal poor ($`x_M=0`$) and the metal dominant ($`x_Mx_e`$). In both cases, the electron fraction is written as

$$x_e=\sqrt{\frac{\zeta }{\beta n}},$$
(9)
where the number density $`n`$ depends on disk models (see eq. (17)) and $`\beta `$ is a recombination rate coefficient. In the metal poor case, only the recombination of electrons with molecular ions (dissociative recombination) is important and $`\beta =\beta _d=2\times 10^6T^{1/2}\mathrm{cm}^3\mathrm{s}^1`$. We use this coefficient throughout the paper except §5. In the metal dominant case, only the recombination of electrons with heavy metal ions (radiative recombination) becomes important and the recombination rate coefficient is replaced by the radiative recombination coefficient $`\beta =\beta _r=3\times 10^{11}T^{1/2}\mathrm{cm}^3\mathrm{s}^1`$ (Fromang et al., 2002). The astrophysical significance of these two cases is related to the density regime of the protostellar disk. We compare the results for these two cases as applied to AA Tau, in §5.

## 3 Ionization Rate of a Self-Consistent Disk Model

The previous section outlined precisely how the dead zone within a disk can be calculated once a disk model and ionization rate can be specified. In this section, we review the three sources of ionization for protostellar disks – X-rays from the central protostar, cosmic rays, and radioactive elements that are mixed with the disk material – and then calculate their respective ionization rates in the context of our chosen disk model. There is one caveat concerning cosmic ray ionization of protostellar disks – cosmic rays are likely to be swept away from an region containing a turbulent MHD jet or outflow (e.g. Skilling & Strong, 1976; Cesarsky & Völk, 1978). Since jets are manifest in Class I and II Young Stellar Objects (YSOs), it is possible that the bulk of the ionization of these sources is achieved solely by their YSO X-ray fluxes. For completeness, we show both sets of results – ionization from purely X-ray as well as combined ionization.

The combined contribution of each of these ionization sources produces the total ionization rate which is their sum;

$$\zeta _{tot}=\zeta _x+\zeta _{CR}+\zeta _{RA}.$$
(10)
We summarize the calculation of each of these three rates below and then also compare the effects of X-rays with that of cosmic rays in our disk model calculations.

### 3.1 Ionization Rates

The X-ray ionization rate was investigated by Krolik & Kallman (1983) and further developed by Glassgold et al. (1997). Glassgold et al. (2000) wrote the secondary electron contribution to the X-ray ionization rate $`\zeta _x`$ as follows:

$$\zeta _x=\left[\left(\frac{L_x}{kT_x4\pi d^2}\right)\sigma (kT_x)\right]\left(\frac{kT_x}{\mathrm{\Delta }ϵ}\right)J(\tau ,x_0),$$
(11)
where $`L_x`$ is the observed X-ray luminosity, and $`\sigma (kT_x)`$ and $`\tau (kT_x)`$ are the total photoelectric absorption cross section and the optical depth at the energy $`E=kT_x`$ respectively. The distance between the X-ray source and some point of the disk surface is denoted by “d” and the energy to make an ion pair is $`\mathrm{\Delta }ϵ`$. In above equation, the first factor between the square brackets corresponds to the primary ionization rate, assuming the same energy $`E=kT_x`$ for all primary electrons. The second factor $`kT_x/\mathrm{\Delta }ϵ`$ then reads as the number of secondary electrons produced by a photoelectron with the energy $`kT_x`$. The last factor $`J(\tau ,x_0)`$ represents the attenuation of X-rays. Using the dimensionless energy parameter $`x=E/kT_x`$, the attenuation factor $`J(\tau ,x_0)`$ is written as

$$J(\tau ,x_0)=_{x_0}^{\mathrm{}}x^ne^{x\tau \left(kT_x\right)x^n}𝑑x.$$
(12)
The optical depth $`\tau (kT_x)`$ measures how opaque the disk is toward the X-ray radiation from the star and is written as

$$\tau (kT_x)=N_H\sigma (kT_x),$$
(13)
where $`\sigma (kT_x)`$ is

$$\sigma (kT_x)=\stackrel{~}{\sigma }_0\left(\frac{kT_x}{keV}\right)^n.$$
(14)
We assumed that heavy elements are depleted onto grains and used $`\stackrel{~}{\sigma }_0=8.50\times 10^{23}\mathrm{cm}^2`$ and $`n=2.810`$ (Glassgold et al., 1997). The surface number density $`N_H`$ is measured along the radiation path from the X-ray source. Letting $`\alpha ^{}`$ be the angle between the radiation path and the radial axis, the surface number density is

$`N_H`$ $`=`$ $`{\displaystyle \frac{N_{}}{\mathrm{sin}\alpha ^{}}}`$ (15)
$`=`$ $`{\displaystyle \frac{_z^{\mathrm{}}n(a,z^{})𝑑z^{}}{\mathrm{sin}\alpha ^{}}},`$
where $`N_{}`$ is the vertical surface density and $`n(a,z)`$ is the number density at the disk radius “a” and the height “z”. Note that $`n(a,z)`$ and $`\alpha ^{}`$ depend on the disk model (see eq. (17) and (21) respectively).

Usually, the ionization rate by cosmic rays is estimated as $`\zeta _{CR,0}10^{17}\mathrm{s}^1`$ (Spitzer & Tomasko, 1968), but it’s not generally known to within better than an order of magnitude (Glassgold et al., 2000). In some cases, however, the cosmic-ray ionization rate is well constrained. As an example, van der Tak & van Dishoeck (2000) obtained the cosmic-ray ionization rate of $`2.6\pm 1.8\times 10^{17}\mathrm{s}^1`$ through $`\mathrm{H}^{13}\mathrm{CO}`$ submillimeter emission lines from massive protostars. This value is in good agreement with Voyager/Pioneer data ($`4\times 10^{17}\mathrm{s}^1`$). Since the attenuation length for cosmic rays is $`\chi _{CR}96\mathrm{g}\mathrm{cm}^2`$ (Umebayashi & Nakano, 1981), we follow Sano et al. (2000) and write the cosmic ray ionization rate as follows

$$\zeta _{CR}=\frac{\zeta _{CR,0}}{2}\left(\mathrm{exp}\left[\frac{\chi (a,z)}{\chi _{CR}}\right]+\mathrm{exp}\left[\frac{\mathrm{\Sigma }(a)\chi (a,z)}{\chi _{CR}}\right]\right),$$
(16)
where $`\chi (a,z)`$ is the vertical column density measured at some height $`z`$ ($`\chi (a,z)=_z^{\mathrm{}}\rho (a,z^{})𝑑z^{}`$).

The radioactive elements in the disk is yet another source of ionization, though they usually have only a minor effect on the total ionization rate. Following Umebayashi & Nakano (1981), we use the ionization rate by radioactive elements of $`\zeta _{RA}=6.9\times 10^{23}\mathrm{s}^1`$.

### 3.2 Disk Model

Chiang & Goldreich (1997), following Kenyon & Hartmann (1987), developed a self-consistent passive disk model. Their model has a two-layer structure – a high temperature surface layer and a lower temperature interior. Dust grains in the surface layer of the disk absorb the flux of ultra-violet (UV) photons from the central star. Half of the emission from the grains in the surface layer escapes into the space and the remaining half heats the disk interior up. The model assumes the vertical hydrostatic equilibrium so that the number density may be expressed as

$$n(a,z)=n_0(a)e^{\frac{z^2}{2h^2}},$$
(17)
where the number density at the midplane $`n_0`$ is

$$n_0(a)=\mathrm{\Sigma }(a)/(\sqrt{2\pi h(a)^2}\mu _g).$$
(18)
Here, $`\mu _g=3\times 10^{24}g`$ is the mean molecular weight of the gas, $`\mathrm{\Sigma }(a)`$ is the surface density and $`h(a)`$ is the pressure scale height (Chiang et al., 2001):

$`\mathrm{\Sigma }(a)`$ $`=`$ $`\mathrm{\Sigma }_0\left({\displaystyle \frac{a}{\mathrm{AU}}}\right)^{3/2}`$ (19)
$`h(a)`$ $`=`$ $`\left({\displaystyle \frac{T_i}{T_c}}\right)^{1/2}\left({\displaystyle \frac{a}{R_{}}}\right)^{1/2}a,`$ (20)
where $`\mathrm{\Sigma }_0`$ is the surface mass density at 1 AU, $`T_cGM_{}\mu _g/kR_{}`$ is a virial temperature that measures the gravitational potential at the surface of the central star, and $`M_{}`$ and $`R_{}`$ are the stellar mass and radius respectively. Note that the disk interior temperature $`T_i`$ (the disk temperature below the pressure scale height “h”) decreases with a disk radius, so $`h(a)`$ increases slightly weaker than the power of three halves (see eq. (20)) and $`n_0(a)`$ decreases slightly weaker than the power of three (see eq. (18)). The pressure scale height $`h(a)`$ is the height measured from the midplane to the interior disk, not the surface layer of the disk.

The ionization of this disk by cosmic rays and radioactive elements is straight-forward to calculate. The X-ray ionization requires that we model the geometry of the source of the X-rays – the central young stellar object. X-ray emission is thought to arise from the reconnection of large-scale magnetized loops, that are possibly but not necessarily related to the magnetopause radius of the stellar magnetic field and the inner edge of the disk. Accordingly, we place the X-ray source in this picture at a fiducial distance of $`12R_{}`$ above the disk and $`12R_{}`$ away from the central star following Glassgold et al. (1997).

The X-rays from this source graze the surface of the flaring disk and penetrate it to produce the X-ray induced ionization. It is convenient to relate this grazing angle, $`\alpha ^{}`$, to the grazing angle defined in the disk model of Chiang et al. (2001) $`\alpha `$ (see eq. (5) in Chiang et al., 2001). Assuming that the average radiation from the central star originates at a distance $`R_{}/2`$ from the disk plane on the stellar surface (Hartmann, 2000), we write the grazing angle from the X-ray source (see eq. (15)) as

$`\alpha ^{}=\alpha \beta +\gamma `$ (21)
$`\beta \mathrm{tan}^1\left({\displaystyle \frac{d\mathrm{lnH}}{d\mathrm{lna}}}{\displaystyle \frac{H}{a}}\right)`$
$`\gamma \mathrm{tan}^1\left({\displaystyle \frac{H\frac{1}{2}R_{}}{a\frac{\sqrt{3}}{2}R_{}}}\right)\mathrm{tan}^1\left({\displaystyle \frac{H12R_{}}{a12R_{}}}\right),`$
where $`H`$ is the height of the surface layer of the disk and defined as $`H/h4`$ in Chiang et al. (2001). Fig. 1 is a schematic figure of this equation. The grazing angle from the X-ray source, $`\alpha ^{}`$ is the angle between the disk midplane and the solid arrow. The first term of the righthand side of the equation, $`\alpha `$ is the angle between the disk surface and the dashed arrow of the UV emission from the star. The second term $`\beta `$ shows the flaring angle of the disk and the third term $`\gamma `$ comes from the geometry of the X-ray source. Note that $`\gamma =0`$ corresponds to an X-ray source that is at the stellar surface.

We calculated the X-ray ionization rate (eq. (11)) by integrating the attenuation factor (eq. (12)) using MATLAB. Instead of an infinite upper limit of the energy, we used $`x=100`$ for the upper limit of the integration. For the lower limit, we used $`x=1\mathrm{keV}/kT_x`$ ($`E=1\mathrm{keV}`$), following Glassgold et al. (1997). The surface number density along the path of the radiation $`N_H`$ found in the optical depth (eq. (13)) is calculated by using the grazing angle $`\alpha ^{}`$ (eq. (21)) and the number density $`n(a,z)`$ that is obtained from the disk model of Chiang et al. (2001) (see eq. (17) and (18)). Note that this number density $`n(a,z)`$ depends on the disk interior temperature $`T_i`$ and the pressure scale height of the disk $`h`$. Both of them are determined by solving the radiative equilibrium equations (3) and (4) in Chiang et al. (2001).

## 4 Results

In this section, we study the extent of dead zones for various disks. We show that X-rays and cosmic rays are comparable ionization sources and determine the size of dead zones from the total ionization rate. Throughout this section, we assume the metal poor disk (set $`\beta =\beta _d`$ in eq. (9)).

### 4.1 Standard Disk Model

First, we show the results of a fiducial model which uses the standard disk model of Chiang et al. (2001), the critical magnetic Reynolds number $`R_{eM,crit}=1`$ with $`\alpha _{turb}=1,0.1,0.01`$, and when X-rays exist, the X-ray luminosity of $`L_x=10^{29}\mathrm{erg}\mathrm{s}^1`$ with the temperature of $`kT_x=1`$ keV. We chose the standard X-ray luminosity of $`L_x=10^{29}\mathrm{erg}\mathrm{s}^1`$ because typical young stars have X-ray luminosities $`L_x10^{28}10^{30}\mathrm{erg}\mathrm{s}^1`$. Though the standard disk of Chiang et al. (2001) has the surface mass density at 1 AU of $`\mathrm{\Sigma }_0=10^3\mathrm{g}\mathrm{cm}^2`$, we don’t call this minimum mass solar nebular model because these two models are qualitatively different from each other. With the standard disk model of Chiang et al. (2001), we calculated the surface number density (eq. (15)) and/or the surface mass density ($`\chi `$ in eq. (16)) to estimate the X-ray and/or cosmic ray ionization rate at each radius and height of the 2D disk. All disk parameters are the same as in Table 1 of Chiang et al. (2001). Other parameters used in this subsection is summarized in Table 2 in this paper. We determined the critical height of the dead zone at each radius of the disk from eq. (7).

The size of dead zones arising solely from X-ray irradiation are plotted in Fig. 2. The figure shows the vertical cross section of the standard disk. The region below the dashed, long-dashed and dot-dashed lines show the dead zones. As the magnetic field becomes weaker (the alpha value changes $`\alpha _{turb}=1,0.1,0.01`$), the dead zone becomes larger. The result follows from the fact that as the magnetic field becomes weak, the growth rate of the MRI instability is reduced with respect to the local damping rate (see eq. (4) and the text therein).

The equivalent result for cosmic rays is shown in Fig. 3. In this case, we just replaced the ionization rate in electron fraction (eq. (9)) with the cosmic ray ionization rate calculated by eq. (16). Again, we find the smaller dead zone as the magnetic field becomes stronger.

These two figures show that cosmic rays dominate 1 keV X-rays in ionizing a disk. The dead zone stretches out to $`57`$ AU in the case of illumination by X-rays and only to $`12`$ AU in cosmic ray case.

The dividing line between X-ray or cosmic ray ionized regions may be found by equating their ionization rates. We call this resulting vertical scale height the Röntgen height. Note that due to the higher number density, the X-ray ionization rate decreases toward the disk midplane as does the electron fraction, while the cosmic ray ionization rate is about the same everywhere in the disk. X-ray ionization dominates cosmic ray ionization above the Röntgen height.

Fig. 4 shows the dead zones estimated by the total ionization rate. Again, we took a surface mass density at 1 AU of $`\mathrm{\Sigma }_0=10^3\mathrm{g}\mathrm{cm}^2`$, the dead zone criterion $`R_{eM}1`$, an X-ray luminosity of $`L_x=10^{29}\mathrm{erg}\mathrm{s}^1`$, and an X-ray energy of $`kT_x=1`$ keV. The Röntgen height suggests that X-rays dominate cosmic rays only in the surface layer – where the disk is optically thin. Comparing Fig. 3 with 4, it is apparent that the size of dead zones are determined by cosmic rays for the standard disk case.

There are two characteristic dead zone radii. One of them is the midplane radius where the dead zone surface cuts the disk midplane. The other one is the fiducial radius where the dead zone surface crosses the pressure scale height of the interior disk. We take this latter case to be our “fiducial” value because planets are likely to form in the interior disk that has much higher density (the order of $`10^3`$) and larger particles ($`0.011000\mu m`$ instead of $`0.011\mu m`$) compared to the surface layers.

### 4.2 Sensitivity of Dead Zone Radius to Model Parameters

Next, we compare cases of different parameters ($`\mathrm{\Sigma }_0,kT_x`$) with one another. This is done to reveal which parameter has the largest effects on the size of dead zone. We also present results for two cases; one for only X-ray induced ionization, and the other for ionization arising from the total ionization rate. We make this distinction in order that the contribution of X-rays can be clearly discerned (cf. discussion at the beginning of §5). Throughout this subsection, we set $`R_{eM,crit}=1`$ and $`\alpha _{turb}=0.01`$ (see Table 2). These are the conditions suggested by, for example, Gammie (1996).

Fig. 5 is the plot of the fiducial dead zone radius estimated by X-ray alone as a function of the X-ray energy ($`kT_x=1,2,3,5,\mathrm{and}10`$ keV). For the X-ray ionization calculation, we used the typical X-ray luminosity of $`L_x=10^{29}\mathrm{erg}\mathrm{s}^1`$. We include the minimum mass solar nebula ($`\mathrm{\Sigma }_0=10^3\mathrm{g}\mathrm{cm}^2`$) as well as a very massive disk of $`\mathrm{\Sigma }_0=10^5\mathrm{g}\mathrm{cm}^2`$ suggested by Murray et al. (1998) for Jupiter to migrate from 5 AU to $`<1`$ AU. It is apparent that the sizes of dead zones decrease as the X-ray energy increases and/or the disk surface mass density decreases.

Note that the figure depends on the choice of the critical magnetic Reynolds number $`R_{eM,crit}`$, the alpha parameter $`\alpha _{turb}`$, and the X-ray luminosity $`L_x`$. For example, if we increase the X-ray luminosity by two orders of magnitude ($`L_x=10^{31}\mathrm{erg}\mathrm{s}^1`$) in Fig. 5, then the resulting curves are pushed down by the factor of 10. This is because of the definition of the magnetic Reynolds number $`R_{eM}=\sqrt{\alpha _{turb}\zeta }[c_sh/(234\sqrt{\beta nT})]`$. Here, the bracket term is determined by a disk model (i.e. independent of the X-ray energy or luminosity) so that the magnetic Reynolds number is changed depending on $`\alpha _{turb}`$ and $`\zeta `$. For X-rays, the ionization rate is written as $`\zeta _x=L_x[\sigma (kT_x)/(4\pi d^2\mathrm{\Delta }ϵ)J(\tau ,x_0)]`$ where the luminosity changes the result by two orders of magnitude in our case (we used $`L_x=10^{29}\mathrm{erg}\mathrm{s}^1`$ and $`10^{31}\mathrm{erg}\mathrm{s}^1`$).

The results for the dead zone radius for a disk undergoing the total ionization rate is shown in Fig. 6. It shows that the radii of dead zones are almost constant for X-ray energy of $`kT_x=13\mathrm{keV}`$ and coincide with the results of Fig. 5 in the higher X-ray energy region.

These figures also show the obvious point that the size of dead zones is not affected by low energy X-rays ($`kT_x13\mathrm{keV}`$) and that the dead zone gets smaller as the surface mass density decreases. The decrease of the surface mass density (from upper lines to lower ones) leads to a smaller optical depth $`\tau `$ and a larger ionization rate $`\zeta `$, so that the dead zone becomes smaller.

We summarize the effect of changing parameters. If the X-ray luminosity is larger ($`10^{31}\mathrm{erg}\mathrm{s}^1`$ rather than $`10^{29}\mathrm{erg}\mathrm{s}^1`$), X-rays penetrate deeper in the disk and therefore the dead zone becomes smaller. If the critical magnetic Reynolds number is larger (100 rather than 1), a very small magnetic field diffusivity can destroy the turbulence – the viscosity becomes negligible and therefore the dead zone tends to be larger. If the alpha parameter is smaller (0.01 rather than 1), then the magnetic field is weaker, turbulent viscosity is smaller and the dead zone becomes larger. In short, to have a smaller dead zone, a disk has to have a smaller surface mass density $`\mathrm{\Sigma }_0`$, a larger X-ray luminosity $`L_x`$ and a larger $`\alpha _{turb}`$ value as well as a smaller critical Reynolds number $`R_{eM,crit}`$.

To draw the reader’s attention to our major result, we present Fig. 7 and 8 and study the effect of the surface mass density at 1 AU $`\mathrm{\Sigma }_0`$. Fig. 7 uses the same data sets as Fig. 5, but here we plot the fiducial dead zone radius as a function of the surface mass density at 1 AU ($`\mathrm{\Sigma }_0=10^3,5\times 10^3,10^4,6\times 10^4,\mathrm{and}10^5\mathrm{g}\mathrm{cm}^2`$) instead of the X-ray energy. It is apparent that if cosmic rays penetrate disks without being swept away, they dominate X-rays in ionizing a disk except for $`kT_x510\mathrm{keV}`$.

Fig. 8 shows the equivalent result for the total ionization rate and uses the same data sets as Fig. 6. The X-rays have a rather small effect on the size of dead zones even when they have high energy ($`kT_x510\mathrm{keV}`$). The most striking result is that dead zone radii vary by at most an order of magnitude, from $`110`$ AU for a range of two orders of magnitude in disk surface density and an order of magnitude in the X-ray energy. The robustness of this result indicates that the region of Jovian-mass planet formation is rather similar for a wide variety of protostellar disk systems.

We can see that the differences between lines in Fig. 7, 8 are smaller than those in Fig. 5, 6. It appears that the size of the dead zone is most sensitive to the surface mass density and, to a lesser extent, the X-ray energy. These figures also show that the dead zone produced by cosmic rays is usually smaller than that by X-rays. X-rays dominate cosmic rays in ionizing the circumstellar disk only when the X-ray source has a high energy ($`kT_x510\mathrm{keV}`$).

Table 3 shows the size range of dead zones for four extreme cases – $`(R_{eM,crit},\alpha _{turb})=(1,1)`$, (1, 0.01), (100, 1), and (100, 0.01). The minimum and maximum dead zone radii correspond to the minimum and maximum surface mass density at 1 AU. It is apparent that the dead zone is sensitive not only to the surface mass density or X-ray energy, but also to the value of the magnetic Reynolds number. The important point of this table is that the dead zone radii for total ionization rate are about the same as the region of terrestrial planets in our solar system except the one extreme case. For this extreme case of $`(R_{eM,crit},\alpha _{turb})=(100,0.01)`$, which is equivalent to $`R_{eM,crit}^{}=1000`$ or $`(R_{eM,crit},\alpha _{turb})=(1,10^6)`$, the dead zone gets pushed out to 5.8 AU.

Finally, we note that the data presented in Figures 5, 7 and 8 corresponding to the variation of dead zone radii can be well fit by approximate linear curves. For completeness, we present these fits in Appendix A - Table 5, 6, 7, and 8.

## 5 Application to the case of AA Tau

We now specialize our model to an observed protostellar disk – AA Tau. We consider this source because of the excellent fit of the Chiang et al. model to the data from this system. The stellar and disk parameters that we used are listed in Table 4 and are taken from Chiang et al. (2001) other than X-ray parameters. The X-ray luminosity $`L_x`$ and the X-ray energy $`kT_x`$ are taken from Neuhäuser et al. (1995). For the dead zone criterion, we used $`R_{eM,crit}=1`$ with $`\alpha _{turb}=1,0.1,0.01`$. These parameters are summarized in Table 2.

Fig. 9 is the estimated dead zone for AA Tau. It shows that X-ray ionization dominates cosmic ray ionization only in the surface layer where the optical depth is very small (see Röntgen height). The dead zone becomes smaller as the magnetic field gets stronger because of the stronger viscous torque that it exerts. In all cases the dead zones by cosmic rays are smaller than those produced by X-rays.

If the disk of AA Tau is metal dominant, the dead zone gets smaller. Fig. 10 shows the dead zones estimated by total ionization rate for both metal poor and metal dominant cases. Since metal ions recombine with electrons slowly, the electron fraction for metal dominant case is larger (compare $`\beta _d`$ with $`\beta _r`$) – $`R_{eM}`$ tends to be larger – and hence the dead zone is smaller. Note that the recombination coefficient of metals is $`5`$ orders of magnitude smaller than that of molecular ions. The “real” dead zone would be calculated by solving eq. (8) directly and would be located somewhere between these two extremes.

## 6 Discussion and Conclusions

We calculated the dead zones for a variety of disk models, as well as the specific case of a Class II source, AA Tau. This kind of calculation has been done by several authors (e.g., Gammie, 1996; Igea & Glassgold, 1999; Sano et al., 2000; Fromang et al., 2002), but all of them use either the minimum mass solar nebula model developed by Hayashi et al. (1985) (e.g., Gammie, 1996; Igea & Glassgold, 1999; Sano et al., 2000) or the $`\alpha `$ disk model developed by Shakura & Sunyaev (1973) (Fromang et al., 2002). All our disk models are obtained by using the self-consistent passive disk model of Chiang et al. (2001). To calculate X-ray ionization rates, we followed the method used by Glassgold et al. (2000), while for cosmic rays, we used the method adopted by Sano et al. (2000). When the ionization rate becomes sufficiently low so that the magnetic Reynolds number is less than some critical value (1 or 100), we determined the size of the dead zone. We compared the result by X-rays with cosmic rays and also obtained the size of the dead zone determined by the total ionization rate.

The fact that our dead zones extend as far as the region of terrestrial planet formation in our solar system is very important. In the accretion picture, considerable gas accretion must occur onto a sufficiently massive protoplanetary core. This must inevitably take place in the turbulent region of the disk beyond the dead zone. The major implication of our work is that the division between Jovian and sub-Jovian or terrestrial planets may be due to the presence of a dead zone. We explore the significance of these results for predicting planetary masses in a forthcoming paper. Our basic results are as follows:

1. Our major finding is that the typical dead zone encompasses, in physical scale, the terrestrial planets in our solar system.

A typical dead zone size estimated from the total ionization rate is $`0.242.7`$ AU for the fiducial disk. For an extreme case of $`(R_{eM,crit},\alpha _{turb})=(100,0.01)`$ ($`R_{eM,crit}^{}=1000`$ or $`(R_{eM,crit},\alpha _{turb})=(1,10^6)`$ equivalently), the dead zone extends out to $`5.8`$ AU. These results suggest that the observed exosolar Jovian-mass planets must have migrated to their observed positions from points of origin farther out in disk radius, beyond these dead zone radii.

2. X-rays dominate cosmic rays in ionizing a disk only in the surface of the disk for most cases.

Glassgold et al. (2000) noted that X-ray ionization dominates cosmic ray ionization out to $`1000`$ AU. We found that the disks are usually too optically thick for this to be generally true. X-rays only dominate cosmic rays in the surface layer of the disk. Cosmic rays determine the size of the dead zone under typical conditions (compare, for example, Fig. 3 with Fig. 4).

3. Sufficiently high energy X-rays could dominate cosmic rays in ionizing a disk.

The X-rays could dominate the cosmic rays if the X-ray energy is high $`kT_x=510\mathrm{k}\mathrm{e}\mathrm{V}`$ (compare Fig. 7 with 8). This is much higher energy for most observed sources however.

4. The size of a dead zone is sensitive to the disk surface mass density $`\mathrm{\Sigma }_0`$ and the X-ray energy $`kT_x`$.

We found that the size of the dead zone is most sensitive to the surface mass density, and to a lesser extent, the X-ray energy (see Fig. 5 & 6 and Fig. 7 & 8). The X-ray luminosity, the critical magnetic Reynolds number, the alpha viscous parameter however, have a significant effect on the size of a dead zone.

5. There is a power law relation between the size of a dead zone estimated by X-rays alone and $`\mathrm{\Sigma }_0`$ or $`kT_x`$ (see Appendix).

The dead zone radius has power-law relations with both the surface mass density and the X-ray energy. The dead zone radii estimated by X-rays obey $`a_d\mathrm{\Sigma }_0^{0.44}\mathrm{\Sigma }_0^{0.59}`$ and $`a_dkT_x^{0.46}kT_x^{1.1}`$, while those estimated by cosmic rays obey $`a_d\mathrm{\Sigma }_0^{0.28}\mathrm{\Sigma }_0^{0.49}`$. The dead zone radii determined by the total ionization rate has the relation of $`a_d\mathrm{\Sigma }_0^{0.28}\mathrm{\Sigma }_0^{0.59}`$.

We thank Norm Murray, Joe Weingartuer and Steve Balbus for stimulating conversations on these topics. We also thank an anonymous referee for useful review of the manuscript. SM is supported by McMaster University, while REP is supported by a grant from the National Science and Engineering Research Council of Canada (NSERC).

## Appendix A

## X-rays and Cosmic rays

We supply tables of X-ray and cosmic ray power-law relations separately for interested readers. These approximate relations are all obtained by using the least square fitting.

Table 5 shows the fitting parameters for the dead zone radius by X-rays as a function of X-ray energy $`kT_x`$. For the same choice of the critical magnetic Reynolds number $`R_{eM,crit}`$ and the alpha parameter $`\alpha _{turb}`$, it is apparent that the curves show roughly the same power. One of the examples is shown in Fig. 5.

Table 6 shows the fitting parameters for the dead zone radius by X-rays as a function of the surface mass density at 1 AU, $`\mathrm{\Sigma }_0`$. Again, we can see the similar power for the same choice of the critical magnetic Reynolds number $`R_{eM,crit}`$ and the alpha parameter $`\alpha _{turb}`$. One of the examples is shown in Fig. 7.

Table 7 shows the fitting parameters for the dead zone radius by cosmic rays as a function of the surface mass density at 1 AU, $`\mathrm{\Sigma }_0`$. One of the examples is shown in Fig. 7 (the gray solid line with crosses).

Table 8 shows the fitting parameters for the dead zone radius estimated by the total ionization rate as a function of the surface mass density at 1 AU, $`\mathrm{\Sigma }_0`$. Again, we can see the similar power for the same choice of the critical magnetic Reynolds number $`R_{eM,crit}`$ and the alpha parameter $`\alpha _{turb}`$. One of the examples is shown in Fig. 8. In particular, $`kT_x=13`$ keV give a very similar power law relation.


format: text

text:
# Theoretical analysis of spectral gain in a THz quantum cascade laser: prospects for gain at 1 THz

Preprint, in press APL

( )

## Abstract

In a recent Letter \[Appl. Phys. Lett. 82, 1015 (2003)\], Williams et al. reported the development of a terahertz quantum cascade laser operating at 3.4 THz or 14.2 meV. We have calculated and analyzed the gain spectra of the quantum cascade structure described in their work, and in addition to gain at the reported lasing energy of $``$ 14 meV, we have discovered substantial gain at a much lower energy of around 5 meV or just over 1 THz. This suggests an avenue for the development of a terahertz laser at this lower energy, or of a two-color terahertz laser.

Following on the initial successful development of quantum cascade lasers in the mid-infrared,FaiGma  the field has continued its advance with the recent appearance of terahertz (far-infrared) quantum cascade lasers. The initial examplesKohRoc  of these terahertz lasers, operating at $`4.5`$ THz, were based on chirped superlattice structures. An alternative design, operating at $`3.4`$ THz, and based on a simple four well (per period) structure was recently reportedWil03  by Williams et al. The simplicity of this structure provides a testbed for a nonequilibrium Green’s function theory we have recently developedLee02a ; Wac02b  to determine the nonequilibrium stationary state of quantum cascade laser structures operating under an applied voltage. This theory enables us to analyze transport and gain properties of these structures, giving us a tool to evaluate current-voltage (I-V) characteristics and nonequilibrium distribution functions, and to estimate level populations and lifetimes. In addition, we can use these parameters as a basis for evaluating and analyzing gain and absorption spectra in these structures. We have applied this theory to the structure reportedWil03  by Williams et al., and we report here the results of our theoretical analysis. We highlight especially the new finding of substantial gain at the low energy of $`5`$ meV, well below the spectral region considered in Ref. Wil03, .

Figure 1 shows the conduction band profile of two periods of the structure reported in Ref. Wil03,  with an applied bias of 64 mV/period or 12.2 kV/cm. The Wannier-Stark wavefunctions (modulus-squared) are also shown. The energetic positions of the wavefunctions shown in the figure are the Wannier-Stark energies renormalized by the mean-field due to electron-electron scattering.Lee02a  Following Ref. Wil03, , we designate levels 4 and 5 as the lower and upper laser levels respectively. In addition, we label level 1 (and 1’ for the neighboring period) the lower collector level, and level 2 (and 2’) the upper collector level.

We calculatedLee02a  the gain spectra by considering a simple two-level model interacting with a classical time-dependent electric field which gives the gain coefficientHau93

$$g(\omega )=\frac{\omega \pi }{cn_Bϵ_0L_p}\underset{\genfrac{}{}{0pt}{}{ij}{(E_j>E_i)}}{}|d_{ij}|^2(n_in_j)_{ij}(\omega ),$$
(1)
where we sum over contributions from each pair of levels $`i`$ and $`j`$. $`L_p=52.4`$ nm is the length of one period of the structure, and $`n_B=\sqrt{13}`$ is the background refractive index. The Lorentzian $`_{ij}(\omega )=(\mathrm{\Gamma }_{ij}/2\pi )/[(\mathrm{}\omega \mathrm{\Delta }E_{ji})^2+(\mathrm{\Gamma }_{ij}/2)^2]`$, with $`\mathrm{\Gamma }_{ij}=\mathrm{\Gamma }_i+\mathrm{\Gamma }_j`$, and $`\mathrm{\Delta }E_{ji}=E_jE_i`$. The populations, e.g., $`n_i`$, and the broadening parameters, e.g., $`\mathrm{\Gamma }_i`$, for level $`i`$ are extracted from the nonequilibrium Green’s functions as described in Ref. Lee02a, . The populations and broadening parameters, as well as the energy differences $`\mathrm{\Delta }E_{ji}`$ and dipole matrix elements $`d_{ij}=ez_{ij}`$ relevant for the following discussion are given in Table 1. In addition to the shift due to the mean-field potential, the energy differences given in Table 1 include also the renormalization due to electron-phonon and interface roughness scattering.Lee02a  We neglect impurity scattering since the doping density per period, $`n_e=2.8\times 10^{10}`$ cm<sup>-2</sup>, for this structure is low.

Figure 2 shows the calculated gain spectra at 30 K for three applied voltages, 64, 66, and 68 mV/period. At each of the two lower voltages, two strong gain features are seen: one lying between 10 and 15 meV, the other between 5 and 8 meV. There is a blue shift of the gain features as the bias increases in agreement with the experimental data.Wil03  At 68 mV/period, the two gain features merge giving a single broad gain feature stretching from $``$ 5 to 17 meV. To see the origin of these features, we look at the contributions to the spectra from transitions between each pair of levels.

Figure 3(a) shows the main transitions contributing to the gain for the applied voltage 64 mV/period. Contrary to the assignment given in Ref. Wil03, , we find that the strong gain feature at around 12.6 meV is due to the transition between the lower collector level (1’) and the lower laser level (4), and not between the upper and lower laser levels. In fact, at this bias (see table 1), the population in the lower laser level (4) is still slightly larger than the population in the upper laser level (5), giving rise to absorption rather than gain between these two levels. The dipole matrix element between levels 1’ and 4 ($`|z_{1^{}4}|=3.2`$ nm) is of the same order, if slightly smaller than that between levels 4 and 5 ($`|z_{54}|=5.14`$ nm). As importantly, however, there is a large population in the lower collector level 1’ (in fact, the most highly occupied level) and this gives rise to a substantial population inversion between levels 1’ and 4.

The strong gain feature at around 5.5 meV arises from the transition between the upper collector level (2’) and the upper laser level (5). This transition is favored by the large dipole matrix element ($`|z_{2^{}5}|=6`$ nm), larger than that between levels 4 and 5. The large population in 2’ again gives rise to a significant population inversion.

Figure 3(b) shows the gain contributions at 66 mV/period. As in the previous case, the main contributions to the gain arise again from the transitions between levels 1’ and 4 (13.6 meV), and 2’ and 5 (6 meV). In addition, there are two smaller gain peaks. One peak at 15.2 meV can be assigned to the transition between levels 4 and 5, i.e., the lasing transition considered in Ref. Wil03, . The other small peak at 7.8 meV is assigned to the transition between the collector levels 1’ and 2’. At the highest bias (68 mV/period) in Fig. 2, the contributions from these two latter transitions (4 – 5 and 1’ – 2’) become stronger leading to the double peaked structure seen at $`15`$ meV, and the shoulder at $`8`$ meV. The blue shift of the overall gain features with increasing bias can therefore be attributed to the appearance of these secondary peaks, as well as to the Stark shift mentioned in Ref. Wil03, .

In addition to the gain features discussed above, the calculations also show a strong gain feature at around 2 meV at the lower voltage of 48 mV/period. This feature disappeared as the bias was increased. However, its reappearance is seen again in the low energy region of the 68 mV/period spectrum in Fig. 2. At a higher voltage of 70 mV/period this feature (originating from the 1’ – 5 transition) becomes stronger.

Thus, a theoretical analysis of this structure shows a complex behavior of the gain spectra, and their dependence on the applied voltage. Several transitions contribute to the gain spectra, with the spectral position and strength of each contribution depending sensitively on the applied bias. This suggests the possibility of generating gain in a spectral region stretching from $`2`$ – 17 meV, with, in addition, some ability to tune the wavelengths at which the gain is enhanced by varying the applied bias.

We discuss next the robustness of our results to changes in the parameters used in the calculation. The values we have used for the broadening parameters $`\mathrm{\Gamma }_i`$ tend to overestimate the intersubband relaxation rate because $`\mathrm{\Gamma }_i`$ also includes the effect of intrasubband scattering processes. We have, however, also used an alternative approachWac02b  to calculate the gain spectra, in which the linear response of the Green’s functions and self-energies to an optical field are evaluated directly. This approach, which does not use the Lorentzian lineshape function, and hence does not require the parameters $`\mathrm{\Gamma }_i`$, reproduces the same gain features, in particular, the strong peaks around 6 meV, and between 10 – 15 meV. We have also repeated the calculations at 5 K, and including one extra level (i.e, considering six levels instead of five), and the main results of the above analysis are reproduced. The band structure parameters we used to describe the superlattice potential \[conduction band offset (66%): $`\mathrm{\Delta }E_c=0.15`$ eV; electron effective masses: $`m_e^w=0.067`$ in well, $`m_e^b=0.07945`$ in barrier\] are extracted from Ref. Vur01, . The resulting theoretical current densities are approximately double the experimental values reported in Ref. Wil03, . If, however, we use an 80% conduction band offset ($`\mathrm{\Delta }E_c=0.1825`$ eV) as reported in Ref. Wil03, , the calculated current densities are reduced, coming in close agreement with the experimental measurements. The main gain features discussed above are still present but slightly red-shifted, e.g., for 0.064 mV/period, the 5.5 meV feature in Fig. 3(a) is shifted down to 4 meV ($``$ 1 THz). The height of this gain feature increases by $`50`$%, while the height of the higher frequency gain feature decreases by $`25`$%.

We have reported here a detailed analysis of the gain spectra of a THz quantum cascade structure described in a recent Letter.Wil03  The gain spectra exhibits a complicated behavior and dependence on the applied voltage, with several different transitions contributing to the gain. Surprisingly, besides the main lasing transition designated in Ref. Wil03, , there are strong contributions to the gain from neighboring transitions, particularly transitions between the collector and upper or lower laser levels. These contributions are in fact larger than that of the designed lasing transition. The origin of these additional contributions to the gain is a sufficiently large overlap of the wavefunctions in the collector levels with those of the upper or lower laser levels to give dipole matrix elements comparable in magnitude to that of the lasing transition, as well as the large accumulation of population in the collector levels giving rise to a substantial population inversion with respect to the upper or lower laser levels. A notable finding is the appearance of a strong gain feature at around 5 meV, suggesting that this cascade structure design could also form the basis for a laser operating at around 1.2 THz ($`\lambda 250`$ $`\mu `$m). There was no consideration or mention of this long wavelength spectral region in Ref. Wil03, , which focused on the lasing emission around 3.4 THz. Thus, the results we report here urge more experimental investigations on this interesting structure, both for the prospect of a THz laser at a longer wavelength than achieved to date for any other quantum cascade laser structure, and as a test and verification of our theoretical approach and its predictive power.

This work is supported by the Deutsche Forschungsgemeinschaft.


text:
# Particle decay branching ratios for states of astrophysical importance in 19Ne

## I Introduction

States in 19Ne corresponding to 15O+$`\alpha `$ (Qαγ=3529.4 keV) and 18F+p (Q=6411.2 keV) resonances play an important role in explosive hydrogen burning. In energetic events known as x-ray bursts, 15O($`\alpha `$,$`\gamma `$)19Ne has long been thought to be a possible pathway from the hotCNO energy-generation bottleneck into the more exoergic hydrogen-burning processes above A=19, such as the NeNa cycle and the rp process Wiescher et al. (1999). In novae, the $`\beta `$+ decay of 18F is expected to be the prime source of the 511 keV $`\gamma `$-ray line, and of the sub-511-keV continuum $`\gamma `$-ray flux, observable by orbiting $`\gamma `$-ray observatories such as INTEGRAL Hernanz et al. (1999). The final yield of 18F and its associated radiation is sensitive to the reaction rates of 18F(p,$`\gamma `$)19Ne and 18F(p,$`\alpha `$)15O Coc et al. (2000); Iliadis et al. (2002).

The level density at relevant energies above the proton and alpha thresholds is low, so that the nuclear reactions mostly take place through isolated and narrow resonances. The temperature-dependent resonant reaction rates may be calculated by summing over resonances Rolfs and Rodney (1988):

$$<\sigma v>=\underset{R}{}\left(\frac{2\pi }{\mu kT}\right)^{3/2}\mathrm{}^2\left(\omega \gamma \right)_Rexp\left(\frac{E_R}{kT}\right),$$
(1)
where $`\mu `$ is the reduced mass, $`k`$ is Boltzmann’s constant, $`T`$ is the temperature, $`E_R`$ is the resonance energy and $`\omega \gamma `$ is its strength. The first factor in the resonance strength, $`\omega `$, is a statistical factor that depends on the spins in the incoming and outgoing channel, while $`\gamma `$ depends on the width and branching ratios for the resonance:

$$\gamma =\frac{\mathrm{\Gamma }_a}{\mathrm{\Gamma }}\frac{\mathrm{\Gamma }_b}{\mathrm{\Gamma }}\mathrm{\Gamma },$$
(2)
where $`a`$ and $`b`$ denote the incoming and outgoing particles ($`p`$, $`\alpha `$ or $`\gamma `$) for the reaction $`I(a,b)F`$. Thus, the reaction rate has a linear dependence on the branching ratios.

## II Experiment

Previous experiments have measured astrophysically-important branching ratios Magnus et al. (1990); Utku et al. (1998); Laird et al. (2002); Rehm et al. (2003); Davids et al. (2003) and resonance strengths Graulich et al. (2000); Rehm et al. (1996); Bardayan et al. (2001, 2002) for 19Ne. The goal of the present experiment was to perform an independent measurement of these branching ratios, and to achieve higher sensitivity in that measurement by using large-area silicon strip detectors Davinson et al. (2000). These detectors provided a large solid angle and an angular resolution in the lab of approximately 2 Visser (2003).

The 19F(3He,$`t`$)19Ne reaction was studied using 25 MeV 3He beams provided by the WNSL ESTU tandem Van de Graaf accelerator. The beam, with a typical current of 19 pnA, was incident on 80 $`\mu `$g/cm2 of CaF2 deposited on a 10 $`\mu `$g/cm2 carbon foil. The Enge spectrograph was placed at $`\theta =0^{}`$ with an aperture solid angle of 12.8 msr, and reaction products were momentum analyzed using a magnetic field of approximately 12.9 kG. The focal plane detector Chen et al. (2001) measured the position (momentum), horizontal angle, $`\mathrm{\Delta }E`$ (energy loss in isobutane) and $`E`$ (energy loss in a plastic scintillator). Particle identification (PID) was accomplished mainly by a gate on energy vs. momentum. Gates on other 2-d event plots further refined the PID. Tritons were easily distinguished from the only other incident particles, deuterons. The horizontal tracking information was used to correct the dominant $`(x|\mathrm{\Theta }^3)`$ aberration in the spectrograph’s focussing.

Protons and alpha particles from the decay of unbound states were detected by the Yale Lamp Shade Array (YLSA), an array constructed from five sectors of silicon strip detectors of the type used in the Louvain-Edinburgh Detector Array Davinson et al. (2000). The segments were arranged in an axially-symmetric 5-sided “lamp shade” configuration covering angles from $`\theta =130^{}`$ to $`\theta =165^{}`$, with approximately 14% lab efficiency. The detectors were approximately 500 $`\mu `$m thick; a bias voltage of 100 V was applied to them. Since it was essential that the beam go through the center of the chamber, two collimators were used for tuning: one at the entrance and one at the exit of the scattering chamber; both were 2 mm in diameter. YLSA was protected from particles scattered off the up- and downstream collimators by the detector mount and an aluminum sheet, respectively Bradfield-Smith et al. (2002). After tuning, these collimators were rotated out of the beam path.

In order to ensure proper functioning of the silicon detectors in the environment of electrons ejected from the target by the relatively intense beam, the sides of the detectors closest to the strips’ pn-junctions, which have exposed SiO2 between them, faced away from the target. Also, the target had strong rare-earth magnets immediately to either side to deflect electrons, and the target ladder, aluminum sheet and detector mount were biased with +500 V to attract electrons away from the detector Bradfield-Smith et al. (2002).

A fast trigger was produced by the scintillator in the spectrograph’s focal plane detector. Uninteresting events were discarded using rough PID based on the $`\mathrm{\Delta }E`$ and $`E`$ signals. The “fast-clear” (FCLR) feature of the CAEN V785/V775 ADC’s and TDC’s was used for this hardware cut, and did not reliably clear all modules at high trigger rates. This was undiscovered until after the measurement was made, resulting in the need for efficiency corrections. Fast discriminators on each of the 80 ($`5\times 16`$) silicon detector channels produced logic pulses which were passed through fixed 300-nsec delays to produce time spectra. Decay protons were distinguished from decay alpha particles by their differing energies.

Once PID was used to select the (3He,t) channel and events with possible decay energies were extracted from the data, the momentum spectrum of coincident decays (see Fig. 1) was determined by requiring the time to lie within the peak in the time spectrum. A background spectrum was produced using candidate events with measured times outside the peak. This background spectrum was used to quantify the rate of random coincidences.

The energy calibration for YLSA was determined using a 228Th alpha-particle source. The exact location of the discriminator thresholds varied by strip and was uncertain. The detectors in YLSA have an estimated 2-$`\mu `$m dead layer; in simulations, deposited energies for alpha particles incident on YLSA ranged from 300 to 600 keV, depending on lab angle and the angle of incidence. Because of this, the lowest $`{}_{}{}^{15}O+\alpha `$ resonance for which a reliable branching ratio could be extracted was at $`E_R`$=1020 keV ($`E_x`$=4549 keV). Likewise, the only 18F+p resonance with a reliable proton branching ratio was $`E_R`$=665 keV ($`E_x`$=7076 keV). Above $`E_x`$=7076 keV, the triton energy calibration was uncertain, and states were not well resolved from one another.

## III Simultaneous Fit of the Data Sets

Since all spherical harmonics with non-zero magnetic substates are suppressed at $`\theta =0^{}`$, the tritons were detected at $`\theta =0^{}`$ in order to simplify the calculation of the angular correlations. In order to have sufficient counts for fitting, the detected decay particles were divided into four angles (4 strips $`\times `$ 5 detectors per angle). As in Magnus et al. (1990), the angular correlation function used for breakup to 15Og.s. was

$$W\left(\theta \right)=\underset{M,m}{}|\frac{1}{2}m;lMm|JMY_l^{Mm}(\theta ,\varphi )|^2,$$
(3)
where $`M`$ are the substates of the 19Ne state with spin $`J`$, $`p_M=p_M`$ are the $`J+\frac{1}{2}`$ probabilities of substate populations, and $`m=\pm \frac{1}{2}`$ are the possible magnetic substates of the 15O Satchler (1980). For $`|M|>\frac{3}{2}`$, $`p_M0`$ because the tritons were detected at $`\theta =0^{}`$.

Four separate sets of spectra were analyzed and fit to produce angular distributions. Each set corresponded to a unique set of experimental conditions. Higher ADC thresholds were used in the first data set (D1). The thresholds were lowered for the second and third data sets (D2,D3), which have slightly different energy calibrations from one another. For the fourth data set (D4), the YLSA electronics channels were given a shorter ADC gate than the focal plane channels, in an attempt to reduce pileup in the spectra.

The geometric efficiency of each strip in the array was determined for each resonance using a Monte Carlo simulation of the formation and breakup of the resonances. The $`E_x`$=5351 keV state is well-populated and has an isotropic c.m. angular distribution for its decay products. It is known to decay almost entirely by alpha-particle emission. (Calculations give $`\mathrm{\Gamma }_\gamma /\mathrm{\Gamma }<10^3`$ Langanke et al. (1986); $`\omega \gamma =1.69\pm 0.14`$ eV and $`\mathrm{\Gamma }=1.3\pm 0.5`$ keV have been measured for the analog in 19F which is 499 keV closer to its threshold Wilmes et al. (2002).) When measuring the branching ratio for this state, it was found that the result was consistently less than one. The specific result was dependent on the details of the electronics gating. The three different electronics configurations resulted in three different “coincidence efficiencies”, ranging from 0.581$`\pm `$0.018 (D4) in the worst case to 0.867$`\pm `$0.027 (D1). For each of the runs, the decay data were divided by this factor.

For Ex=4379 keV, the lab energies of the alpha-particles extended below the thresholds for the detector. The exact energy of these thresholds varies with detector strip and is not precisely known. If one (falsely) assumes that all decay alpha particles were detected for this state, a 90% confidence interval for $`\mathrm{\Gamma }_\alpha /\mathrm{\Gamma }`$ of (0.0027,0.0143) is found. This is interpreted as meaning $`\mathrm{\Gamma }_\alpha /\mathrm{\Gamma }>`$0.0027 with a probability greater than 90%. Lower limits of 0.003 and 0.007 for $`\mathrm{\Gamma }_p/\mathrm{\Gamma }`$ of the Ex=6742 and 6861 keV states, respectively, were determined in the same way. In these latter two cases, only data set D4 supports a non-zero result, which is cause for skepticism, especially since calculations based on the measured mirror state in 19F Utku et al. (1998) and a direct measurement of $`\mathrm{\Gamma }_p`$=2.22$`\pm `$0.69 eV Bardayan et al. (2002) estimate $`\mathrm{\Gamma }_p/\mathrm{\Gamma }`$=8.2$`\times 10^4`$ for the Ex=6742 keV state Utku et al. (1998). The ratio extracted from D4 is significantly different only for these two states.

As can be seen in Fig. 1, the energy (Ex) resolution of the triton spectrum was about 40 keV, so the peaks for the states at Ex=4600 and 4635 keV were not well resolved. Since the Ex=4635 keV state is not expected to emit significant numbers of alpha-particles (due to its $`l=7`$ angular momentum barrier), in our analysis of the coincidence data for these two states we attributed all of the measured decay alpha-particles to the Ex=4600 keV state.

For the 330-keV proton resonance at $`E_x`$=6742 keV, we measured the spin and parity using the angular distribution of its decay alpha-particles. The 20Ne(3He,$`\alpha `$)19Ne angular distribution measured for this state shows $`l=1`$ character, implying negative parity and $`J=\frac{1}{2}`$ or $`\frac{3}{2}`$ Garrett et al. (1970). $`J^\pi =\frac{1}{2}^{}`$ would imply an isotropic correlation function for the decay alpha-particles, but Fig. 2a shows that $`J^\pi =\frac{3}{2}^{}`$ results in a much better fit to the data.

A previous high-resolution measurement has shown that the (3He,t) reaction populates the 450 keV proton resonance Utku et al. (1998), and not an $`l=1`$ resonance (isospin mirror of the Jπ=$`\frac{3}{2}^{}`$, Ex=6891 keV state in 19F) predicted to lie near 430 keV Shu et al. (2003). The decay alpha-particle angular distribution in Fig. 2b shows that a spin assignment of $`\frac{7}{2}^{}`$ is the appropriate choice for the peak at this location in our spectrum.

The measurement of $`\mathrm{\Gamma }_p/\mathrm{\Gamma }=0.387\pm 0.016`$ for the 665-keV proton resonance is in agreement with the previously measured value Utku et al. (1998). It disagrees with the branching ratio implied by the direct measurement performed at Louvain-la-Neuve Graulich et al. (2000), but agrees with the measurements made at Oak Ridge National Laboratory Bardayan et al. (2001).

## IV Conclusions

The branching ratio measurements presented here are compared with previous measurements in Tables 1 and 2. In cases where other measured values exist, weighted averages are given. The measurement provides support for the $`J^\pi `$=$`\frac{3}{2}^{}`$ assignment for the 330-keV 18F(p,$`\gamma `$) resonance, which dominates the reaction rate above 450 MK Utku et al. (1998). The 450-keV 18F+p resonance is also confirmed to have $`J^\pi `$=$`\frac{7}{2}^{}`$, re-affirming its weak role in astrophysical processes Utku et al. (1998). The measured branching ratio for the 665-keV 18F+p resonance, which strongly affects both 18F(p,$`\gamma `$) and 18F(p,$`\alpha `$) astrophysical rates Utku et al. (1998), differs from one recent direct measurement Graulich et al. (2000) and agrees with another Bardayan et al. (2001). Therefore, these reaction rates can now be calculated with greater confidence.

## V Acknowledgements

This work was supported under US Department of Energy Grant No.DE-FG02-91ER-40609.


text:
# Efficiency of stripping mechanisms

## 1. Introduction: mechanisms to remove matter from galaxies

Among the various possibilities to explain the stripping of matter (gas and stars) from galaxies in clusters, the principal actors can be classified in three groups:
1. Tidal forces:
– interaction with a companion, merger: in this case, a correlation between morphological type (T) and density ($`\mathrm{\Sigma }`$) should be expected, (T-$`\mathrm{\Sigma }`$ relation)
– interaction with the cluster; then a correlation between type and radius in the cluster is expected (T-R relation)
– harassment due to numerous interactions at high velocity and density
2. ICM-ISM interactions:
ram pressure stripping, but also thermal evaporation, turbulent, viscous stripping; these are purely hydrodynamical mechanisms, and should affect only the diffuse gas. However, they are acting simultaneously with the others, and relative roles are hard to disentangle. Since they are efficient only when the cluster is formed, and the ICM gathered, have they enough time to act? Or have tides acted before?
3. Outflows due to violent events:
– starbursts and winds
– AGN jets and outflows

All these processes result in morphological type changes for galaxies, and stripping of their gas, therefore star formation quenching, or ”starvation” as is observed in clusters. The delicate issue is that many mechanisms are able alone to account for the stripping/quenching, and very specific tests have to be found to disentangle what is happening.

## 2. Observational clues

### 2.1. Intra-cluster diffuse light

One of the clear evidence of tidal interactions and stripping is the existence of intra-cluster diffuse light (ICL): these intergalactic stars, stripped from their parent galaxies by tidal interactions, represent a large fraction of the total stellar mass of the cluster, between 10-40% (cf Figure 1, Feldmeier et al 2003). Cluster images at low luminosity levels show evidence of tidal debris in the form of plumes and arclike structures (example of the Centaurus cluster, Calcáneo-Roldán et al 2000). The quantity of ICL does not appear to depend on cluster radius, but more on the surface density of galaxies ($`\mathrm{\Sigma }`$), which favors the interactions between galaxies.

Although CCD images are now able to reveal ICL in most clusters, a large sensitivity for this diffuse component is gained from Planetary Nebulae tracers, without the problems of flat fielding, etc, since they are detected by emission lines (Feldmeier et al 1998, Arnaboldi et al 2002) The intra-cluster stars have moderate metallicity (Durrell et al 2002), which supports the scenario of their stripping from intermediate mass galaxies. These tidal debris and plumes are expected from simulations of galaxy clusters (cf Dubinski 1998), even more prominent than what is observed. However, the background noise dilutes the weaker features, explaining the difficulty to observe them clearly (e.g. Mihos 2003).

### 2.2. Larger fraction of blue galaxies at z=0.4 (BO effect)

It has been known for a long time that there exists in clusters a larger fraction of blue galaxies at increasing redshift (Butcher & Oemler 1978, 1984). These blue galaxies indicate much more star formation in the recent past, and correspond to irregular shapes in the clusters. The existence in z=0.4 clusters of sign of tidal interaction/mergers also confirm that clusters have evolved very recently: in the last few Gyrs, there was a much larger fraction of perturbed galaxies, late-types and starbursts, as if the cluster had relaxed only since then. Rings of star formation were much more frequent than 2-arms spirals, contrary to what is found today (Oemler et al 1997). These rings could be due to bars triggered in tidal interactions. Part of them could also be due to fast encounters, expected in galaxy clusters, that lead to head-on collisions like the Cartwheel. Alternatively galaxies, through harassment, could be stripped at this epoch of their dark halos, de-stabilising disks. and triggering more violent star formation. These tidal interactions visible at z=0.4, must have profoundly and rapidly modified the galaxy morphologies, since at z=0.2, the evolution is almost terminated. Milder effects are observed by Balogh et al (1999) in an X-ray selected sample of clusters (CNOC1), who suggest a more gradual decline of star formation.

### 2.3. Star formation rate versus environment and density

In an H$`\alpha `$ line study of 11000 galaxies in the 2dF survey, over 17 galaxy clusters, Lewis et al. (2002) find the star formation rate (SFR) increasing gradually from low values at the cluster centers, towards the field value at about 3 virial radii. They find a strong correlation between SFR and local projected density, as soon as the density is above 1 galaxy/Mpc2, independent of the size of the structure (i.e. also valid in groups). Gómez et al (2003) find also a strong SFR-$`\mathrm{\Sigma }`$ relation with the early data release of the SDSS, the SF-quenching effect being even more noticeable for strongly star-forming galaxies. The same break of the SFR-$`\mathrm{\Sigma }`$ relation is observed at 1 galaxy/Mpc2. This relation is somewhat linked to the morphological type-density (T-$`\mathrm{\Sigma }`$) relation, but cannot be reduced to it, since at any given type, the SFR-$`\mathrm{\Sigma }`$ relation is still observed. This strong relation valid even outside cluster cores is a precious clue to derive the dominant mechanisms.

### 2.4. Morphological segregation

From the morphological segregation in nearby clusters, drawn by Dressler et al (1980), it is now possible to see the evolution from about 5 Gyrs ago, at z=0.4 (Dressler et al 1997, Figure 2): at z=0, there was the same T-$`\mathrm{\Sigma }`$ correlation for relaxed or non-relaxed clusters, but it is no longer true at z=0.4. As main lines of evolution, there is at z=0.4 the same fraction of ellipticals than at z=0, but a much smaller fraction of S0s; at z=0.5, the fraction of lenticulars is 3 times lower than now. This suggests that ellipticals form early, before the cluster virialisation. In the hierarchical scenario, clusters form out of loose groups mergers, and it is likely that ellipticals are the result of mergers in groups, before the formation of the cluster.

S0’s are transformed from spirals in virialised clusters, in a few Gyrs time-scale. The study of stellar populations, and the spectral classes of cluster galaxies at z=0.4 reveals that star formation is quenched with respect to the field (Poggianti et al 1999). At z=0.4, passive and post-starburst (E+A, or k+a) spirals are much more frequent than in the field. It appears that the mechanism reponsible for that must act on shorter time-scales than the mechanism reponsible for the transformation into S0s.

### 2.5. Density correlation

In cluster regions where the density is not centrally symmetric, it is possible to compare the morphology-radius (T-R) and morphology-density (T-$`\mathrm{\Sigma }`$) relations. The latter (T-$`\mathrm{\Sigma }`$) appears always better than the T-R relation (Treu et al. 2003, example of Cl0024+16 at z=0.4). Galaxies are more aware of their local density than cluster location.

The morphological segregation as a function of radius is quite significant for radii lower than 200kpc (cf Figure 3). The fraction of early-type galaxies drop steadily until 1 Mpc, or nearly the virial radius. The correlation with radius is then weak, while several over-dense regions have galaxies with morphology typical of their high density. It seems that gas depletion and morphological transformation are already well advanced in groups, before forming the cluster. The mild gradient in the morphological mix outside the virial radius could be due to harassment and starvation (ICM interactions are not operating there). It is only upon arrival in the central regions (R $`<`$ 200 kpc) that substructures are erased, as indicated by the tight correlation between cluster radius and $`\mathrm{\Sigma }`$.

### 2.6. HI deficiency

Solanes et al. (2001) have recently made a data compilation on 1900 galaxies, and conclude that about 2/3 of galaxy clusters are HI deficient in their centers. Many observations demonstrate that the interaction with the hot gas (ICM-ISM interactions) might be responsible for this HI stripping in cluster galaxies: large deficiencies (up to a factor 100), deficiency as a function of distance from the central X-ray peak, radial orbits of the stripped galaxies, that allow them to explore the dense hot center. In her review, van Gorkom (2003) describes convincing individual cases proving ram pressure stripping, like galaxies in the Virgo cluster with perturbed and reduced-size HI disks while the stellar disk is normal. Ram pressure stripping appears quite efficient and rapid, playing a role in only one cluster crossing-time.

However, the correlation between HI deficiency and X-ray properties of the cluster is not observed (LX, TX), which is surprising (Solanes et al 2001). Instead, there are correlations with galaxy properties, early-type and probably dwarf spirals are more easily stripped than the intermediate spiral types. That early-types are more stripped is not only due to their position in the cluster, but their deficiency is larger than from late spirals at each cluster radius, until 4 Mpc from the center.

## 3. Simulation clues

### 3.1. Tidal interactions and harassment

Tidal interactions in galaxies were thought marginal because of high velocities, and un-resonant interactions. However, the large number of interations can accumulate perturbations and truncation effects: this has been called “harassment” by Moore et al. (1996), i.e. frequent high-velocities close encounters.

Gnedin (2003) has recently simulated the formation of galaxy clusters in the frame of a hierarchical cosmological scenario, varying the cosmological parameters. Tidal interactions determine the galaxy evolution, and are intensified by the density irregularities, either the presence of massive galaxies, or the infalling groups of galaxies, still not relaxed. These substructures favor the interaction, the typical frequency being estimated at 10 interactions at 10kpc impact parameter per galaxy. Mergers occur essentially at the cluster formation, and are very rare today. This means that elliptical galaxies predate the cluster, as already found by Merritt (1984). Later the tidal interactions can transform spirals to lenticulars, and explain their large fraction increase in the last Gyrs. Tidal interactions truncate massive dark matter halos and thicken stellar disks, increasing disk stability and quenching star formation. Dwarf galaxies can be totally disrupted. The collision rate per galaxy strongly decreases with time, from 8 to 2 per Gyr along the cluster life-time.

### 3.2. Cluster tidal field

The tidal field of the cluster itself has strong dynamical influence on galaxies, in particular in their extended halos. First the dark matter haloes are stripped, and form a common halo, but also the gas reservoirs that replenish the interstellar medium of field galaxies all along their lives, by accretion along gas filaments, is stripped also, and this could easily explain the starvation, and gradual decline in star formation of cluster galaxies. This phenomenon was first invoked by Larson et al (1980), and simulated by Bekki et al (2001). The latter assume an accretion rate of 1 M/yr for a normal field galaxy, and show that the tidal field of the cluster efficiently removes the gas reservoir from a galaxy, and consequently its fueling of star formation. This tidal truncation does not depend very strongly on the orbit of the galaxy in the cluster, and the resulting SF-quenching is widespread through the cluster (contrary to what is expected from ICM interactions). Once their gas reservoir is stripped, spiral galaxies will slowly be transformed into lenticulars by harrasment, thickening and shortening their stellar disk.

### 3.3. Ram pressure and ICM-ISM interactions

A large variety of models have been simulated, since the phenomena associated with ICM interactions depend on many physical assumptions about the small-scale structure of the gas, instabilities, equivalent viscosity, temperatures etc.. With an isothermal SPH gas model, Abadi et al (1999) show that HI gas is effectively stripped in the core of rich clusters, for disks oriented perpendicular to the wind. Galaxies can lose 80% of their gas, the final disk being restricted to 4kpc radius, in a time-scale of 107 yrs. However, in the outer parts of the cluster, or for inclined disks with respect to the wind, the process is much less efficient.

With a finite-difference code, Quilis et al. (2000) show that viscous coupling could favour the stripping; they show that a hole in the center of the galaxy (mimicking the frequent HI depletion in the central regions), could fragilize the gas disk, and enhance considerably the stripping efficiency. Holes can have an influence at several scales. If star formation has already formed shells through supernovae and winds, ram pressure can enlarge the holes in the disk (Bureau & Carignan 2002). Supernovae and winds alone are not efficient enough, except may be in small dwarfs (Dekel & Silk 1986, Martin 1998).

Vollmer et al (2001) use sticky particles for the gas, and follow galaxies on their orbit through the cluster: they show that the ICM has only a significant action in the cluster core, and the stripped gas then falls back on the galaxy, once its orbit gets out of the core. Schulz & Struck (2001) show that the stripping is a multi-step process: the outer gas is quickly stripped, while the inner gas is compressed, and forms a ring. The compressed gas could give rise to triggered star formation in a small starburst.

It is also possible that the gas reservoir required to replenish the ISM of galaxies is hot and diffuse, as assumed by Bekki et al (2002), who show through ram pressure simulations, that the halo will be efficiently stripped, if its density is typically lower than 3 10-5 cm-3. In that case, the global tidal stripping from the cluster and the ram pressure compete to strip the gas reservoirs, and contribute to form passive or anaemic spirals, that will slowly be transformed into S0s.

## 4. Outflows: stellar and nuclear activity

Enhanced stellar activity, when spiral galaxies infall into the cluster, is observed (Kenney & Yale 2002). If the galaxies have been stripped by the cluster global tide of their halo, this can favor the escape of the winds. In a wide sample of galaxies, Kauffmann et al (2003) have shown that star formation efficiency is a strong function of surface density. Low surface density dwarf objects (LSB) are unevolved objects, in which stellar feedback have prevented rapid star formation, by ejection of their gas. The energy of supernovae is enough to disperse the gas, when the mass of the galaxy has fallen below a threshold of 3 1010 M, the observed transition between low and high surface density galaxies (Dekel & Woo 2003).

Nuclear activity (AGN) could also be invoked to provoke gas outflows, and remove gas from galaxies. AGN have been found less frequent in cluster environment, by a factor 5, with a frequency of only 1 percent, compared to 5 percent in the field (Dressler et al. 1985). Recent observations show however that they might be more frequent in X-rays, suggesting an obscuration effect (Martini et al 2002). The mechanical and heating energy of AGN has a strong feedback effect to reverse and self-regulate the gas cooling in the centers of elliptical galaxies, groups and clusters (e.g. Ciotti & Ostriker 2001). AGN feedback has been invoked to account for recent X-ray observations incompatible with the old quasi-state cooling flow model: absence of extremely cool gas ($`<`$ 1kev) in the center of cooling flow clusters, presence of cold bubbles related to AGN radio lobes, etc… The amount of effectively cooling gas has been revised downwards. The cooling flow could be intermittent, with alternate periods of AGN activity, inflow and outflow. A cold gas phase is observed, as shown in the Abell 1795 cooling flow (figure 4). This cluster center is not yet relaxed, with a cooling wake triggered by the central cD motion (Fabian et al. 2001). The X-ray data from the inner 200kpc indicate a mass deposition rate of about 100 M/yr (Ettori et al. 2002).

## 5. Discussion

### 5.1. Chronology of events

To find the different roles of the various stripping mechanisms discussed above, it is essential to go back to the formation of the cluster, and the chronology of events. The dense ICM required for efficient ram pressure is in place only after the formation of the cluster. Also the phenomenon is proportional to the square of velocity dispersion, so virialisation should have occured, for all galaxies to acquire their high velocities. In the hierarchical scenario of structure formation, groups form first, with low velocity dispersion, and a large frequency of resonant tidal interactions and frequent mergers. Giant ellipticals are formed at this stage by spiral galaxy mergers in groups. This is part of the morphological segregation, since the fraction of spirals decreases.

Progressively, the gas in galaxy haloes is stripped by tidal interactions and heated by shocks to the virial temperature of the growing structure, and the importance of the ICM interactions will grow. The ICM itself is formed through gravitational forces. A large fraction of the ICM gas comes from already processed galaxy disks, since its metallicity is important (1/3 solar).

At the present epoch, when clusters are virialised and still relaxing, the merger rate has fallen to zero (at least major mergers), and ellipticals are only passively evolving. Ram pressure and the global cluster tide are almost equally efficient to strip gas from infalling galaxies (i.e. Bekki et al. 2002).

The fact that the T-$`\mathrm{\Sigma }`$ relation is even tighter than the T-R relation, combined to the extension of the T-$`\mathrm{\Sigma }`$ relation from clusters down to groups (Ramella et al. 1999, Helsdon & Ponman, 2003) suggest that two-body gravitational interactions are dominating. The transformation of spirals into S0s in the outer parts of clusters, depending on the local density, does not pleade in favor of the ICM interactions as the main mechanism of the morphological segregation.

Low density galaxies (dwarfs and LSB) entering a cluster can be entirely disrupted by tidal shocks, they form streams of stars contributing to the intra-cluster light. Tidal phenomena lead similarly to a morphology-density relation in loose groups, like that observed in the local group (Mayer et al 2001).

### 5.2. Metallicity considerations

In rich clusters, the mass of the hot gas can be much larger than the mass of baryons in galaxies. Most of the baryons are in the hot ICM, since it represents nearly the baryon fraction of the matter in the universe $`f_b0.16`$. Given its relatively high metallicity, it is then not surprising that most of the metals in a cluster comes from the ICM: there is about twice more Fe in the ICM than in galaxies (Renzini 1997, 2003).

Since the metals are synthetized in galaxies, this means that either part of the hot gas comes directly from galaxies (by stripping, or disruption), or that stellar winds and supernovae have enriched the ICM. In fact, both sources of metals should be there, since metals expelled by SNe would not be sufficient.

Figure 5 shows that the iron mass-to-light ratio is about constant as a function of the virial temperature of the structure (below 2 kev, the estimations are less simple to derive, and there could be biases). These observations suggest that metals are ejected via winds (SNe or AGN), not ram pressure, since there is no dependence on richness, or cluster velocity dispersion, but only the dispersion of individual galaxies (Renzini 2003). There is the same MFe/LB in clusters and galaxies, implying the same processing in clusters than in the field.

Ellipticals in the field or in clusters have the same properties, confirming that these galaxies have been formed before the clusters. Since the stellar mass is essentially in elliptical galaxies today (3/4 of the mass in spheroids, 1/4 in disks, less than 1% in Irr, Fukugita et al. 1998, and Es are dominating even more in clusters), most cluster stars are formed before cluster formation.

Stars and corresponding metals have likely been made at z$``$ 2-3, at the peak of stellar activity in the universe. Part of the iron has been made through SNII quickly, then by SNIa, 1Gyr after the main episodes of star formation. Clusters have not lost iron, nor accreted pristine material since the ratio between $`\alpha `$ elements (made in SNII) and iron is about solar in the ISM of all clusters (with no or little variation from cluster to cluster). This means that the metals come from normal stellar nucleosynthesis, with similar ratio of type Ia to type II SNe, as well as the same global IMF, etc..

Also, at the present time SNIa continue to enrich in Fe the medium essentially near the cD at the center; it has been shown that clusters have almost negligible metallicity gradients, except those with cooling flows, and bright central galaxies (Böhringer, this symposium).

## 6. Conclusion

The various clues brought by observations and simulations allow to draw the following conclusions, about the stripping mechanisms and the recycling of matter in clusters:
– Tidal forces are at play from the beginning of structure formation, and are the dominant factor in the T-$`\mathrm{\Sigma }`$ relation; they need long time-scales (1Gyr) in groups that will finally (recently, z$``$1) merge in a cluster
– ICM interactions have entered in action only recently (after the virialisation of rich clusters); they are very efficient in HI stripping from galaxy disks infalling today into the cluster in nearly radial orbits. The corresponding stripping time-scales are very short (107-108 yr).
– Starburst and AGN winds have played a role continuously (with a peak at z=2-3), and especially in the metallicity distribution
– E’s have formed in groups before the cluster, S0’s have been transformed from spirals in rich clusters in the last 5 Gyr; their gas reservoir has been stripped through both the global tide of the cluster and the ICM interaction, and their star formation quenched. Harassment progressively truncates and thickens their disks.

The difficulty to derive exactly the cluster evolution is that the mechanisms may act simultaneously, and also they show some duality: for instance, on one hand tidal interactions trigger some star formation, leading to the BO-effect, and observed larger fraction of blue galaxies at z=0.4 in clusters, but also the tidal truncation of gas reservoirs implies a gradual decline in star formation in the last Gyrs.

## References

Abadi, M.G., Moore, B., Bower, R.G.: 1999, MNRAS 308, 947

Arnaboldi, M., Aguerri, J., Napolitano, N. et al.: 2002, AJ 123, 760

Balogh, M., Morris, S. L., Yee, H. K. C. et al.; 1999, ApJ 527, 54

Balogh, M., Navarro, J., Morris, S. L. et al.: 2000, ApJ 540, 113

Bekki, K., Couch, W.J., Shioya, Y.: 2001 PASJ 53, 395

Bekki, K., Couch, W.J., Shioya, Y.: 2002 ApJ 577, 651

Bureau, M., Carignan, C.: 2002 AJ 123, 1316

Butcher, H., Oemler, A.: 1978, ApJ 219, 18

Butcher, H., Oemler, A.: 1984, ApJ 285, 426

Calcáneo-Roldán, C., Moore, B., Bland-Hawthorn, J., et al.: 2000 MNRAS 314, 324

Ciotti, L., Ostriker, J.: 2001 ApJ, 551, 131

Dekel, A., Silk, J.: 1986 ApJ 303, 39

Dekel, A., Woo, J.: 2003, in press (astro-ph/0210454)

Dressler, A., Thompson I., Shectman S.: 1985 ApJ 288, 481

Dressler, A., Oemler, A., Couch, W.J., et al. : 1997 ApJ 490, 577

Dressler, A., Smail, I., Poggianti, B. et al.: 1999 ApJS 122, 51

Dubinski, J.: 1998, ApJ 502, 141

Durrell, P., Ciardullo, R., Feldmeier, J. et al.; 2002 ApJ 570, 119

Ettori, S., Fabian, A. C., Allen, S. W., Johnstone, R. M. : 2002 MNRAS 331, 635

Fabian, A. C., Sanders, J. S., Ettori, S. et al.: 2001, MNRAS 321, L33

Feldmeier, J., Ciardullo, R., Jacoby, G.; 1998, ApJ 503, 109

Feldmeier, J., Mihos, J.C., Morrison, H., et al.: 2003 in Clusters of Galaxies: Probes of Cosmological Structure and Galaxy Evolution, Carnegie Observatories Symposium III. (astro-ph/0303340)

Fukugita, M., Hogan, C.J., Peebles, P.J.E.: 1998, ApJ 503, 518

Gnedin, O.: 2003, ApJ 582, 141 and ApJ 589, 752

Gómez, P., Nichol, R., Miller, C. et al.: 2003, ApJ 584, 210

Helsdon, S. F., Ponman, T. J.: 2003 MNRAS 339, L29

Kauffmann, G., Heckman, T., White, S.D.M. et al: 2003, MNRAS 341, 54

Kenney, J.D.P., Yale, E.E.: 2002 ApJ 567, 865

Larson, R. B., Tinsley, B. M., Caldwell, C. N.: 1980 ApJ 237, 692

Lewis, I., Balogh, M., de Propris, R. et al. 2002: MNRAS 334, 673

Martin, C.: 1998 ApJ 506, 222

Martini, P., Kelson, D. D., Mulchaey, J. S., Trager, S. C.: 2002 ApJ 576, L109

Mayer, L., Governato, F., Colpi, M. et al. 2001 ApJ 547, L123

Merritt, D.: 1984 ApJ 276, 26

Mihos, J.C.: 2003 in Clusters of Galaxies: Probes of Cosmological Structure and Galaxy Evolution, Carnegie Observatories Symposium III. (astro-ph/0305512)

Moore, B., Katz, N., Lake, G. et al.: 1996, Nature 379, 613

Moore, B.: 2003, in Clusters of Galaxies: Probes of Cosmological Structure and Galaxy Evolution, Carnegie Observatories Symposium III. (astro-ph/0306596)

Oemler, A., Dressler, A., Butcher, H.: 1997 ApJ 474, 561

Poggianti, B., Smail, I., Dressler, A. et al. : 1999 ApJ 518, 576

Quilis, V., Moore, B., Bower, R. 2000 Sci. 288, 1617

Ramella, M., Zamorani, G., Zucca, E. et al.: 1999 A&A 342, 1

Renzini, A.: 1997, ApJ 488, 35

Renzini, A.: 2003, in Clusters of Galaxies: Probes of Cosmological Structure and Galaxy Evolution, Carnegie Observatories Symposium III. (astro-ph/0307146)

Salomé, P., Combes, F.: 2003, in prep

Schulz, S., Struck, C. : 2001 MNRAS 328, 185

Solanes, J., Manrique, A., García-Gómez, C. et al.: 2001, ApJ 548, 97

Treu, T., Ellis, R. S., Kneib, J.-P., et al.: 2003, ApJ 591, 53

van Gorkom, J.H.: 2003, in Clusters of Galaxies: Probes of Cosmological Structure and Galaxy Evolution, Carnegie Observatories Symposium III. (astro-ph/0308209)

Vollmer, B., Cayatte, V., Balkowski, C., Duschl, W. J. : 2001, ApJ 561, 708