text:
# RELATIVISTIC JETS IN COLLAPSARS

## 1 INTRODUCTION

Growing evidence connects GRBs to the death of massive stars. Analysis by Frail et al. (2001) and others of the radio afterglows of “long-soft” GRBs suggests that, despite great diversity in apparent brightnesses, these events have a total kinetic energy in relativistic matter tightly clustered around $`3\times 10^{51}\mathrm{erg}`$ (Frail et al., 2001) \- a supernova-like energy. “Bumps” resembling the light curves of Type I supernovae have also been seen in the optical afterglows of at least four GRBs (GRB 980326, Bloom et al. 1999; GRB 970228, Reichart 1999, Galama et al. 2000; GRB 011121, Bloom et al. 2002, Garnavich et al. 2002; and GRB 020405, Price et al. 2002) and GRB 980425 has been associated with an optical supernova, SN 1998bw (e.g., Galama et al. 1998; Iwamoto et al. 1998; Woosley, Eastman, & Schmidt 1999). The observational evidence that (long-soft) GRBs are associated with regions of star formation has become overwhelming (Bloom, Kulkarni, & Djorgovski, 2002). Given the evidence for beaming and relativistic motion, it is probable that at least one major subclass of GRBs is a consequence of massive stars that, in their explosive deaths, produce relativistic jets.

Here we examine the passage of relativistic jets through a collapsing massive star and their breakout. We begin with a rotating star that has already experienced 10 seconds of collapse. Ten seconds is the nominal time for a disk to form around a black hole and the polar region to become sufficiently evacuated for a polar jet to propagate (MacFadyen & Woosley, 1999). The initial progenitor is a helium core of $`15M_{\mathrm{}}`$ evolved by Heger & Woosley (2002) to iron core collapse with approximations to all (non-magnetic) forms of angular momentum transport included. The missing 10 seconds is followed using a non-relativistic two-dimensional code as in MacFadyen & Woosley (1999). We pick up the calculation when the jet, which presumably began in a region $`30\mathrm{km}`$ in size, has already reached a radius of $`2000\mathrm{km}`$ and do not consider what has gone on inside. For present purposes, details of whether the jet was made by black hole angular momentum, MHD processes in the disk, or neutrinos do not concern us. The jet is initiated in a parametric way based upon its power, opening angle, Lorentz factor, and internal energy. Its propagation to the stellar surface at $`800,000\mathrm{km}`$, and its interaction with the stellar mantle is then followed (§ 3.1). Additional calculations (§ 3.3) also examine what happens to the jet immediately after it escapes the star and converts its residual internal energy into additional relativistic motion.

We find, in agreement with Aloy et al. (2000), that the passage of the jet through the star leads to its additional collimation. We also find that instabilities along the beam’s surface lead to mixing with the nearly stationary stellar material and cocoon. The mixing produces variations in the mass loading, and therefore the Lorentz factor of the jet. The opening angle of the jet also varies with time, gradually, if irregularly, growing as the star is blown aside. These results, discussed in § 4, have important implications for the observed light curves and energies of GRBs and imply that what is seen may vary greatly with viewing angle. In particular, we predict the existence of a large number of low energy GRBs with mild Lorentz factors (§ 4.3) that may be related to GRB 980425/SN 1998bw and to the recently discovered “hard x-ray flashes” (Heise et al., 2001).

Finally, we consider the breakout of the jet. As the shock breaks out a small amount of material being pushed ahead by the jet head is accelerated to relativistic speeds. Interaction of this material with the stellar wind of the progenitor will produce a transient of some sort (Woosley, Eastman, & Schmidt, 1999; Tan, Matzner & McKee, 2001). We speculate that this is the origin of a hard precursor to GRBs that, at least in some cases, might be characterized as a “short-hard GRB” in isolation(§ 4.4).

## 2 COMPUTATIONAL PROCEDURE AND ASSUMPTIONS

### 2.1 The Relativistic Code Employed

In order to study relativistic jets in collapsars, we developed a special relativistic, multiple-dimensional hydrodynamics code similar to the GENESIS code (Aloy et al., 1999). Our program employs an explicit Eulerian Godunov-type method with numerical flux calculated using an approximate Riemann solver: Marquina’s flux formula. The PPM algorithm is used to reconstruct the variables so that high-order spatial accuracy is achieved. Second- or third-order Runga Kutta methods are employed for the time integration. This method has previously been used successfully to study relativistic jets in collapsars and many other ultra-relativistic flows (Aloy et al., 1999). The code has been extensively tested and the results are comparable with those found in the literature (Martí & Müller, 1999, and references therein). To model collapsars more accurately, the code was revised to include a realistic equation of state (Blinnikov, Dunina-Barkovskaya, & Nadyozhin, 1996), that, in particular, allows a Fermi gas of electrons and positrons with arbitrary relativity and degeneracy. Radiation and an ideal gas of nuclei are also included. The code was also modified to include the advection of nuclear species. Nuclear burning, which is not important in the present problem, is not considered. Because we are interested in phenomena far from the central black hole, gravity is treated in the Newtonian approximation. We do not include the self gravitational force of the mass on our grid, but only a radial force from the mass inside our inner boundary condition.

### 2.2 The Stellar Model

Our initial model is derived from a 15 $`M_{\mathrm{}}`$ helium star of 0.1 solar metallicity with an initial surface rotation rate of 10% Keplerian on the helium burning main sequence. This Wolf-Rayet star might be the result of a star with initially 40 $`M_{\mathrm{}}`$ on the main sequence. The rotating helium star is evolved to iron core collapse (defined as when any part of the star is collapsing at $`10^8\mathrm{cm}\mathrm{s}^1`$) including angular momentum transport, but neglecting mass loss and magnetic fields (Heger, Langer, & Woosley 2001; Heger & Woosley 2002). It is presumed that this helium star has lost its hydrogen envelope either to a companion star or by stellar winds. The collapsing, but still initially spherical star is mapped into a two-dimensional non-relativistic code (MacFadyen & Woosley, 1999) and the inner 2 $`M_{\mathrm{}}`$ of iron core, which is assumed to collapse to a black hole, replaced by an absorbing inner boundary condition at $`10^7\mathrm{cm}`$. Poisson’s equation is then solved using a post-Newtonian gravitational potential. Neutrino cooling by thermal processes and pair capture are included, as well as the photodisintegration of nuclei into nucleons in the inner part of the disk, once one has formed. The calculation is run until a low density channel ($`\rho 10^6\mathrm{g}\mathrm{cm}^3`$) has cleared along the rotational axis. This takes approximately $`10\mathrm{s}`$ after the initial core collapse. The collapsing rotating star is then mapped into our two-dimensional special relativistic code(§ 2.1), the inner boundary moved out to $`2\times 10^8\mathrm{cm}`$, and a $`5.0M_{\mathrm{}}`$ point mass placed at its center.

During this $`10\mathrm{s}`$, rotation and infall appreciably modify the core structure, but only inside of $`5000\mathrm{km}`$. Outside $`10,000\mathrm{km}`$, one might just as well use a non-rotating presupernova model. The present study with an inner boundary at $`2000\mathrm{km}`$ is thus moderately sensitive to the post-collapse evolution of the rotating progenitor. For example, the density in the equatorial plane at $`2000\mathrm{km}`$ and $`10\mathrm{s}`$ is ten times greater than that at the same radius along the axis. However, this non-spherical density structure does not critically influence jet propagation. Our two-dimensional starting model is available to others and might be the basis for comparison calculations that are more sensitive to the structure of the inner star.

### 2.3 Initial Conditions

A two-dimensional spherical grid ($`r,\theta `$) is employed consisting of 480 logarithmically spaced radial zones, 120 uniform angular zones in the polar region ($`0\mathrm{°}\theta 30\mathrm{°}`$) and 80 logarithmically spaced angular zones in the range $`30\mathrm{°}\theta 90\mathrm{°}`$. The initial model is remapped onto this grid with an inner boundary at $`2\times 10^8\mathrm{cm}`$, and an outer boundary at $`9\times 10^{10}\mathrm{cm}`$.

As the initial condition for the calculation, a highly relativistic jet is assumed to have formed interior to the inner boundary of the computational grid. The jet is implemented in the calculation as an inner boundary condition. An axisymmetric outflow is injected in a purely radial direction through the inner boundary at $`2\times 10^8\mathrm{cm}`$ within a half-angle $`\theta \theta _0`$, where $`\theta _0`$ is a free parameter. The processes producing the jet are uncertain, but are presumed to have endowed it with a large energy per baryon. At the origin this may have been in the form of a large internal energy (as in the neutrino-version of the collapsar model or a jet energized by magnetic reconnection) or the jet may have been born with low entropy and highly directed motion (as in some versions of MHD accelerated jets). As we shall see, processes deep inside the star tend to erase the distinction.

The jet is specified by its power, $`\dot{E}`$, initial opening half-angle, $`\theta _0`$, initial Lorentz factor, $`\mathrm{\Gamma }_0`$, and the ratio of its kinetic energy111The kinetic energy density in the lab frame is defined as $`\rho \mathrm{\Gamma }(\mathrm{\Gamma }1)`$, here $`\rho `$ is the local rest mass density, and $`\mathrm{\Gamma }`$ is the Lorentz factor in the lab frame to its total energy, $`f_0`$. Of these four parameters, the power is the easiest to address. First, observationally, we know that the total energy in relativistic ejecta is of order a few times $`10^{51}\mathrm{erg}`$ for the entire event (e.g., Frail et al. 2001). The energy in subrelativistic ejecta, i.e., the supernova, is probably comparable or greater (Woosley et al. 1999; Iwamoto et al. 1998). A total energy of order $`10^{52}\mathrm{erg}`$ is thus reasonable. From another perspective, roughly one to several solar masses of material will accrete into the black hole during the principal phase of the GRB. If a few tenths of a percent of the rest mass energy is used to power jet-like outflows, one again obtains $`10^{52}\mathrm{erg}`$. In the lab frame, the duration of a (long-soft) GRB is $`10\mathrm{s}`$, so a jet power of $`10^{51}\mathrm{erg}\mathrm{s}^1`$ is reasonable and probably accurate to an order of magnitude. We employed this energy for each jet so that the total is $`2\times 10^{51}\mathrm{erg}\mathrm{s}^1`$ (or $`2\times 10^{52}\mathrm{erg}`$ in $`10\mathrm{s}`$). In light of recent developments that have tended to downgrade the total energy of a GRB to a few times $`10^{51}\mathrm{erg}`$ (Frail et al., 2001) and even suggest that the kinetic energy after the GRB is over is only a few times $`10^{50}\mathrm{erg}`$ (Panaitescu & Kumar, 2001), this may be too large. So we also ran one set of calculations with an energy 30% as great. This is still large compared with the estimate by Panaitescu & Kumar (2001).

The ratio of kinetic energy to total energy, $`f_0`$, depends on how the jet was born and its history before it entered the computational grid and is more difficult to estimate. We assume, at $`2000\mathrm{km}`$, that our jet has comparable internal and kinetic energies in the laboratory frame. Specifically, we adopt an initial Lorentz factor, $`\mathrm{\Gamma }`$, of 50 for Models JA and JB, and $`f_0=1/3`$. These values are reasonable for a jet which began at say $`30\mathrm{km}`$ consisting initially of a small amount of matter at rest with about $`150`$ times its rest mass in internal energy. That is, in the notation of Mészáros, Laguna, & Rees (1993), $`\eta =150`$. For a jet of constant opening angle, $`\theta _0`$, the Lorentz factor of an adiabatically expanding jet scales as $`r`$. If $`\theta _0`$ decreases with $`r`$ inside of $`2000\mathrm{km}`$ as is likely, $`\mathrm{\Gamma }`$ grows more gradually. For the parameters chosen in Table 1, the terminal Lorentz factors of the jets if they expanded freely to infinity would be 150 for Models JA and JB. We also explored, for comparison, Model JC which had a lower initial Lorentz factor and a smaller total power, but a larger fraction of internal energy (Table 1). Fortunately, as we shall see, the initial partition of energy is not an important parameter of the model - so long as the total energy per baryon is large.

## 3 RESULTS

Because of the difficulty of carrying too large a range of radii on the computational grid, each calculation was carried out in two stages: 1) propagation inside the star, and 2) propagation in the stellar wind after breakout. The calculations, although two-dimensional, did not include the dynamical effects of rotation in the jet.

### 3.1 Jet Propagation Inside the Star

In all models the jet initially propagates along the polar axis. After a short time, it consists of a supersonic beam; a cocoon consisting of shocked jet material and shocked medium gas; a terminal bow shock; a working surface; and backflows (Blandford & Rees, 1974). Some snapshots of the jet for the three models JA, JB and JC, as defined in Table 1, are shown in Figures 1, 2 and 3. In all models, the jet is narrowly collimated and its beam is very thin, but other aspects of jet morphology vary from model to model. In Model JA, because of its large initial opening angle, most of the jet material passes, early on, through a strong shock, which deflects the flow toward the axis. In Model JB, which has a smaller initial opening angle and denser beam, the momentum flux is larger than for Model JA. The jet head thus advances faster and the beam is almost naked with a very thin cocoon. In Model JC, the thickness of the cocoon and beam and the velocity of the jet head are intermediate between JA and JB. At early times, the bow shock has a narrow head and wide tail (Fig. 3a), reflecting the jet’s acceleration. At late times however, especially as the jet nears the surface of the star, the morphologies of all three models are very similar (Fig. 1b, 2b and 3b). This suggests that the initial conditions tend to be forgotten as the jets propagate and interact with the star. It also implies that the properties of the subsequent GRB may not be sensitive to details of the operation of the central engine.

Because of the ram pressure of infalling stellar material at $`2000\mathrm{km}`$, a jet cannot start immediately. It must build up some pressure on the grid. In Model JA, 0.6 second elapses (in addition to the $`10\mathrm{s}`$ since core collapse) before the jet starts to propagate appreciably along the polar axis. The corresponding delays are 0.1 second and 0.3 second for Models JB and JC, respectively. Once they get going, the times for the jets to traverse the star (outside of $`2000\mathrm{km}`$) and break out also vary: 6.9, 3.3 and 5.5 seconds for Models JA, JB and JC, respectively (Table 2). This is reasonable since the jet in Model JB has a smaller cross section area and a larger momentum flux than in Models JA and JC.

Along the polar axis, the beam in each model is divided into two regions: an unshocked region and a shocked region (Fig. 4, 5 and 6). Because of its low density, low pressure, and relativistic velocity, the jet initially creates a low-density, low-pressure funnel along the axis. Within this funnel, the jet is steady and expands adiabatically without internal shocks. The high Lorentz factor helps to suppress the Kelvin-Helmholtz instability.

Farther out though, the jet head makes its way through the star at a sub-relativistic speed of roughly $`c/2`$ (this speed increases as the jet head moves outward in radius). Cocoon material flows back from the working surface. Relativistic material moving along the axis runs into comparatively stationary stellar matter at the jet head and pressure grows in this backed up material until it becomes comparable to that in the cocoon, which is in turn somewhat greater than the surrounding stellar medium. Then lateral expansion can occur. The need for deceleration propagates inwards in jet mass - though its location still moves outward in radius - as a reverse shock. At this inner shock the adiabatically expanding jet has a ram pressure, $`\rho h\mathrm{\Gamma }^2c^2`$, which approximately balances the pressure in the shocked jet, here $`h`$, the specific enthalpy, is $`(e+p)/(\rho c^2)`$, where, $`e`$, the total energy density, includes rest mass density.

To summarize, as seen by a piece of jet starting at the origin: 1) internal energy is converted to expansion until the reverse shock is encountered; 2) kinetic energy is then converted mostly back into internal energy at the reverse shock, though the motion remains moderately relativistic ($`\mathrm{\Gamma }510`$); and 3) the hot jet is further decelerated to subrelativistic speeds at the jet head. In what follows we shall refer to jet material in stage 1 as the “unshocked jet” and stage 2 as the “shocked jet”.

The interaction of the shocked jet with the star is especially interesting since Kelvin-Helmholtz instabilities and oblique shocks inside the cocoon can imprint time structure on its Lorentz factor (see also Martí et al. 1997 and references therein). The consequences differ in the three models. Because of its larger initial opening angle, the jet in Model JA experiences a stronger deceleration at the reverse shock, so the average Lorentz factor of the shocked jet in Model JA is lower than for Models JB and JC. In Model JA, the head of the unshocked jet (i.e., the location of the reverse shock) moves outward with a gradual acceleration. In Model JC the speed is almost constant, but, in Model JB, the location of this shock even moves backwards sometimes.

A possible explanation is that the unshocked jet in Model JB is unstable to pinching modes of the Kelvin-Helmholtz instability. It has been shown numerically (Martí et al., 1997) and analytically (Ferrari, Trussoni, & Zaninetti, 1978) that ultra-relativistic jet beams tend to be unaffected by the Kelvin-Helmholtz instability. However, small scale perturbations might still be able to grow if the jet beam is very thin (Ferrari, Trussoni, & Zaninetti, 1978). It is also important that these other works assumed equal pressures in the beam and external medium. Here the jet is under-pressurized compared to the medium (Fig. 7). When the mixing between jet beam and nearly stationary stellar material occurs, the beam will be decelerated. This behavior imprints time structure on the Lorentz factor which consequently turns out to be more variable Model JB than Models JA and JC. Some implications of this will be discussed in § 4.

As the jets pass through the star, they are narrowly collimated by the external pressure (Fig. 8). Near the origin, opening angles of jet beams are close to their initial values, $`20\mathrm{°}`$, $`5\mathrm{°}`$, and $`10\mathrm{°}`$ for Models JA, JB, and JC, respectively. In the unshocked region, the jet beam forms a low-pressure funnel that is quickly focused by the high pressure of the external stellar matter. In Model JA at $`t=12.0\mathrm{s}`$ for example, the half opening angle of the jet beam has already decreased to $`10\mathrm{°}`$ at $`r=8.2\times 10^9\mathrm{cm}`$, only $`1/10`$ of radius of the star (Fig. 8).

In the shocked jet the pressure is larger than in the external medium. The jet beam is surrounded by an over-pressurized cocoon which is in turn in balance with the star. There are no longer dramatic changes in $`\theta _0`$ in this shocked region. (Fig. 8). At the surface of the star, the half opening angles of the jet beams have decreased to very similar values, $`4\mathrm{°}`$, $`2\mathrm{°}`$, and $`3\mathrm{°}`$, for Models JA at $`t=12.0\mathrm{s}`$, JB at $`t=8.0\mathrm{s}`$, and JC at $`t=10.0\mathrm{s}`$, respectively. Still there are interesting differences among the models. Model JA which started with the largest initial angle still has the greatest angle at breakout. On the other hand, the energy flux per solid angle in Model JB is higher. Thus along the jet it would appear brighter. The distribution of energy with angle is also interesting. In Model JA, the jet is “hollow”, but in Models JB and JC, peak flux occurs on axis. Fig. 9 also shows that the opening angle in Model JA is increasing as the jet blows off the star. All of these effects, different intensities per solid angle, more rapidly variable Lorentz factor in Model JA, changing opening angle with time, etc. will have important consequences for the GRB light curves and its afterglows.

### 3.2 The Emergence of the Jet

Eventually, the jet accelerates and breaks free of the star. In Model JA, the average velocity of the head of the jet just prior to breakout is $`1.2\times 10^{10}\mathrm{cm}\mathrm{s}^1`$. For Models JB and JC the speeds are $`2.5\times 10^{10}\mathrm{cm}\mathrm{s}^1`$ and $`1.5\times 10^{10}\mathrm{cm}\mathrm{s}^1`$ respectively. Of course the speeds behind the jet head and in the beam are higher, $`\mathrm{\Gamma }`$ 10, 50, and 20 for Models JA, JB and JC, respectively (Fig. 10). These would be regarded by most as too low to produce healthy GRBs (e.g., Lithwick & Sari 2001). However, these jets still have a large ratio of internal energy to total energy: $``$ 80-95%, $``$ 50-75%, and $``$ 80-90% for Models JA, JB and JC. That is the jets themselves are still fireballs with $`\eta `$ 3 to 10. After all internal energy is converted to kinetic energy, the terminal Lorentz factor will still be $``$ 100-200. This implies that despite being shocked deep inside the star and experiencing Kelvin-Helmholtz instabilities, the matter that emerges has almost the same energy per baryon that it started out with at $`2000\mathrm{km}`$. The baryonic entrainment has been small. Of course this result needs to be tested in a three dimensional calculation, but it is encouraging.

After breakout, a channel is open through which the jet can propagate without losing much energy (though it is still shocked and modulated while passing through the star). The total power carried by matter with Lorentz factor greater than 10 is $`0.9\times 10^{51}\mathrm{erg}\mathrm{s}^1`$, $`1\times 10^{51}\mathrm{erg}\mathrm{s}^1`$ and $`3.0\times 10^{50}\mathrm{erg}\mathrm{s}^1`$ just after breakout for Models JA, JB and JC, respectively (Fig. 11). Those are very close to the initial values at the inner boundary (Tab. 1).

By the time of breakout, the total energy injected as twin jets was 13.8, 6.6, and $`3.3\times 10^{51}\mathrm{erg}`$ in Models JA, JB, and JC respectively. The energy on the grid in various forms is given at that time in Table 2. In particular we can estimate the available energy to produce a supernova, $`E_{\mathrm{sn}}`$, from the total energy in matter moving at less than 10% the speed of light ($`\mathrm{\Gamma }<1.005`$) by considering corrections due to initial energy and energy lost at the computational inner boundary. The jet propagates very efficiently through the star after breakout (Fig. 11) so that there is very little further “wasted” energy available for powering a supernova in any of the models. However, one should keep in mind a) that some non-relativistic energy may still be available to power the supernova as the jet starts to die off, and b) the wind off the accretion disk may also contribute significant energy to the explosion (MacFadyen & Woosley 1999). The supernova energies in Table 2 are thus lower bounds for the (admittedly high) jet energies assumed.

We followed Models JA, JB, and JC for a total of $`21.8\mathrm{s}`$, $`10.0\mathrm{s}`$, and $`15.0\mathrm{s}`$ respectively (Table 3) and their total energies were these times multiplied by the assumed injection rates (Table 1). Before breakout some energy is “wasted” because the head of the jet moves much slower than the jet itself. The amount of “wasted” energy is given by, $`E_\mathrm{c}=2\dot{E}(t_\mathrm{b}R_\mathrm{s}/c)`$, where $`\dot{E}`$ is the power at the base for each jet, $`t_\mathrm{b}`$ is the time for the jet to break out of the star, and $`R_\mathrm{s}`$ is the stellar radius. Using this expression, we estimate the sums of the energies in the supernova and jet cocoon (Table 3). Most of this energy is in non-relativistic hot material near the jet which can power the supernova; the remainder is in the cocoon which emerges from the star (Ramirez-Ruiz, Celotti, & Rees, 2002).

The history of Lorentz factor at the edge of our computational grid on the polar axis is shown in Fig. 10. This figure also shows the estimated terminal Lorentz factor if all internal energy is converted into kinetic energy and the process is adiabatic. The fluctuations of the Lorentz factor are not so large as to make efficient internal shocks, but our grid is relatively coarse and numerical viscosity may wash out most fluctuations, especially those due to small scale instabilities. Because of the large range of radii covered, our two-dimensional spherical grid consists of logarithmically spaced radial zones, and both uniform and logarithmically angular zones. This gives good resolution at the center, but unfortunately the resolution near the edge of the star is too low to study the details of time structure. The Lorentz factor along the polar axis at $`4.0\mathrm{s}`$ in Model B contains many significant fluctuations (Fig. 5), but these are smoothed quickly by numerical viscosity. In the near future we plan further calculations using cylindrical geometry and finer resolution to study the mixing better.

When the jet breaks out, it is the hot shocked beam that emerges first. If the engine lasts long enough, the unshocked beam may eventually emerge, e.g. at $`20`$ s in Model JA (Fig. 4; remember that $`10\mathrm{s}`$ must be added to all our times to get the time when the core actually collapsed). In Model JB, the mixing instability caused by the interaction among the jet beam, cocoon and star slows the jet more. About 6.7 seconds after breakout, the head of the unshocked jet is still only at at $`2.5\times 10^9\mathrm{cm}`$ while in Model JC it is at $`7.9\times 10^9\mathrm{cm}`$. Because more than half of the total energy in the shocked jet material is in form of internal energy, it is inevitable that the jet will experience acceleration and, possibly, sideways expansion. What is the final opening angle of the jet? What is the final Lorentz factor of the jet as the internal energy is converted into kinetic energy? In order to power a GRB, the jet has to have a Lorentz factor of more than 100. Can the shocked mildly relativistic jet ($`\mathrm{\Gamma }1050`$) acquire high enough Lorentz factor ($`\mathrm{\Gamma }>100`$) after breakout?

### 3.3 Jet Propagation in the Near Stellar Environment

After it breaks out of the star, the relativistic jet, which has a high internal energy loading, accelerates and spreads. It also encounters the stellar wind surrounding the star. In a second series of calculations we explored, in a preliminary way, this expansion. The jet is once more specified by its energy deposition rate, $`\dot{E}`$, initial opening angle, $`\theta _0`$, initial Lorentz factor, $`\mathrm{\Gamma }_0`$, and the ratio of kinetic energy to total energy, $`f_0`$. Two models were considered (Table 4): (W1) $`\dot{E}=8\times 10^{50}\mathrm{erg}\mathrm{s}^1`$, $`\theta _0=3\mathrm{°}`$, $`\mathrm{\Gamma }_0=10`$, $`f_0=0.06`$; (W2) $`\dot{E}=8\times 10^{50}\mathrm{erg}\mathrm{s}^1`$, $`\theta _0=3\mathrm{°}`$, $`\mathrm{\Gamma }_0=50`$, $`f_0=0.33`$. These conditions are representative of the actual conditions in the emerging jets for Models JA, JB, and JC.

Because, for this set of calculations, we are not so interested in what happens close to the stellar surface, a two-dimensional cylindrical grid ($`r,z`$) is employed consisting of 2400 zones in $`z`$-direction from $`10^{11}\mathrm{cm}`$ to $`1.3\times 10^{12}\mathrm{cm}`$ and 720 zones in $`r`$-direction from $`0`$ to $`3.6\times 10^{11}\mathrm{cm}`$. Due to the relativistic velocities and large radius we can ignore gravity. A gamma-law equation of state with $`\gamma =4/3`$ also suffices. The background density, which might come from stellar mass loss, has been assumed to decline as $`\rho R^2`$, and the density at $`R=10^{11}\mathrm{cm}`$ is set to $`3.7\times 10^{11}\mathrm{g}\mathrm{cm}^3`$, which roughly corresponds to a mass loss of $`10^5M_{\mathrm{}}\mathrm{yr}^1`$ with a velocity of $`1000\mathrm{km}\mathrm{s}^1`$. Actually, the details of environment are unimportant here because we have not yet followed the jet to such large radii that the interaction with the wind becomes important. The jet is injected at $`z=10^{11}\mathrm{cm}`$ with an opening angle of $`3\mathrm{°}`$ with momentum that is purely radial. Though potentially quite important, to keep things simple, we do not include a cocoon in these calculations. Also for simplicity, the jet is given uniform initial properties within its opening angle.

In our simulations, the jets are injected with the constant properties specified above for the first 10 seconds. Eventually though, we expect the central engine to decline in power and produce a less energetic jet. To mimic this decline, the jet energy and Lorentz factor were both turned down during the interval 10 to $`20\mathrm{s}`$ while the pressure and density of the jet were kept constant. The energy deposition rate decreased linearly from $`8\times 10^{50}\mathrm{erg}\mathrm{s}^1`$ at $`t=10\mathrm{s}`$ to 0 at $`t=20\mathrm{s}`$ (Fig. 12). After 20 seconds, we used an outflow (zero gradient) boundary condition for the lower $`z`$ boundary. The total energy for each jet is thus $`1.2\times 10^{52}\mathrm{erg}`$ for both Models W1 and W2.

Some results are shown in Figs. 13 and 14. In both models, the Lorentz factor increases to more than 100 as the jet propagates in the stellar wind. There are very significant differences between Models W1 and W2. Although the initial opening angle in both is $`3\mathrm{°}`$, the opening angle at the head of the jet at $`t=35\mathrm{s}`$ is $`15\mathrm{°}`$ for Model W1 and $`5\mathrm{°}`$ for Model W2 (Fig. 15). In Model W1 the jet initially had a Lorentz factor of 10 and a large internal energy (93% of the total energy) leading to a large lateral expansion. In Model W2, the internal energy is less (67% of the total energy) and the Lorentz factor greater (50). One might argue that the final opening angle of the jet should be $`\theta _0+1/\mathrm{\Gamma }_0`$ (see, e.g., Matzner 2002) in which case the final opening angles would be $`8.7\mathrm{°}`$ and $`4.1\mathrm{°}`$ for Models W1 and W2, respectively. Our larger angles correctly reflect the dynamics of a differentially accelerating relativistic jet.

As the jet expands sideways, the velocity along its edge does not increase as rapidly in the radial direction as at the center. A small amount of material actually ends up, briefly, with a slower radial component. The conversion of internal energy to kinetic energy during this phase causes sideways expansion. It should be noted that the radial component of momentum still increases when the radial component of velocity decreases because the total Lorentz factor is increasing during the sideways expansion. The combined effect of these two processes is that a small amount of material ends up at larger angles than expected from $`\theta _0+1/\mathrm{\Gamma }_0`$. In Model W1, because of its lower initial Lorentz factor and higher internal energy loading, this effect is more pronounced and the opening angle increases from $`3\mathrm{°}`$ to $`>6\mathrm{°}`$ during the first second. The effect becomes less important as the Lorentz factor grows. In contrast, the final opening angle in Model W2 is only slightly larger than $`\theta _0+1/\mathrm{\Gamma }_0`$ because of its higher initial Lorentz factor and lower internal energy loading. In order to test whether this large sideways expansion is due to low numerical resolution, we ran two higher resolution tests for Model W1. Finer grids (up to 3.3 times the original resolution of Model W1) were run for the same conditions. In one of the tests, the density of stellar wind was also increased by two orders of magnitude. These studies demonstrated that the resolution was already high enough to ensure the numerical viscosity contributed very little to the sideways expansion but are somewhat sensitive to the density of stellar wind (Fig. 16). Because of sideways expansion, the energy flux far off axis in Model W1 is small but not zero. In Model W2, the profile of energy flux is more uniform and has a sharper edge. In both cases the Lorentz factor profile is flat inside the jet in both models.

At later times, 10 to 20 seconds after the jet was initiated outside of the star, the energy deposition rate and Lorentz factor were gradually turned off as described previously. As a result, the tail of the jet experiences much more lateral expansion. Although the details depend on exactly how the central engine shuts down, it is natural that the jet should go a period of gradual decline rather than abruptly shutting off. Fig. 15 shows the energy flux and Lorentz factor at different locations for the two models. Mildly relativistic material is ejected at large angles. At the end of our simulations, the jets have two components: a core with high Lorentz factor and a wide tail with low Lorentz factor. In Model W1, the core has a rather flat energy flux profile inside $`3\mathrm{°}`$ and a rapidly decreasing energy flux profile extended to $`15\mathrm{°}`$. There is also a moderate relativistic ($`\mathrm{\Gamma }10`$) tail outside the high Lorentz factor core. The Lorentz factor of the core is about $`60100`$, but it has not reached its saturated value and is expected to become $`150`$ because of the remaining internal energy. In Model W2, the jet has a narrow core of $`5\mathrm{°}`$ with a very flat energy flux profile. It also contains moderately relativistic material outside the high Lorentz factor core. The implications for GRBs will now be discussed.

## 4 OBSERVATIONAL IMPLICATIONS

### 4.1 Beaming

Our special relativistic calculations show, as did Aloy et al. (2000), that a jet originating near the center of a collapsing massive star will emerge with an opening angle $`5`$ degrees. This result is relatively insensitive to the initial opening angle. Further expansion occurs after the jet has already emerged from the star. Depending on its Lorentz factor and internal energy loading at breakout, the final opening can be smaller than 5 degrees, or more than 15 degrees. This suggests that gamma-ray bursts can have a wide range of opening angle, and can be beamed to less than 1% of the sky with obvious implications for the requisite GRB energy and event rate.

### 4.2 GRB Light Curves

For a typical total energy in relativistic matter ($`\mathrm{\Gamma }>`$ a few) of $`3\times 10^{51}\mathrm{erg}`$ (Frail et al., 2001; Freedman & Waxman, 2001), and an assumed Lorentz factor, $`\mathrm{\Gamma }100`$ (Lithwick & Sari, 2001), the mass within the GRB-producing beam is $`10^5M_{\mathrm{}}`$, or about 1000 times less than that within the same solid angle in the precollapse star. In order to make a GRB using a jet launched from the stellar center, this stellar matter must be pushed aside. The matter shoved aside has two consequences. First, its lateral motion initiates a shock wave that moves round the star converging on the equator and exploding at least the outer part of the star as a supernova. Second, some of the material is swept back creating the cocoon of the jet. Instabilities between the cocoon and jet beam modulate the flow leading to a variable Lorentz factor. It is even possible that the jet might be temporarily clogged and reopened by the flow.

In this fashion, the variable Lorentz factor needed for the internal shock model for GRBs is generated. Light curves having very complex temporal behavior can result in situations where the central engine produced only a jet with constant power. The converse is also true. Short term time structure (less than the jet transit time through the star or about $`5\mathrm{s}`$) of the central engine may be erased as the jet passes through the reverse shock. If the variability of the engine itself is visible at all, it may only be at very late times when the star has exploded.

Our calculations show significant variation in the Lorentz factor of the emergent jet (Figs. 4, 5, 6 and 10). We suspect these variations may be larger in calculations with higher resolution, and in three dimensions.

Our model also suggests that the variability of GRB light curves may evolve with time. Though it takes awhile, of order a minute, for the star to explode away from the jet opening its channel appreciably, the Kelvin-Helmholtz instabilities should decrease with time. This implies that GRB variability may also decrease with time. A narrower jet experiences more instabilities and therefore the GRB is more variable. For a given total energy input at the bottom, a narrower jet also carries more energy per solid angle. Perhaps this might help to explain the observed correlation between luminosity and variability (Fenimore & Ramirez-Ruiz, 2000; Reichart et al., 2001).

### 4.3 A Unified Model For High Energy Transients

According to the “Unified Model” for active galactic nuclei (e.g., Antonucci 1993), one sees a variety of phenomena depending upon the angle at which a standard source is viewed. These range from tremendously luminous blazars, thought to be jets seen on axis, to narrow line radio galaxies and Type 2 Seyferts thought to be similar sources seen edge on. Given that an accreting black hole and relativistic jet may be involved in both, it is natural to seek analogies with GRBs.

In the equatorial plane of a collapsar - the common case - probably little more is seen than an extraordinary supernova. Along the axis, one sees an ordinary GRB, but there may also be interesting phenomena at intermediate angles. The calculations presented here clearly show that the edges of jets are not discontinuous surfaces. Moving off axis, one expects and calculates a smooth decline in the Lorentz factor and energy of relativistic ejecta. These low energy wings with moderate Lorentz factor come about in four ways. First, the jet that breaks out still has a lot of internal energy. Expansion of this material in the co-moving frame leads to a broadening of the jet. As a result a small amount of material with low energy ends up moving with intermediate Lorentz factors - say 10 - 30 and at angles up to several times that of the main GRB-producing jet. Second, as the star explodes from around the jet, the emerging beam opens up. Even at constant total power, the power per solid angle will decline as this channel broadens. Third, the jet is surrounded by a hot mildly relativistic cocoon. This material has low energy and low Lorentz factor. It can expand to large angles. Fourth, and potentially most important, the jet can continue with a declining total power for a long time. It is natural that the Lorentz factor of the emerging jet decline as well. A hot jet emerging with lower Lorentz factor will expand to larger angles (§ 3.3). As a consequence of the spreading of the jet, a large region of the sky, much larger than that which sees the main GRB, will see a hard transient with less power, lower Lorentz factor, and perhaps coming from an external shock instead of internal ones.

We have speculated for some time now (Woosley, Eastman, & Schmidt, 1999; Woosley & MacFadyen, 1999; Woosley, 2000, 2001) that these ordinary bursts seen off axis might appear as hard x-ray transients of one sort or another. We have identified them with GRB 980425 and with the class of hard x-ray flashes reported by Heise et al. (2001). These are not ordinary high $`\mathrm{\Gamma }`$ jets seen just beyond a sharp edge (Ioka & Nakamura, 2001). The events are made by matter moving toward us.

### 4.4 Breakout Transients and Short Hard Bursts

The breakout of the jet and its interaction with the stellar wind will surely lead to some sort of “precursor” event. In fact, it is possible that this event could, in some cases, dominate the display.

The total energy in a short hard GRB (SHB) is, on the average, about 1/50 that of a long-soft burst (LSB). If the total gamma-ray energy of the latter averages $`3\times 10^{50}\mathrm{erg}`$, then SHBs at the same distance have gamma-ray energy, including the effects of beaming, of $`10^{49}\mathrm{erg}`$, or $`3\times 10^{51}\mathrm{erg}`$ in equivalent isotropic energy.

Our present calculations lack adequate resolution at the surface of the star to follow shock breakout with any accuracy. They also have not been followed sufficiently long to see the interaction of the leading head of the jet with the circumstellar medium. However, we do focus approximately $`5\times 10^{51}\mathrm{erg}`$ on only 1/800 of the stellar surface (for a fiducial opening half-angle of 4 degrees and a jet that takes $`5\mathrm{s}`$ to exit the star). This corresponds to an equivalent isotropic energy of $`4\times 10^{54}\mathrm{erg}`$. For their fiducial model with total isotropic energy $`5\times 10^{54}\mathrm{erg}`$, Tan, Matzner & McKee (2001) calculate about $`3\times 10^{51}\mathrm{erg}`$ is accelerated to Lorentz factor $`\mathrm{\Gamma }>40`$. It may be that a larger amount of material is carried along by the leading edge of the jet as it expands and accelerates outside the star (see also Waxman & Mészáros 2002) if the lateral expansion of the “plug” cannot keep pace with the radial expansion of the jet.

Material having $`3\times 10^{51}\mathrm{erg}`$ in kinetic energy and $`\mathrm{\Gamma }40`$, has a mass $`4\times 10^5`$ $`M_{\mathrm{}}`$ which will decelerate after running into $`10^6M_{\mathrm{}}`$. If the mass loss rate of the star during its last days was $`10^4M_{\mathrm{}}\mathrm{yr}^1`$ and the velocity $`10^8\mathrm{cm}\mathrm{s}^1`$, this deceleration will occur within $`3\times 10^{13}\mathrm{cm}`$ or $`1000\mathrm{s}`$ in the lab frame. The duration of the burst is shortened by $`\mathrm{\Gamma }^2`$ to be less than a second, that is it may be a short one. But why then don’t we see such short hard precursors in all GRBs and how can one get short hard bursts in isolation?

The absence of SHBs in LSBs could, to some extent, be semantic, If the burst continues, a spike at the beginning may be counted as part of a long burst. But what about the converse problem - SHBs seen in isolation? The relative brightness of the breakout burst to the subsequent burst depends on many factors - the circumstellar density, the total energy of the jet, the variation of its Lorentz factor (especially critical in the internal shock model), and the mass of the star. This is especially so since SHBs could be made by external shocks and LSBs by internal shocks. So it may be, occasionally, that the breakout burst is an order of magnitude brighter (in flux, if not fluence) than the ensuing long burst. Of course, not all SHBs need to be produced this way and other models like merging neutron stars remain viable.

This possibility suggests that a careful search for postburst activity might be productive in SHBs. Preliminary results by Connaughton (2002) are suggestive of continuing postburst emission on a time scale not unlike that of LSBs, but are not sufficiently sensitive for a detailed comparison.

### 4.5 Post-burst Phenomenology

After the main burst is over, accretion will continue at a decaying rate. The lateral shock launched by the jet starts at the pole and wraps around the star, but does not reach into the origin at the equator. One may envision an angle-dependent “mass cut”. Consequently, some reservoir of matter remains to be accreted at late time. This accretion occurs at a rate governed by the viscosity of the residual disk and the free fall time of material farther out not ejected in the supernova. MacFadyen, Woosley, & Heger (2001) and Chevalier (1989) estimate the accretion rate from fall back to be $`3\times 10^6t_4^{5/3}M_{\mathrm{}}\mathrm{s}^1`$. Here $`t_4`$ is the elapsed time since core collapse in units of $`10^4\mathrm{s}`$. Given the slow rate, the disk that forms is not neutrino dominated and there may be considerable high velocity flow from its surface (MacFadyen & Woosley, 1999; Narayan, Piran, & Kumar, 2001). That outflow will still be jet-like in nature since the equatorial plane is blocked by the disk and its energy will be $`5\times 10^{46}t_4^{5/3}ϵ_2\mathrm{erg}\mathrm{s}^1`$ where $`ϵ_2=ϵ/10^2`$ is the efficiency for converting rest mass into outflow kinetic energy measured in percent. This is comparable to the energy in x-ray afterglows and might be important for producing the emission lines reported in some bursts (Rees & Meszaros, 2000; Mészáros & Rees, 2001; McLaughlin et al., 2002) and for providing an extended tail of hard emission in the GRB itself.

As a consequence of this continuing outflow, the polar regions of the supernova made by the GRB remain evacuated and the photosphere of the object resembles an ellipse seen along its major axis but with conical sections removed along the axis. An observer can see deeper into the explosion than they could have without the operation of the jet’s “afterburner”.

## 5 FUTURE WORK

The present study, and its precursors by MacFadyen & Woosley (1999) and Aloy et al. (2000), have shown the critical role of the jet-star interaction in producing the observable properties of GRBs. Further work will be needed though to confirm the validity of the present results and add strength to some of the model’s more speculative aspects (§ 4).

First, it is important to repeat our calculations in three dimensions. Use of two dimensions imposes a rotational symmetry that makes it difficult for instabilities to deflect the beam. It may be that jets which maintain a stable focus here in 2D will waiver and even disperse in 3D. In any case, the Kelvin-Helmholtz instabilities may behave differently. Fortunately the computational requirements are not too extreme and such a calculation can be carried out in the near future.

Second, our jet needs to be studied with higher resolution. One can do that using cylindrical coordinates and a smaller range of radii. Again, we need to know the nature of the relativistic Kelvin-Helmholtz instability better.

The emergence of the jet and its immediate interaction with the circumstellar medium needs further examination, especially in a calculation that finally resolves the stellar surface at breakout and includes the jet cocoon as well as its core (see also Ramirez-Ruiz et al. 2002).

Finally, we are interested in the long term evolution of the jet and the star it is exploding. How much of the star remains to accrete even after the initial shock from the jet has wrapped around the star? How will the object appear, not seconds, but hours and days after the burst? How does this relate to the possibility of x-ray lines being detected in GRB afterglows.

Calculations to address all of these issues are underway.

We appreciate many helpful conversations on the subject of GRBs with Chris Fryer, Alex Heger, Chryssa Kouveliotou, Chris Matzner, and Martin Rees and are particularly indebted to Dr. Heger for providing the unpublished precollapse model used in this calculation. We thank Enrico Ramirez-Ruiz and M. A. Aloy for helpful comments on the manuscript. We also thank the anonymous referee for useful comments. This research has been supported by NASA (NAG5-8128, NAG5-12036, and MIT-292701) and the DOE Program for Scientific Discovery through Advanced Computing (SciDAC; DE-FC02-01ER41176).


format: text

text:
# Modeling Brown Dwarfs, L Dwarfs, and T Dwarfs

## 1. Introduction

The subject of brown dwarfs (and substellar-mass objects in general) is entering a new phase of rapid expansion and discovery. More than 200 L dwarfs are now known, and they are joined by $``$30 T dwarfs. This development requires a corresponding expansion in theoretical effort, involving at its core evolutionary, spectral, and compositional modeling. Crucial to spectral modeling are molecular and atomic opacities, many of which have not been addressed before with the degree of seriousness and completeness that standard stellar atmospheres studies have long enjoyed. However, the pace of relevant spectral and opacity calculations is accelerating and in this spirit I summarize in this paper a few such topics of recent interest. In §2., I show new proof of the “rainout” and settling of refractory elements in T dwarfs. Rainout leaves as one of its consequences the lower-temperature reaches of T dwarf atmospheres enhanced in sodium and potassium atoms. In §3., I touch on some new calculations of the alkali line wings to many 1000’s of Å detunings. The red wing of the K I line centered at 0.77 $`\mu `$m in particular defines much of the T dwarf continuum between 0.77 $`\mu `$m and 1.0 $`\mu `$m . This is followed in §4. by a short discussion of the new CrH opacities recently generated by our group. Finally, and a bit tentatively, in §5. I finish with a short discussion on nomenclature, a subject that continues to exercise the brown dwarf and extrasolar planet communities.

## 2. A New Proof of Rainout/Depletion in T Dwarfs

As explained in Burrows, Marley, and Sharp (2000, BMS), rainout depletes the atmospheres of the refractory elements Ca, Al, Mg, Fe, and Si. In particular, rainout suppresses the formation of alkali feldspars and enables atomic Na and K to survive to lower temperatures and pressures, at which point they form Na2S(c) and KCl(c) (Lodders 1999). As a consequence, their influence on the emergent spectrum in the “optical” between 0.5 $`\mu `$m and 1.1 $`\mu `$m is enhanced. This effect is seen dramatically for the first time in Fig. 1. At $`0.7\mu \mathrm{m}`$, the data for Gliese 570D (Burgasser et al 2001), the coolest known T dwarf, clearly favor the “rainout” model. This fact reiterates the conclusion of BMS that the wings of the Na D and K I (7700Å) resonance lines and aggressive rainout of heavy metals (with the resulting enhancement of the sodium and potassium abundances at altitude) are required to fit the T dwarf data shortward of 1.0 $`\mu `$m . It is shortward of $``$1.0 $`\mu \mathrm{m}`$ that the alkali chemistry can be most readily probed (Burrows et al. 2002a).

## 3. Alkali Metal Line Profiles

To continue on the alkali theme, in standard stellar atmospheres, atomic lines are superposed on a background continuum and the concepts of equivalent width and curve of growth make conceptual and practical sense. An individual line is but a perturbation on the local spectrum. However, due to the rainout of metals in cool atmospheres and the consequent paucity of continuum and alternate opacity sources between 0.4 $`\mu `$m and 1.0 $`\mu `$m , the wings of the strong resonance doublets centered at $``$7700 Å (K I) and 5890 Å(Na-D) assume the role of the continuum throughout most of this spectral range (BMS; Tsuji, Ohnaka, and Aoki 1999). In particular, the red wing of the $`4s^2S_{1/2}4p^2P_{3/2}`$ transitions of K I provides the pseudo-continuum in cool molecular atmospheres all the way from 0.77 $`\mu `$m to $``$1.0 $`\mu `$m . The Na-D doublet, centered as it is in the middle of the visible, determines the true color of brown dwarfs (magenta/purple; Burrows et al. 2001).

Hence, whereas in traditional stellar atmospheres the Lorentzian core and Gaussian wings of a line are not generally of relevance beyond $``$20 Å detunings ($`\mathrm{\Delta }\lambda `$ from the line core), in cool substellar atmospheres the relevant reach of the Na I and K I resonance lines can be thousands of Å. Given this, to achieve accurate spectral fits for brown dwarf, L dwarf, T dwarf, and hot giant planet atmospheres, the shapes of the far wings of these alkali lines as a function of pressure and temperature must be ascertained.

Recently, Burrows and Volobuyev (2002c) performed ab initio calculations of the energy shifts of the ground and excited states of sodium and potassium immersed in H2\- and helium-rich atmospheres and have obtained the opacity profiles of the red and blue wings of the Na-D and K I resonance lines using the Unified Franck-Condon formalism (Szudy and Baylis 1975,1996). The absorption spectra as a function of photon wavelength for the D2 and D1 lines of Na and K centered near 5890 Å and 7700 Å, respectively, are depicted in Fig. 2 at a temperature of 1000 K and a total pressure of one atmosphere. The contributions at various orientation angles of H2, the angle-integrated spectra (solid), and the results for the Na+He and K+He systems are shown. The Lorentzian cores are not included on this plot.

As seen in Fig. 2, the cutoff on the red wing of the potassium feature due to an exponential term containing the ground-state interaction potential of the K+H2 system is seen to be situated between 0.95 $`\mu `$m and 1.0 $`\mu `$m , close to the 0.98 $`\mu `$m used in Burrows et al. (2002c). The cutoff of the red wing of the Na feature is near 0.8 $`\mu `$m , but is a bit more gradual. The corresponding cutoffs for the alkali-He systems are more abrupt, but of less importance due to the lower abundance of helium. Importantly, unlike in the algorithm of BMS, there are no free parameters for this theory of line profiles; the cutoff on the red wing of the K I feature appears naturally. Fig. 2 also makes clear that the red and blue wings are asymmetrical.

## 4. New CrH Opacities

The 1-0 and 0-0 bands of the A$`{}_{}{}^{6}\mathrm{\Sigma }_{}^{+}`$ – X$`{}_{}{}^{6}\mathrm{\Sigma }_{}^{+}`$ transition of CrH are used as primary markers for the new L dwarf spectral class (Kirkpatrick et al. 1999a,b). Hence, accurate CrH line lists and oscillator strengths are needed to calculate the CrH opacities now used to model the spectral energy distributions of these transitional and substellar objects.

Burrows et al. (2002b) have recently created a program for generating CrH cross sections for the six band sequences with band heads from $``$0.7 $`\mu `$m to $``$ 1.4 $`\mu `$m . Figure 3 depicts the resulting CrH cross sections for representative pressures and temperatures at which CrH is typically found in substellar atmospheres. Specifically, opacity spectra for $`T/P`$ pairs of 1500 K/10 bars, 2000 K/10 bars, and 2000 K/100 bars are portrayed. The corresponding opacities using the older database for the 0–0 transition of CrH of Ram, Jarman, and Bernath (1993), as calculated by R. Freedman (private communication), are two orders of magnitude weaker. The lion’s share of the difference between the old and the new CrH opacities can be traced to an increase by a factor of $``$13.5 in the oscillator strength and to a previously inappropriate division by 6, the electronic spin degeneracy factor. The A–X 0–0 band at $``$0.86 microns (just shortward of the neighboring FeH feature at $``$0.87 microns) and the A–X 0–1 CrH band near $``$0.997 microns (just longward of the classic Wing-Ford band of FeH) are prominent and, in principle, diagnostic features in measured L dwarf spectra (Kirkpatrick et al. 1999a,b).

Since CrH is a defining molecule of the L dwarf spectroscopic class, accurate opacities as a function of temperature and pressure are necessary for spectral syntheses and to extract CrH abundances for L dwarf atmospheres. Burrows et al. (2002b) have used the new theoretical opacity data to obtain an abundance for the L5 dwarf 2MASSI J1507038-151648 (Kirkpatrick et al. 1999b). Figure 4 depicts this comparison. The CrH/H2 number ratio they find is $`24\times 10^9`$, in reasonable agreement with expectations.

## 5. What’s in a Name?

The deuterium burning cut was used originally by some to distinguish planets from brown dwarfs. This was partly because of a figure in Burrows et al. (1997) that arbitrarily distinguished evolutionary trajectories in luminosity/age space on this basis, calling one class that burned tens of percent or more of their deuterium brown dwarfs and the others that did not “planets” (as I recall, the original figure had the quotes as well). Also, one has the unfounded prejudice that a planet should not have a thermonuclear phase.

One should distinguish giant planets from brown dwarfs by their mode of formation, by their origin. Since an object’s origin is not easily “detected” or discerned, the situation is ready-made for confusion. In fact, there is really no need to define and codify the nomenclature nor the name to be applied to any particular object or class of objects at this time; people just want (perhaps a bit too rigidly and with unjustified Procrustean zeal) to name things, or to stuff objects into cubbyholes.

There is ambiguity in the provenance of these objects and this ambiguity will be with us for a while. If one uses a particular term (planet or brown dwarf), one should always include the definition one is employing. Not to do so implies that a standard, acceptable definition has been arrived at; despite what some have written and the IAU, such is not the case.

Hence, formation history is a key, but not yet known. One would expect a planet (or most of them) to be orbiting a star (or brown dwarf). There may be metallicity (composition) differences between objects with the same mass, but different origins or modes of formation. Eccentricity is probably not a good discriminant (though this is not known). Very few “planets” are expected to be flung out of their cradles. Brown dwarfs are just the tail of the stellar mass distribution, but the mass distributions of brown dwarfs and giant planets may overlap.

When classifying a newly-discovered substellar-mass object, one can use a variety of “reasonable” criteria, but the tentative classification is just that, tentative. The usefulness of a classification scheme at this stage in the development of the two related subjects of brown dwarfs and giant planets is inversely proportional to its rigidity. A flexible and open-minded philosophy towards nomenclature is best, which more data and information will progressively guide towards a more and more reasonable naming scheme. It may be that the mass distributions of planets and brown dwarfs don’t overlap, that their compositions will enable one to distinguish one from the other, etc. Until then, I prefer a bit of ambiguity to the illusion of clarity and the pretense of certainty. When there is fog, to acknowledge its presence is more honest than to try to ignore it.

### Acknowledgments.

The author thanks Bill Hubbard, Jonathan Lunine, David Sudarsky, Ivan Hubeny, Christopher Sharp, Drew Milsom, Maxim Volobuyev, Curtis Cooper, Adam Burgasser, Davy Kirkpatrick, and Jonathan Fortney for informative exchanges and collaboration, as well as NASA for its financial support via grants NAG5-10760 and NAG5-10629.

## References

Anders, E. and Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53, 197

Burgasser, A.J., et al. 2002, ApJ, 564, 421

Burrows, A., Marley M., Hubbard, W.B. Lunine, J.I., Guillot, T., Saumon, D. Freedman, R., Sudarsky, D. and Sharp, C.M. 1997, ApJ, 491, 856

Burrows, A., Marley, M.S., and Sharp, C.M. 2000, ApJ, 531, 438

Burrows, A., Hubbard, W.B., Lunine, J.I., and Liebert, J. 2001, Rev. Mod. Phys., 73, 71

Burrows, A., Burgasser, A.J., Kirkpatrick, J.D., Liebert, J., Milsom, J.A., Sudarsky, D., and Hubeny, I. 2002a, ApJ,, 573, 394

Burrows, Ram, R, Bernath, P, Sharp, C., and Milsom, J.A. 2002b, ApJ, in press

Burrows, A. and Volobuyev, M. 2002c, submitted to ApJ

Kirkpatrick, J.D., Reid, I.N., Liebert, J., Cutri, R.M., Nelson, B., Beichman, C.A., Dahn, C.C., Monet, D.G., Gizis, J., and Skrutskie, M.F. 1999a, ApJ, 519, 802

Kirkpatrick, J.D., Allard, F., Bida, T., Zuckerman, B., Becklin, E.E. Chabrier, G., Baraffe, I. 1999b, ApJ, 519, 834

Lodders, K. 1999, ApJ, 519, 793

Ram, R.S., Jarman, C.N. and Bernath, P.F. 1993, J. Mol. Spectrosc., 161, 445

Szudy, J. and Baylis, W. 1975, J. Quant. Spectrosc. Radiat. Transfer, 15, 641

Szudy, J. and Baylis, W. 1996, Physics Reports, 266, 127

Tsuji, T., Ohnaka, K., and Aoki, W. 1999, ApJ, 520, L119


text:
# Size dependent tunneling and optical spectroscopy of CdSe quantum rods

## Abstract

Photoluminescence excitation spectroscopy and scanning tunneling spectroscopy are used to study the electronic states in CdSe quantum rods that manifest a transition from a zero dimensional to a one dimensional quantum confined structure. Both optical and tunneling spectra show that the level structure depends primarily on the rod diameter and not on length. With increasing diameter, the band-gap and the excited state level spacings shift to the red. The level structure was assigned using a multi-band effective-mass model, showing a similar dependence on rod dimensions.

Colloidal semiconductor nanocrystals are a class of nanomaterials that manifest the transition from the molecular limit to the solid state Alivisatos (1996); Millo et al. (2001), with significant potential for serving as building blocks of nano-devices in applications ranging from lasers Klimov et al. (2000); Kazes et al. (2002) and opto-electronic devices Tessler et al. (2002) to biological fluorescence tagging Bruchez et al. (1998). Shape control of such colloidally prepared nanostructures has been recently achieved by modifying the synthesis to obtain rod shaped particles \- quantum rods Peng et al. (2000). Quantum rods (QRs) exhibit electronic and optical properties different than QDs. For example, unlike the spherical dots, QRs have linearly polarized emission as demonstrated recently by fluorescence measurements on single rods Hu et al. (2001), leading to polarized lasing Kazes et al. (2002). In this letter we combine optical and tunneling spectroscopies to investigate the electronic level structure of CdSe quantum rods, and study its dependence on rod length and diameter. The levels are assigned with the use of a multi-band effective-mass model. The study provides significant insight on the evolution of the electronic structure from zero dimensional QDs to one-dimensional quantum wires.

The combination of scanning tunneling spectroscopy with PLE has proven to be a powerful approach to decipher the level structure of nanocrystal QDs Millo et al. (2001); Banin et al. (1999). While in the optical spectra, allowed valence band (VB) to conduction band (CB) transitions are detected Norris and Bawendi (1996); Banin et al. (1998), in tunneling spectroscopy the CB and VB states can be separately probed yielding complimentary information on the level structure Banin et al. (1999); Bakkers and Vanmaekelbergh (2000); Bakkers et al. (2001); Katz et al. (2001). These data provide an important benchmark for testing theoretical models developed for describing the level structure in strongly quantum-confined nanostructures. This was demonstrated for spherical nanocrystal QDs where multi-band effective mass approaches Norris and Bawendi (1996); Banin et al. (1998); Ekimov et al. (1993) and atomistic pseudo-potential theory Fu and Zunger (1998) were both used to describe measured size-dependent level structure. In a recent work by Hu et al. Hu et al. (2002), the pseudo-potential approach was also applied for the calculation of energy levels in rods of small diameter and aspect ratios up to six.

The CdSe system studied here is highly developed synthetically, to a level where it serves as a model system for colloidal semiconductor nanostructures. Fig. 1 presents transmission electron microscopy (TEM) images of typical samples of CdSe quantum rods grown using the well developed methods of colloidal nanocrystal synthesis Peng et al. (2000); Manna et al. (2000); Peng and Peng (2001). For the study of length and diameter dependence we prepared six different samples with dimensions that can be divided into two groups- three samples with radii $`1.8`$ nm and lengths ranging from 11 to 31 nm, and three samples, with radii $`3.2`$ nm and lengths ranging from 11 to 60 nm. For the low temperature optical experiments, optically clear free-standing polymer films containing the rods were prepared. The cooled samples showed absorption spectra with several transitions and had a distinct photoluminescence (PL) peak assigned to band-gap emission. The positions of the absorption onset and PL peak red shifted with increasing diameter and showed no significant variation with rod length \[Fig. 4(a)\] Li et al. (2001).

The PLE spectra presented in Fig. 1, measured on four samples, were obtained by opening a narrow detection window at the blue edge of the inhomogeneously broadened PL peak Norris and Bawendi (1996); Banin et al. (1998). The spectra are less structured as compared with those measured on QDs due to the increased sources of inhomogeneous broadening in rods and due to the intrinsically less discrete level structure. A striking feature is the nearly identical level structure observed for the two QR samples of small diameter (traces a,b), which differs significantly from the spectra measured on the thicker rods (traces c,d), although in each group the aspect ratios vary from 3 to about 10. This shows clearly that not only the band-gap of the rods depends mainly on the diameter, but also the excited optical transitions.

The optical measurements were accompanied by tunneling measurements performed at $`4.2`$ $`K`$ on $`single`$ QRs deposited on graphite, using a cryogenic scanning tunneling microscope (STM) Katz et al. (2001). A topographic image of a nanorod 25 nm in length and 2 nm in radius is presented in Fig. 2(a). After identifying such an isolated QR, tunneling $`dI/dV`$ vs. $`V`$ spectra were obtained in a double barrier tunnel junction (DBTJ) configuration Banin et al. (1999); Katz et al. (2001); Bakkers and Vanmaekelbergh (2000); Bakkers et al. (2001), as shown schematically in Fig. 2(b). The spectra were measured with the tip retracted from the QR to a distance where the bias voltage is dropped largely across the tip-QR junction, and charging effects are avoided. As a result, CB (VB) states appear at positive (negative) sample bias, and the peak separations are close to the real QR level spacings. Hence, the $`dI/dV`$ spectra yield direct information on the QR level spectrum Katz et al. (2001); Bakkers and Vanmaekelbergh (2000); Bakkers et al. (2001).

The tunneling spectra in Fig. 2(c) further demonstrate, in accordance with the PLE data, that the QR level structure depends primarily on the diameter of the QRs, not on their length. Most significantly, the region of suppressed tunneling conductance (null density of states) around zero bias, associated with the quasi-particle energy gap, is red shifted upon QR thickening, from $`2.4`$ eV in the upper two curves ($`r2`$ nm) to $`2.2`$ eV in the lower curves ($`r2.7`$ mn). This trend was not observed as clearly for the spacing between the CB ground state (CB1) and the first excited state (CB2), as evident also from Fig. 4(b). Level CB3 appearing in the lower trace was observed in $`50`$% of the measurements, where the current did not reach the saturation limit before resolving this level. The VB is considerably more dense and complex (see below) and its level structure was not reliably resolved in our tunneling spectra; we thus only denote the first (ground state) peak as VB1 \[Fig. 2(c)\]. More insight into the QR level structure, including the excited VB levels, is gained by correlating the tunneling spectra with allowed optical VB to CB transitions and comparing these data with model level-structure calculations.

To calculate the electronic structure of the rods we employ the formalism developed by Sercel and Vahala for quantum wires Sercel and Vahala (1990). Similar basis functions were used for the envelope wavefunctions, with quantum numbers reflecting the cylindrical symmetry of the QR: $`n`$, the principal number (associated with the number of nodes along the radius), and $`F_z`$ and $`L_z`$, the projections of the total and the envelope angular momenta along the z (symmetry) axis, respectively. Here, $`F_z`$ = $`J_z`$ \+ $`L_z`$, where $`J_z`$ is the projection of band-edge Bloch angular momentum along the z axis. The finite length $`L`$ of the QRs is expressed by introducing an additional quantum number, denoted by $`m`$, and a related factor, $`exp(i\pi mz/L)`$, to the envelope wavefunctions, replacing the continuum wavevector $`K_z`$ in Ref. Sercel and Vahala, 1990.

We first consider the CB level structure. Here, due to the large band-gap of CdSe (1.84 eV), the uncoupled one-band effective-mass model with quantum numbers ($`n,L_z,m`$) is adequate for rods of large enough aspect ratio Sercel and Vahala (1990); Hu et al. (2002). Figure 3(a) presents CB energy levels calculated using this model for a particle in a finite cylindrical potential-well ($`U=5`$ eV), as a function of radius and length (inset). The CB states with $`L_z=0,1,`$ and 2 are denoted CB1, CB2 and CB3, respectively. This figure clearly depicts the strong dependence of the energy levels on the QR radius, due to the strong confinement associated with this ’narrow dimension’, as opposed to the near lack of length dependence shown in the upper inset. The length dependence becomes significant only for small aspect ratio (around 2), or for high $`m`$ values, effectively shortening the QR, since the number of nodes along the z axis is $`m1`$, in agreement with Ref. Hu et al., 2002.

The VB states were calculated using the four-band effective-mass model that couples heavy and light holes Sercel and Vahala (1990), where the quantum numbers are now ($`n,F_z,m`$). Results of this calculation, performed for an infinite cylindrical potential-well, are presented in Fig. 3(b), where states with $`F_z=1/2,3/2,`$ and 5/2 (and $`m=1`$) are shown. Again, the near lack of dependence of the level structure on the rod length is manifested by the inset. It is evident that the VB is significantly more dense as compared to the CB. Yet, the VB levels near the band-edge can be divided into three groups: VB1, denoting the two upper levels, VB2, the three subsequent ones and VB3, the next group. Each of these groups is composed of various $`F_z`$ levels, where, in turn, each $`F_z`$ level has contributions from both heavy and light holes ($`J_z=3/2`$ and 1/2, respectively), and thus from wavefunctions of different $`L_z`$ values. This VB mixing leads to very weak selection rules for the optical transitions.

In Fig. 4 we compare the measured optical transitions and tunneling spectra with our theoretical calculations. The gap extracted from the optical measurements (solid circles) and the gap identified in the tunneling data (empty circles), are plotted in Fig. 4(a) along with the calculated energy gap, VB1-CB1. In order to compare with the quasi-particle (tunneling) gap, the measured excitonic (optical) gap was corrected for the electron-hole Coulomb interaction. As a first approximation, we modified the expression given in Ref. Norris and Bawendi, 1996 for QDs, $`1.8e^2/kr`$ (where $`k`$ is the dielectric constant), to take into account the fact that in our QRs strong confinement holds only for the radial dimension. We thus replace the radius, $`r`$, by $`(r^2a_o)^{1/3}`$, where $`a_o`$ is the Bohr radius, 5.7 nm. A relatively good agreement between the calculated and measured energy gaps is found, for both tunneling and optical experiments. However, the tunneling gaps are higher in energy and show a stronger dependence on radius as compared to the optical data, as was previously observed for QDs Banin et al. (1999); Norris and Bawendi (1996); Banin et al. (1998). This can be ascribed to the non-vanishing voltage drop on the QR-substrate junction Katz et al. (2001); Bakkers and Vanmaekelbergh (2000); Bakkers et al. (2001) and on the uncertainty involved in the determination of the QR radii Banin et al. (1999). Concerning the optical gap, the transition from both $`F_z=1/2`$ and 3/2 VB1 states are allowed. In the absorption they are too close to be resolved, but we note that the state with $`F_z=1/2`$ is the top-most state, from which theory predicts a transition polarized along the symmetry axis Sercel and Vahala (1991), consistent with PL polarization measurements for QRs Hu et al. (2001).

Turning now to the excited levels, Fig. 4(b) shows spacings of the PLE transitions I and II with respect to the band-gap transition, together with levels CB2 and CB3 measured with respect to CB1, detected by tunneling. The data are presented as a function of the band-gap measured in each experiment, thus eliminating the possible problem of QR radius estimation mentioned above. Most interesting is the good agreement of the spacing between PLE transition II and the band-gap transition (solid squares) with the tunneling data for the CB2-CB1 level separation (open squares), allowing a clear assignment of optical transition II. Both the optical and tunneling data correlate with the theoretical curve for the CB2-CB1 spacing. It is also clear now that PLE transition I (solid stars) can only take place between an excited VB state and the CB ground state. However, due to the intricate VB level structure and the broad PLE features, it cannot be unambiguously assigned, and relatively good agreement is found for both calculated VB2-VB1 and VB3-VB1 level spacing (dashed lines). Transition III was resolved only for the three samples with larger diameters, and thus cannot be assigned. Considering again the tunneling data, Fig. 4(b) shows a weaker size (energy-gap) dependence for CB2 and CB3 as compared to theory. This may be partly due to our oversimplified model, where we assume a sharp (square-well) and somewhat large (5 eV) confining potential. These approximations should affect mainly the higher levels.

The contribution of high $`m`$ states was not clearly resolved in our spectra. In the PLE data, such a contribution possibly manifests itself in the broadening of the peaks, as compared with QDs. As for the tunneling spectra, due to the small wavector component parallel to the surface associated with the tunneling electron, coupling to higher $`m`$ values may be suppressed. We note, however, that our spectra often exhibit a broad background on top of which the peaks are superimposed (e.g., the middle curves in Fig. 2).

In conclusion, optical and tunneling spectroscopic data for CdSe QRs show that the level structure is dominated by radius rather than by the length. This behavior, manifested also in the theoretical calculations, reveals the quasi one-dimensional nature of the QR even for aspect ratios as small as three. This allows one to select a desired QR length and tune its optical or electrical properties at will, using the diameter. This ability is of significant importance for future nanotechnology applications of QRs Kazes et al. (2002); Huynh et al. (2002).

###### Acknowledgements.

This work was supported in parts by the US-Israel Binational Foundation, the DIP foundation and by the Israel Science Foundation.


text:
# Submitted Landau damping of partially incoherent Langmuir waves

## Abstract

It is shown that partial incoherence, in the form of stochastic phase noise, of a Langmuir wave in an unmagnetized plasma gives rise to a Landau-type damping. Starting from the Zakharov equations, which describe the nonlinear interaction between Langmuir and ion-acoustic waves, a kinetic equation is derived for the plasmons by introducing the Wigner-Moyal transform of the complex Langmuir wave field. This equation is then used to analyze the stability properties of small perturbations on a stationary solution consisting of a constant amplitude wave with stochastic phase noise. The concomitant dispersion relation exhibits the phenomenon of Landau-like damping. However, this damping differs from the classical Landau damping in which a Langmuir wave, interacting with the plasma electrons, loses energy. In the present process, the damping is non-dissipative and is caused by the resonant interaction between an instantaneously-produced disturbance, due to the parametric interactions, and a partially incoherent Langmuir wave, which can be considered as a quasi-particle composed of an ensemble of partially incoherent plasmons.

PACS number(s): 52.35.Mw, 52.35.Fp, 52.25.Dg, 03.40.Kf

Keywords: quasi-particles, Langmuir envelopes, Landau damping, modulational instability, nonlinear Schroedinger equation, Zakharov equations.

More than fifty years ago, Landau  discovered a damping of electron plasma waves whose linear dynamics is governed by the Vlasov-Poisson system of equations. This damping was caused by resonant wave-electron interactions. The Landau theory clearly shows that the decay rate of the wave energy due to the interaction with the resonant electrons of the system (satisfying $`p=\omega /k`$, where $`p`$ denotes the electron velocity and $`\omega `$ and $`k`$ are the frequency and wave number of the electron plasma wave ) is proportional to the first derivative of the equilibrium distribution function $`\rho _0(p)`$ of the electrons. Typically the shape of $`\rho _0(p)`$ is such that $`d\rho _0(p=\omega /k)/dp<0`$, which implies that there are more particles with $`p<\omega /k`$ (which gain energy from the wave) than with $`p>\omega /k`$ (which give energy to the wave). Consequently, the net result is that the wave is damped. However, finite amplitude Langmuir waves in plasmas can be created when some free energy sources, such as electron and laser beams, are available in the system. In the past, several authors - have considered the nonlinear coupling between high-frequency Langmuir and low-frequency ion-acoustic waves. Physically, this coupling occurs because a large amplitude plasma wave, interacting with small density ripples, produces a current which becomes the source for an envelope of Langmuir waves. The low-frequency ponderomotive force of the latter reinforces the density modulation in such a way as to produce a space-charge electric field, which varies on a time scale longer than the electron plasma period. Thus, within the slowly-varying envelope approximation, the dynamics of the Langmuir wavepacket is governed by a nonlinear Schrödinger equation (NLSE), while the non-resonant density ripples follow an ion-acoustic wave equation driven by a ponderomotive force created by the Langmuir wave). In his classic paper, Zakharov  demonstrated that these nonlinearly coupled equations admit a class of modulationally unstable stationary solutions and also exhibit the phenomenon of wave collapse.
In this paper, we present a Landau damping of partially coherent Langmuir waves whose dynamics is governed by the Zakharov equations. In fact, we show that the latter can be reduced to a kinetic-like equation for a Wigner-Moyal-like transform for the amplitude of the electric field associated with the Langmuir wave. Based on this kinetic equation we consider the dynamics of small perturbations on a partially coherent background solution consisting of a constant amplitude and a stochastically varying phase. The concomitant dispersion relation depends strongly on the spectral energy density of the Langmuir wave and a Landau-type damping of the perturbations is found due to the resonant interaction between the Langmuir quasi-particle and the instantaneously created ripples.

As it is well known, the nonlinear propagation of a one–dimensional Langmuir wave–packet in an unmagnetized plasma with unperturbed density $`n_0`$, is governed by the Zakharov equations

$$i\frac{E}{t}+\frac{3v_{te}^2}{2\omega _{pe}}\frac{^2E}{x^2}\omega _{pe}\frac{\delta n}{2n_0}E=0,$$
(1)
$$\left(\frac{^2}{t^2}c_s^2\frac{^2}{x^2}\right)\frac{\delta n}{2n_0}=\frac{c_s^2}{2}\frac{^2}{x^2}\left(\frac{|E|^2|E_0|^2}{4\pi n_0T_e}\right),$$
(2)
where $`E`$ is the slowly varying complex electric field amplitude, $`\delta n`$ is the density fluctuation, $`v_{te}`$ is the electron thermal velocity, $`\omega _{pe}`$ is the electron plasma frequency, $`T_e`$ is the electron temperature, $`c_s`$ is the sound speed, and $`E_0`$ is the unperturbed electric field amplitude. Moreover, $`x`$ and $`t`$ play the role of configurational space coordinate (longitudinal coordinate with respect to the wavepacket centre) and time, respectively. We will show that the nonlinear propagation of a Langmuir wavepacket, governed by the above Zakharov system, can be suitably described in phase space by a kinetic–like equation within a framework similar to the one of the Vlasov-Poisson system. A similar approach has recently been used to study the modulational instability of a large-amplitude electromagnetic wavepacket as well as the coherent instability analysis of intense charged–particle beam dynamics . Within this framework, we show that a phenomenon similar to the classical Landau damping occurs for partially coherent plasma waves. In particular, this phenomenon tends to suppress the modulational instability. To this end, the above system of equations (1) and (2) can be cast in the following form

$$i\alpha \frac{\mathrm{\Psi }}{s}+\frac{\alpha ^2}{2}\frac{^2\mathrm{\Psi }}{x^2}U\mathrm{\Psi }=0,$$
(3)
$$\left(\frac{^2}{s^2}\mu ^2\frac{^2}{x^2}\right)U=\frac{^2}{x^2}\left(|\mathrm{\Psi }|^2|\mathrm{\Psi }_0|^2\right),$$
(4)
where the angle brackets account for the statistical ensemble average ála Klimomtovich  due to the partial incoherence of the waves, $`s=\sqrt{3}\lambda _{De}\omega _{pe}t`$ ($`\lambda _{De}v_{te}/\omega _{pe}`$ is the electron Debye length), $`\alpha =\sqrt{3}\lambda _{De}`$, $`U=\delta n/2n_0`$, $`\mu =\sqrt{m_e/3m_i}`$ ($`m_e`$ and $`m_i`$ are the electron and the ion masses, respectively), and

$$\mathrm{\Psi }=\frac{\mu E}{\sqrt{8\pi n_0T_e}},$$
(5)
$$\mathrm{\Psi }_0=\frac{\mu E_0}{\sqrt{8\pi n_0T_e}}.$$
(6)
Note that (6) can be written as

$$\rho _0|\mathrm{\Psi }_0|^2=\mu ^2\frac{ϵ}{ϵ_T},$$
(7)
where $`ϵE_0^2/8\pi `$ and $`ϵ_Tn_0T_e`$ are the electric energy density and the electron thermal energy density, respectively. According to the theory modelled by the Zakharov system (3) and (4), $`ϵ/ϵ_T>>1`$ (i.e, the ponderomotive effect is larger than the thermal pressure effect). Note also that Eq. (4) implies that $`U`$ is a functional of $`|\mathrm{\Psi }|^2`$, namely $`U=U\left[|\mathrm{\Psi }|^2\right]`$. Consequently, once (4) is combined with (3), the latter becomes a nonlinear Schrödinger equation (NLSE) with the nonlinear term $`U\left[|\mathrm{\Psi }|^2\right]`$. This implies that, as in Quantum Mechanics, one can introduce the density matrix  as a sort of ”two-points correlation function” and then the Wigner-Moyal transform  which is solution of the von Neumann-Weyl equation . In case of surface gravity waves, whose dynamics and instability is described in deep water by a nonlinear Schrödinger equation, such a kind of approach has been given in the pioneering works by Alber , Crawford et al.  and Janssen . More recently, a similar approach has been developed for the propagation of electromagnetic wavepackets in nonlinear media where an appropriate kinetic equation is able to show a random version of the modulational instability . In particular, this approach can be also seen as an application of the Klimontovich statistical average method . The transition to the phase-space $`xp`$ ($`pdx/ds`$ is the conjugate variable of $`x`$) allows us to write a von Neumann equation  for $`w(x,p,s)`$. In fact, although, we are dealing with a classical system, by following Alber , one can introduce the correlation function (whose corresponding meaning in Quantum Mechanics is nothing but the density matrix for mixed states)

$$\rho (x,x^{},s)=\mathrm{\Psi }^{}(x,s)\mathrm{\Psi }(x^{},s),$$
(8)
where the statistical ensemble average takes into account the incoherency of the waves. The technique was successfully introduced in the field of statistical quantum mechanics to describe the dynamics of a system in the classical space language -. Thus one can define the following Wigner-Weyl transform

$$w(x,p,s)=\frac{1}{2\pi \alpha }_{\mathrm{}}^{\mathrm{}}\rho (x+\frac{y}{2},x\frac{y}{2},s)\mathrm{exp}\left(i\frac{py}{\alpha }\right)𝑑y,$$
(9)
which allows us to transit from the configuration space to the phase space. In fact, $`w(x,p,s)`$ is governed by the following kinetic-like equation (von Neumann-like equation  )

$$\frac{w}{s}+p\frac{w}{x}\underset{n=0}{\overset{\mathrm{}}{}}\frac{(1)^n}{\left(2n+1\right)!}\left(\frac{\alpha }{2}\right)^{2n}\frac{^{2n+1}U}{x^{2n+1}}\frac{^{2n+1}w}{p^{2n+1}}=0.$$
(10)
On the other hand, since (9) implies that

$$|\mathrm{\Psi }|^2=_{\mathrm{}}^{\mathrm{}}w(x,p,s)𝑑p,$$
(11)
the functional $`U\left[|\mathrm{\Psi }|^2\right]`$ can be expressed as a functional of $`w`$, namely $`U=U\left[w\right]`$. Consequently, (4) can be written as

$$\left(\frac{^2}{s^2}\mu ^2\frac{^2}{x^2}\right)U=\frac{^2}{x^2}\left(_{\mathrm{}}^{\mathrm{}}w𝑑p_{\mathrm{}}^{\mathrm{}}w_0𝑑p\right),$$
(12)
where $`w_0`$ is the Wigner-Moyal transform of $`\mathrm{\Psi }_0`$ ($`\rho _0=|\mathrm{\Psi }_0|^2=_{\mathrm{}}^{\mathrm{}}w_0𝑑p`$). The above scheme allows us to carry out a stability analysis of the stationary Langmuir wave by considering the linear dispersion relation of small amplitude perturbations. We start from the equilibrium state: $`w=w_0(p)`$, $`U=U_0=0`$, and perturb the system according to

$$w(x,p,s)=w_0(p)+w_1(x,p,s),$$
(13)
$$U(x,s)=U_0+U_1(x,s)=U_1(x,s),$$
(14)
where $`w_1(x,p,s)`$ and $`U_1(x,s)`$ are first-order quantities. Consequently, (10) and (12) can be linearized as follows

$$\frac{w_1}{s}+p\frac{w_1}{x}=\underset{n=0}{\overset{\mathrm{}}{}}\frac{(1)^n}{\left(2n+1\right)!}\left(\frac{\alpha }{2}\right)^{2n}\frac{^{2n+1}U_1}{x^{2n+1}}w_0^{(2n+1)},$$
(15)
$$\left(\frac{^2}{s^2}\mu ^2\frac{^2}{x^2}\right)U_1=\frac{^2}{x^2}\left(_{\mathrm{}}^{\mathrm{}}w_1𝑑p\right),$$
(16)
where $`w_0^{(2n+1)}d^{2n+1}\rho _0/dp^{2n+1}`$. By introducing the Fourier transform of $`U_1(x,s)`$ and $`w_1(x,p,s)`$, i.e.

$$U_1(x,s)=_{\mathrm{}}^{\mathrm{}}𝑑k_{\mathrm{}}^{\mathrm{}}𝑑\omega \stackrel{~}{U_1}(k,\omega )\mathrm{exp}\left(ikxi\omega s\right),$$
(17)
$$w_1(x,p,s)=_{\mathrm{}}^{\mathrm{}}𝑑k_{\mathrm{}}^{\mathrm{}}𝑑\omega \stackrel{~}{w_1}(k,p,\omega )\mathrm{exp}\left(ikxi\omega s\right),$$
(18)
and substituting (17) and (18) into Eq.s (15) and (16), we readily obtain the following dispersion relation

$$\omega ^2\mu ^2k^2=k^2_{\mathrm{}}^{\mathrm{}}\frac{w_0\left(p+\alpha k/2\right)w_0\left(p\alpha k/2\right)}{\alpha k}\frac{dp}{p\omega /k}.$$
(19)
We consider a Lorenzian spectrum of the form:

$$w_0(p)=\frac{\rho _0}{\pi }\frac{\mathrm{\Delta }}{\mathrm{\Delta }^2+p^2},$$
(20)
where $`\mathrm{\Delta }`$ is the width of the spectrum. Dispersion relation (19) becomes:

$$\left(\omega ^2\mu ^2k^2\right)\left(\omega ^2\alpha ^2k^4/4+2ik\mathrm{\Delta }\omega k^2\mathrm{\Delta }^2\right)=\rho _0k^4,$$
(21)
As the simplest case, we take the limit of $`\mathrm{\Delta }0`$. This case corresponds to a coherent monochromatic wave-train and the dispersion relation reduces to

$$\left(\omega ^2\mu ^2k^2\right)\left(\omega ^2\alpha ^2k^4/4\right)=\rho _0k^4,$$
(22)
which predicts the well known result of a purely growing instability .

In the following, we concentrate on the limit of small $`\alpha k<<1`$ and $`\omega /k>>1`$. We will show that during its nonlinear evolution the Langmuir wavepacket can exhibit a phenomenon, similar to the Landau damping , between parametrically-driven non-resonant density ripples and Langmuir quasi-particles. It is easy to verify that the condition $`\alpha k<<1`$ (i.e., $`\lambda _{De}k<<1`$) implies that

$$\frac{w_0\left(p+\alpha k/2\right)w_0\left(p\alpha k/2\right)}{\alpha k}dw_0/dpw_0^{^{}}.$$
(23)
Consequently, (19) becomes

$$\omega ^2=k^2\frac{w_0^{^{}}}{p\omega /k}𝑑p,$$
(24)
where, using the assumption $`\omega /k>>1`$, we can neglect the term involving $`\mu ^2`$ because $`\mu <<1`$.
Note that (24) is formally identical to the linear dispersion relation that holds for warm plasma waves  or for non-monochromatic charged-particle bunches in circular accelerating machines for an arbitrary complex coupling impedance . Consequently, the dispersion relation (24) allows us to predict a sort of Landau damping , as described in plasma physics for linear plasma waves as well as in charged-particle beam physics. A similar effect, which has been called quantum-like Landau damping (QLLD), has recently been predicted for the nonlinear propagation of a large-amplitude electromagnetic waves in nonlinear media, such as optical fibers and plasmas as well as for the quantum-like description of charged-particle beams in an accelerating machine . In order to investigate this phenomenon also for large-amplitude Langmuir waves, we apply the standard procedure  used for the Landau damping of a linear plasma wave, which involves integrating the above dispersion relation in the complex plane, evaluating both the Cauchy principle value (PV) and the semi-residue in $`p=\omega /k`$. The result is

$$\omega ^2=k^2_{PV}\frac{w_0^{^{}}}{p\omega /k}𝑑p+i\pi k^2w_0^{^{}}(\beta ),$$
(25)
For a symmetric sufficiently smooth stationary background distribution, $`w_0(p)`$, we have

$$\omega ^2=\left(1+\frac{3\mathrm{\Delta }^2}{\rho _0}\right)k^2+i\pi k^2w_0^{^{}}(\beta ),$$
(26)
where $`\mathrm{\Delta }\left(p^2\right)^{1/2}`$ denotes the r.m.s. width of $`w_0(p)`$, and $`\beta \rho _0^{1/4}\left(1+3\mathrm{\Delta }^2/4\rho _0^2\right)`$. Note that, according to the above hypothesis, we have $`\beta >>1`$ and $`\mathrm{\Delta }^2/\rho _0<<1`$. This implies that $`\rho _0=|\mathrm{\Psi }_0|^2>>1`$, and from (7) it follows that $`ϵ/ϵ_T>>1/\mu ^2`$. Eq. (26) clearly shows that the damping rate is proportional to the derivative of the Wigner distribution $`w_0`$. This is formally similar to the expression for the Landau damping of a linear plasma wave in a warm unmagnetized plasma. In fact, writing $`\omega `$ in the complex form $`\omega =\omega _R+i\omega _I`$, substituting it in (25) and then separating the real and imaginary parts, we obtain

$$\omega _R^2\omega _I^2=k^2_{PV}\frac{w_0^{^{}}}{p\omega /k}𝑑p,$$
(27)
and

$$\omega _I=\frac{\pi }{2}\frac{k^2}{\omega _R}w_0^{^{}}(\beta ).$$
(28)
The Landau damping rate of the Langmuir wave in the appropriate unities is $`\gamma =\sqrt{3}\omega _{pe}\lambda _{De}\omega _I`$. For example, for a Gaussian wavepacket spectrum, i.e. $`w_0(p)=\left(\rho _0/\sqrt{2\pi \mathrm{\Delta }^2}\right)\mathrm{exp}\left(p^2/2\mathrm{\Delta }^2\right)`$, one obtains

$$\gamma =\left(\frac{3\pi }{8}\right)^{1/2}\frac{\mu ^{5/2}}{\mathrm{\Delta }^3}\left(\frac{ϵ}{ϵ_T}\right)^{5/4}\omega _{pe}\lambda _{De}k\mathrm{exp}\left[\frac{\mu }{2\mathrm{\Delta }^2}\left(\frac{ϵ}{ϵ_T}\right)^{1/2}\right],$$
(29)
where higher-order terms have been neglected. The present damping mechanism differs from the standard Landau damping in that here $`w_0(p)`$ does not represent the equilibrium velocity distribution of the plasma electrons, but it can be considered as the ”kinetic” distribution of all Fourier components of the partially incoherent large-amplitude Langmuir wave in the warm plasma. In terms of plasmons, we realize that $`w_0(p)`$ represents the distribution of the partially incoherent plasmons that are distributed in $`p`$-space (i.e. in $`k`$-space) with a finite ”temperature”. Consequently, the QLLD described here is due to the partial incoherence of the wave which corresponds to a finite-width Wigner distribution of the plasmons.

In this paper, we have carried out an analysis for a Landau-type damping in plasmas containing partially incoherent Langmuir waves. Starting from the Zakharov system of equations and employing the Wigner-Moyal quasi-distribution function for the complex Langmuir wave electric field amplitude, it has been shown that the Zakharov system of equations can be converted into a pair of coupled evolution equations in phase space. For a monochromatic and coherent wave ($`w_0(p)=\rho _0\delta (p)`$), we have obtained a dispersion relation which exhibits a purely growing instability . Furthermore, for a broad-band Langmuir wave spectrum, corresponding to a partially incoherent constant amplitude wave, we have investigated the properties of small perturbations for $`\alpha k<<1`$ (i.e., $`\lambda _{De}k<<1`$), and very large dimensionless phase velocity (i.e., $`\omega /k>>1`$). Our dispersion relation for this case exhibits a sort of weak Landau damping, very similar to the one predicted for the Vlasov-Poisson system for linear plasma waves. The physical origin of this phenomenon is attributed to the ”non-monochromatic” behaviour of the Wigner spectrum of the Langmuir wave, which can be considered as an ”ensemble of incoherent plasmons”. Actually, a broad-band spectrum of incoherent plasmons, forming the ”wave packet”, interact individually with the self-excited density ripples. Similar to the standard Landau damping, where the electrons interact individually with a linear plasma wave and statistically produce a net transfer of energy from the wave to the particles, the plasmons interact individually with the ripples and produce a change of energy which is more significant around the resonance, which is determined by the condition $`p=\omega /k`$. Finally, we observe that since the system can be modulationally unstable , this phenomenon acts in competition with the modulational instability.