format: text

text:
# Abstract

## Abstract

We show that and how the Coulomb potential $`V(x)=Ze^2/x`$ can be regularized and solved exactly at the imaginary coupling $`Ze^2`$. The new spectrum of energies is real and bounded as expected, but its explicit form proves totally different from the usual real-coupling case.

## 1 Introduction

Quantum mechanics often works with the exactly solvable simplified models. For the precise fits of data or for some more subtle quantitative analyses, unfortunately, the number of solvable models is too limited. In $`D`$ dimensions, the only really useful and easily tractable interactions are harmonic oscillators and/or the central Coulomb well $`V^{(Z)}(\stackrel{}{r})=Ze^2/|\stackrel{}{r}|`$. A new way out of this deadlock emerges within the framework of the alternative, so called “$`𝒫𝒯`$ symmetric” quantum mechanics. With its complex Hamiltonians $`H`$ breaking both the parity $`𝒫`$ and the time-reflection symmetry $`𝒯`$ and commuting only with their product $`𝒫𝒯`$, this formalism was proposed by Bessis  and by Bender et al  as a possible way towards weakening of the standard requirements of Hermiticity.

Several new exactly solvable $`𝒫𝒯`$ symmetric models have been proposed recently . This is a promising development with possible applications ranging from field theories  to supersymmetric models  and from quasi-classical methods  to perturbation theory .

Even the solvable harmonic oscillator itself acquires a richer spectrum after its consequent $`𝒫𝒯`$ symmetric regularization in $`D`$ dimensions . The detailed structure of spectrum of this prominent example does not in fact offer any really serious surprise. A manifest violation of the parity $`𝒫`$ is compensated by an emergence of the so called quasi-parity $`q=\pm 1`$ tractable as a signature of two equidistant subspectra. The new quantum number $`q`$ degenerates back to the eigenvalue of parity after a return to the standard Hermitean and one-dimensional oscillator.

No immediate surprise emerges also for the quartic anharmonic oscillator . The situation only becomes less clear after one moves towards the asymptotically vanishing models. They exhibit several counterintuitive properties and open new mathematical challenges . In particular, the popular Coulomb potential did not not even seem particularly suitable for any immediate $`𝒫𝒯`$ symmetric regularization . A psychological barrier has been created by the numerical and semiclassical studies of the general power-law forces $`V(x)(ix)^\delta `$. They may be well defined everywhere near the harmonic exponents $`\delta =2`$, $`\delta =6`$ etc . At the same time, the related analyses hinted that it is apparently difficult to move beyond the Herbst’s singularity located at $`\delta =1`$ .

In what follows, we intend to employ a slightly different strategy and try to study the Coulomb problem directly, via its well known correspondence to the harmonic oscillator. This correspondence is based on an elementary change of variables. Its background dates back to the nineteenth century mathematics and, in particular, to the work of Liouville . The Newton’s monograph  cites also Fivel  as a newer source of the idea. In the contemporary literature (cf., e.g.,  for further references) people usually speak about the Kustaanheimo - Steifel (KS) transformation . In all the implementations of this idea the parameters appearing in the Coulombic and oscillator problems are interrelated of course. Details will be mentioned below. Preliminarily, let us only warn the reader that all the KS-type mappings can also change the dimensions and angular momenta and that the energies of one problem are related to the coupling constants of the other one and vice versa. Within the ”normal” quantum mechanics, all this has already been thoroughly discussed elsewhere: In ref.  for $`D=3`$ and in ref.  for the continuous transformation between Coulomb problems and harmonic oscillators in various dimensions.

## 2 Liouvillean changes of variables

The change-of-variable approach to the Coulombic bound-state problem enables us to start directly from the harmonic oscillator potential $`W(r)=r^2`$ or, in the present less traditional context, from its $`𝒫𝒯`$ symmetric radial Schrödinger equation

$$\left[\frac{d^2}{dr^2}+\frac{l(l+1)}{r^2}+W(r)\right]\chi (r)=\epsilon ^2\chi (r)$$
(1)
of ref. , using the complex coordinate $`r=xic`$ with real $`x(\mathrm{},\mathrm{})`$ and with, say, positive $`c>0`$. This means that the integration path has been shifted down from the position where it would cross the strong centrifugal singularity. Such a regularization preserves the asymptotic decrease of the normalizable solutions. Only in the limit $`c0`$ and beyond the trivial one-dimensional case one has to omit all the so called irregular solutions  defined by their “physically unnacceptable” $`\chi (r)r^l`$ behaviour near the origin.

Within the framework of the general Liouville method the change of variables mediates a transition to the different potential $`V(t)`$. It is easy to show that once we forget about boundary conditions one merely has to demand the existence of an invertible function $`r=r(t)`$ and its few derivatives $`r^{}(t),r^{\prime \prime }(t),\mathrm{}`$. Then, the explicit correspondence between the two bound state problems may be explicitly given by the elementary formulae. From our original eq. (1) (i.e., in our case, harmonic oscillator) one obtains the new (i.e., in our case, Coulombic) Schrödinger equation

$$\left[\frac{d^2}{dt^2}+\frac{L(L+1)}{t^2}+V(t)\right]\mathrm{\Psi }(t)=E\mathrm{\Psi }(t)$$
(2)
with the new wave functions

$$\mathrm{\Psi }(t)=\chi [r(t)]/\sqrt{r^{}(t)}$$
(3)
and with the new interaction and the new energies ,

$$\frac{L(L+1)}{t^2}+V(t)E=\left[r^{}(t)\right]^2\left\{\frac{l(l+1)}{r^2(t)}+W[r(t)]\epsilon ^2\right\}+\frac{3}{4}\left[\frac{r^{\prime \prime }(t)}{r^{}(t)}\right]^2\frac{1}{2}\left[\frac{r^{\prime \prime \prime }(t)}{r^{}(t)}\right].$$

Thus, it only remains for us to re-analyse the boundary conditions.

## 3 $`𝒫𝒯`$ symmetric KS transformation

Without any serious formal difficulties let us extend the scope of the present considerations to all the singular forces $`\widehat{W}(r)=W(r)+f/r^2`$ and and/or $`\widehat{V}(t)=V(t)+F/t^2`$. Both these central forces may act in the respective $`d`$ and $`D`$ dimensions. This means that

$$l(l+1)=\left\{[j+(d3)/2][j+(d1)/2]+f\right\}$$

or rather

$$(l+1/2)^2(=\alpha ^2)=[j+(d2)/2]]^2+f$$

(with Re $`\alpha >0`$) and

$$(L+1/2)^2(=A^2)=[J+(D2)/2]]^2+F$$

(with Re $`A>0`$) where the partial waves are numbered by the respective integers $`j=0,1,\mathrm{}`$ and $`J=0,1,\mathrm{}`$.

An important simplification of our effort is provided by our knowledge of the complete harmonic oscillator solution as derived in ref. . Its two equidistant subsets of energies

$$\epsilon ^2=\epsilon _{(n,q)}^2=4n+22q\alpha ,q=\pm 1,n=0,1,\mathrm{}$$

correspond to the two families of the Laguerre-polynomial wave functions

$$\chi _{(n,q)}(r)=𝒩r^{1/2q\alpha }e^{r^2/2}L_n^{(q\alpha )}(r^2).$$

Integration path $`r=xic`$ lies in the lower half of the complex plane and does not change after the subsequent $`𝒫`$ and $`𝒯`$ transformations $`r=r(x)r`$ and $`rr^{}=r(x)`$.

In the spirit of the above-mentioned KS mapping of harmonic oscillators on Coulombic bound states we now have to define a complex variable $`t`$ as a re-scaled square of $`r(x)`$ such that the resulting path $`t(x)`$ remains $`𝒫𝒯`$ invariant. Our requirement implies that the complex plane of $`r`$ will cover twice the complex plane of $`t`$. In such an arrangement, our lower half plane of $`r`$ should cover the Riemann sheet given as a whole plane of $`t`$ which is cut upwards from the origin. In the polar representation, one has $`r\mathrm{exp}(i\phi )`$ mapped upon $`t\mathrm{exp}(2i\phi )`$ with $`\phi (0,\pi )`$. Once we introduce a suitable free parameter $`\kappa >0`$ we can put, say, $`r^2=2\kappa ^2z`$ and then rotate the $`z`$plane (which is cut, by construction, along the real and positive semi-axis) by the angle $`\pi /2`$ giving $`t=iz`$. Our final recipe

$$r^2=2i\kappa ^2t$$
(4)
maps the above-mentioned straight line $`r(x)=xic`$ upon a curve $`t(x)=u+iv`$. Its real part $`u=u(x)=xc/\kappa ^2`$ and imaginary part $`v=v(x)=(x^2c^2)/2\kappa ^2`$ form, as required, a $`𝒫𝒯`$ symmetric and upwards-oriented parabola $`v=c^2/2\kappa ^2+(\kappa ^2/2c^2)u^2`$ in complex plane. Obviously, a small asymptotic deformation of our original curve $`r(x)`$ with the modified shifts $`c=c(x)1/x^{1+\eta }`$ would transform the parabola to a pair of lines which are parallel to the cut in the asymptotic domain of $`|x|1`$.

Having achieved a $`𝒫𝒯`$ symmetry in the complex plane of $`t`$, we may move to the (trivial) insertions and conclude that all the above-mentioned harmonic oscillator bound-state solutions are in a one-to-one correspondence to the solutions of the Coulombic Schrödinger equation (2),

$$\left[\frac{d^2}{dt^2}+\frac{L(L+1)}{t^2}+i\frac{Ze^2}{t}\right]\mathrm{\Psi }(t)=E\mathrm{\Psi }(t),t=u(x)+iv(x),xIR.$$
(5)
The underlying assignment of constants is such that $`\alpha =2A`$ while $`\kappa `$ itself becomes $`n`$dependent, $`\kappa ^2=2Ze^2/\epsilon ^2=Ze^2/(2n+12qA)`$. In full detail one gets the new, Laguerre-polynomial wave functions

$$\mathrm{\Psi }_{(n,q)}(t)=t^{1/2qA}e^{i\kappa ^2t}L_n^{(2qA)}(2i\kappa ^2t)$$
(6)
and their energy spectrum specified by the elementary formula

$$E_{(n,q)}=\kappa _{(n,q)}^4=\frac{Z^2e^4}{(2n+12qA)^2}q=\pm 1,n=0,1,\mathrm{}.$$
(7)
This is our main result.

## 4 Discussion

### 4.1 Consequences of the curvature of our integration path

The latter two formulae exhibit several unusual features. The first concerns the asymptotics of the wave functions which are determined by the decreasing exponential $`\mathrm{exp}(i\kappa ^2t)`$. Its form re-confirms the correctness of the above, slightly counter-intuitive KS-dictated choice of our $`𝒫𝒯`$ symmetric integration path. Asymptotically, it encircles more or less closely the positive imaginary axis in $`t`$ plane. This clarifies the apparent paradox.

The second unexpected result is the positivity and unusual $`n`$dependence of the energies. This can be related to the choice of the KS integration path $`t(x)`$ again. In the very vicinity of the origin, one can visualize this path as a circle with radius $`\sigma `$,

$$u^2(x)+v^2(x)=\sigma ^2,|x|1.$$

From the appropriate definitions we get the formula

$$\sigma =c^2(0)/2\kappa _{(n,q)}^2+𝒪(x^2)$$

and see that this radius is $`n`$dependent and increases with the growth of this principal quantum number, $`\sigma nc^2(0)/Ze^2`$, $`n1`$. As a consequence, an “effective charge” of our $`𝒫𝒯`$ symmetric Coulomb potential appears to decrease with $`n`$ since

$$\left|\frac{iZe^2}{t}\right|\frac{Ze^2}{\sigma }+𝒪(t)=𝒪(1/n).$$

This offers a “rule-of-thumb” guide to the unusual and certainly counterintuitive $`n`$dependence of the energy levels (7). Of course, in practice, a preferred integration path will be $`n`$independent. In such a case, the $`n`$dependence re-appears in the small-$`x`$ deformation of the initial harmonic-oscillator path with $`c=c_n(x)=𝒪(1/n)`$. Such a flexible transfer of the excitation-dependence throws also a new light on the complexified KS transformation itself.

### 4.2 “Flown-away” energies and unavoided crossings

Let us return in more detail to the $`A`$dependence of our energies (7). Firstly, we notice their power-law dependence on $`n`$ and $`A`$ (with exponent = -2) as somewhat similar to the spectra in $`𝒫𝒯`$ symmetric oscillator well  (with exponent = +1) and in the Morse potential of ref.  (with exponent = +2).

In the present case, obviously, we have to distinguish between the two separate families $`E_{(n,q)}`$ with $`q=+1`$ (cf. Figure 1) and $`q=1`$ (cf. Figure 2). The latter set is, up to its sign, analogous to the ordinary Coulombic spectrum. By far not so the former one. Its energies enrich and dominate the spectrum. The $`n_{div}`$th energy “flies away” and disappears from the spectrum at $`A_{div}=n_{div}+1/2`$. Moreover, at all the positive integers and half-integers $`A`$ one encounters the unavoided level crossings. In contrast to the harmonic case, they appear at both the opposite and identical (viz., positive) quasiparities $`q`$. The former case takes place at $`A=A_{crit}=(nn^{})/2`$ while in the latter case we must fulfill the condition $`A=A_{crit}=(n+n^{}+1)/2`$. A sample of this phenomenon is given in Figure 3.

Formally the unavoided crossings generate certain identities which connect different Laguerre polynomials (cf. their sample in ref. ). In applications, these “critical” cases are not exceptional at all. For $`F=0`$ forces without a spike, the critical integer or half-integer coordinates $`A_{crit}=J1+D/2`$ correspond precisely to the physical (namely, integer) dimensions $`D`$ and partial waves $`J`$.

In the conclusion, let us not forget about many open questions. Pars pro toto, one could mention a not yet clear possibility of re-interpretation of our bound states, say, in the limit $`treal`$, i.e., beyond the mathematical and apparently natural boundaries of our present approach. Moreover, we must keep in mind that via our complexification of the coordinates we of course broke their immediate connection to any standard $`D`$dimensional problem. In this sense, our $`F=0`$ and $`F0`$ Hamiltonians differ just in an inessential way. One could even prefer the latter, Kratzer-like option as a model which is formally simpler, due to the generic absence of the puzzling unavoided crossings. At $`F0`$ the structure of the spectrum of our present Coulomb model becomes also richer and, in this sense, more interesting.

## Acknowledgements

M. Z. thanks for the hospitality of the Institute of Nuclear Research of the Hungarian Academy of Sciences, in Debrecen, and appreciates the support by the grant Nr. A 1048004 of GA AS CR. G. L. acknowledges the OTKA grant no. T031945.


text:
# References

IUHET-423

Supersymmetric left-right model and scalar potential
Biswajoy Brahmachari

Physics Department, Indiana University, Bloomington IN-47405, USA

Abstract

We construct a scalar potential of supersymmetric left-right model in the limit when supersymmetry is valid.

Left-right symmetric model is a natural extension of standard model. The symmetry breaking chain can be written as

$`G_{LR}[SU(2)_L\times SU(2)_R\times U(1)_{BL}]`$ $`{\displaystyle \genfrac{}{}{0pt}{}{\mathrm{\Delta }_R}{}}G[SU(2)_L\times U(1)_Y]`$ (1)
$`{\displaystyle \genfrac{}{}{0pt}{}{H_1H_2}{}}G_0[U(1)_{em}].`$
Symmetry breaking mechanism is normally understood in the following way. One uses symmetries of $`G_{LR}`$ to write down allowed terms of a classical scalar potential. A stable breaking of symmetry is obtained when field configuration (values of fields) is such that the potential is at its minimum at all space-time points. We know that these values of fields are the VEVs and we also know that the VEVs obey residual symmetries. Hence electric charge and color symmetries remain.

Mininum energy state of Higgs scalars need not have the symmetry properties of the gauge bosons. So the symmetry of the combined system of fermions, gauge bosons and Higgs scalars can be broken. In this paper we will first state a miminal Higgs scalar spectrum of left-right symmetric model. Then we will write them in the form of square matrices and apply Hamiltion-Cayley theorem to these matrices to see what happens. It will lead to polymonial equations (quadratic for this case) which are satisfied by these matrices. We will physically interpret these equations as minimization conditions.

Minimal Higgs choice for a model based on $`G_{LR}`$ is

$$\mathrm{\Delta }_R=(1,3,2)\mathrm{\Delta }_L=(3,1,2)\varphi =(2,2,0).$$

One can write these VEVs as matrices in a $`SU(2)_L\times SU(2)_R`$ basis where rows are $`SU(2)_R`$ multiplets and columns are $`SU(2)_L`$ multiplets. In writing so one can suppress abelian $`U(1)_{BL}`$ as well as non-abelian color degrees of freedoms. Then one gets

$$\varphi =\left(\begin{array}{cc}\kappa _1& 0\\ 0& \kappa _2\end{array}\right)\mathrm{\Delta }_L=\left(\begin{array}{cc}0& 0\\ v_L& 0\end{array}\right)\mathrm{\Delta }_R=\left(\begin{array}{cc}0& v_R\\ 0& 0\end{array}\right).$$
(2)
We can read-off that $`\varphi `$ breaks both $`SU(2)_L`$ and $`SU(2)_R`$ symmetries. However, $`\mathrm{\Delta }_R`$ breaks only the $`SU(2)_R`$ symmetry. From experiments one knows that $`SU(2)_R`$ symmetry is broken at a higher scale. This is because gauge bosons corresponding to broken $`SU(2)_R`$ symmetries are yet to be observed. Even though mass splitting between fermions within a $`SU(2)_L`$ multiplet parametrize $`SU(2)_L`$ breaking, values of parameters $`\kappa _1,\kappa _2`$ remain unknown up to magnitudes of Yukawa couplings of corresponding fermions. Similarly we can think of similar manifestations of $`SU(2)_R`$ and $`SU(2)_L\times SU(2)_R`$ symmetries in the masses of extra fermions. Thus we must find a way to study allowed values of $`\kappa _1,\kappa _2,v_L`$ and $`v_R`$ from theory. This is a motivation to further study the scalar potential.

A number of studies of minimizing scalar potential exist. Typically, one writes the most general potential using gauge symmetries and then it is minimized by taking derivatives of potential with respect to VEVs and equating them to zero. This set of equations constrain parameter space of VEVs. We consider a reverse situation. We ask whether given only matrix forms of VEVs is it possible to construct a minimal potential which leads to desired symmetry breaking for all possible values of $`\kappa _1,\kappa _2,v_L`$ and $`v_R`$ ? In other words, instead of writing the complete potential and study the ranges of $`\kappa _1,\kappa _2,v_L`$ and $`v_R`$, is it possible to obtain the unique subset of terms of the potential which allow all possible values of $`\kappa _1,\kappa _2,v_L`$ and $`v_R`$ ? This what we answer below.

The Hamilton-Cayley theorem states: Every square matrix must satisfy its own characteristic equation. That is, if

$$\mathrm{det}(𝐀\lambda 𝐈)=c_n\lambda ^n+c_{n1}\lambda ^{n1}+\mathrm{}+c_2\lambda ^2+c_1\lambda +c_0$$
(3)
then

$$c_n𝐀^n+c_{n1}𝐀^{n1}+\mathrm{}+c_2𝐀^2+c_1𝐀+c_0=0.$$
(4)
Using equations (2) (3) and (4) we get

$$\varphi ^2(\kappa _1+\kappa _2)\varphi +\kappa _1\kappa _2=0\mathrm{\Delta }_L^2=0\mathrm{\Delta }_R^2=0.$$
(5)
These equations are satisfied for any value of $`\kappa _1,\kappa _2,v_L`$ and $`v_R`$.

Suppersymmtry allows only trilinear terms in the superpotential. Thus minimization conditions are at most quadratic. It is a nice coincidence that $`2\times 2`$ matrices lead to quadratic characteristic equations. The term $`\kappa _1\kappa _2`$ is present in the minimization condition. We must have a linear term in superpotential. However there is no singlet scalar in the model. So we must have

$$\mathrm{Either}:\kappa _10\mathrm{or}:\kappa _20$$
(6)
Then the superpotential reads as,

$`W=W_1+W_2`$
$`W_1={\displaystyle \frac{1}{3}}\varphi ^3{\displaystyle \frac{1}{2}}\kappa _1\varphi ^2`$
$`W_2=M_L\mathrm{\Delta }_L^2+M_R\mathrm{\Delta }_R^2`$ (7)
$`W_2`$ vanishes identically, for all values of $`v_L,v_R,M_L,M_R`$. It does not contribute to energy if the VEVs are of the form of (2). $`W_1`$ however needs to be minimized and all possible values of $`\kappa _1`$ or $`\kappa _2`$ are not allowed. We had to chose either $`\kappa _1`$ or $`\kappa _2`$ to vanish. So we have got a negative answer to our question. This is our result.

Thus we have constructed a superpotential of supersymmetric left-right model using Hamilton-Cayley theorem. We had to chose either $`\kappa _1`$ or $`\kappa _2`$ to vanish. We have chosen $`\kappa _2=0`$. This means that down sector of fermions remains massless.

This research was supported by U.S Department of Energy under the grant number DE-FG02-91ER40661.


text:
# Protostellar Collapse with Various Metallicities

## 1 Introduction

In the standard theory of cosmic structure formation, galaxies are believed to have originated from primeval density fluctuations. Many authors have been studying the galaxy formation process in this context and have made remarkable progress in understanding such processes as the gravitational growth of the density fluctuation, its decoupling from the cosmic expansion and virialization, and the subsequent radiative cooling of the baryonic gas (e.g., Padmanabhan 1993). However, processes after the cooling of baryonic gas, namely, the transformation of a pregalactic cloud into a cluster of stars, are relatively poorly known. It is inevitably necessary to investigate them in studying the galaxy formation process since the very existence of stellar components is one distinct feature of galaxies. These processes can be viewed as successive fragmentations of a cloud into fragments, or in other words, protostellar clouds, and their contraction into stars (e.g., Hayashi 1984).

Whether a cloud fragments or collapses into a single object depends on the the central flatness, that is, the axial ratio of the isodensity contour in the central region (Tsuribe & Inutsuka 1999). In particular, warm clouds do not fragment before adiabatic core formation, whereas clumps sufficiently more massive than the Jeans mass first collapse disklike and next fragment most likely into filamentary clouds (Miyama, Narita, & Hayashi 1987a, 1987b). The filamentary clouds fragment again into protostellar cloud cores after further contraction. Uehara et al. (1997) investigated the gravitational collapse of metal-free filamentary clouds using one-zone approximation and found that the minimum mass of fragments, or protostellar cores, is essentially Chandrasekhar mass, i.e., $`1M_{\mathrm{}}`$. Nakamura & Umemura (1999) confirmed Uehara et al. (1997)’s result by performing one-dimensional hydrodynamical calculations. In this paper, we will discuss the collapse of protostellar cores into stars.

The distinction between the star formation in the galaxy formation epoch and that in present-day star-forming regions mainly resides in the difference of the composition of gas that stars are made from. Metallicity, and therefore the gas to dust ratio, is lower for earlier star formation. In present-day star forming regions, protostellar clouds remain roughly 10K in a wide range of densities owing to efficient dust emission (e.g., Hayashi & Nakano 1965) in the course of contraction.

On the other hand, many authors have studied the thermal and chemical evolution of collapsing primordial clouds with one-zone approximation (e.g., Matsuda, Sato, & Takeda 1969; Yoneyama 1972; Calberg 1981; Palla, Salpeter, & Stahler 1983; Izotov & Kolesnik 1984; Lahav 1986; Puy & Signore 1997) or by hydrodynamical calculations (Matsuda et al. 1969; Villere & Bodenheimer 1987; Haiman, Thoul, & Loeb 1996; Omukai & Nishi 1998; Oliveira et al. 1998). The collapse of metal-free protostellar clouds, which is relevant for the first star formation in the universe, is induced by radiative cooling owing to molecular hydrogen lines, and their temperatures are about 1000K.

Above are two extremes in composition of protostellar gas, namely, the former for the high-metallicity end and the latter for the low-metallicity end.

Then, how does the collapse of slightly metal-polluted clouds proceed? Even in the galaxy formation epoch, significant metal enrichment can occur since the lifetime of massive stars ($`10^6`$ yr) is shorter than the free-fall time of pregalactic clouds ($`10^8`$ yr). In fact, both Ly$`\alpha `$ forest clouds and Population II stars are slightly polluted by heavy elements. Therefore, it is important to study star formation from slightly metal-polluted gas. In this paper, we aim to fill the gap between the protostellar collapse of primordial and present-day interstellar gas.

This topic has been investigated in the pioneering works of Low & Lynden-Bell (1976) and Yoshii & Sabano (1980). However, the following points remain to be fulfilled; in their works, (i) chemical reactions had not been solved, (ii) molecular coolants, which become important at high densities, had not been considered, and (iii) the evolution of clouds as dense as the stellar density ($`>10^{22}\mathrm{cm}^3`$) had not been studied. With these points in mind, we investigate in this paper the thermal and chemical evolution of collapsing spherical clouds as a function of metallicity.

The outline of this paper is as follows. In §2, we describe the method of our calculations. In §3, results of our calculations are presented. We summarize our work in §4.

## 2 Model

### 2.1 Basic Equations

We consider a spherical cloud with mean metallicity $`Z`$. The helium concentration is assumed to be $`y_{\mathrm{He}}=0.0972`$ for all clouds.111The concentration of element X is defined by
$$y_X=n_X/n,$$
(1)  where $`n`$ and $`n_X`$ are the number densities of hydrogen nuclei and nuclei of element X. Similarly, we write for each atomic, molecular, or ionic species
$$y(x)=n(x)/n,$$
(2)  where $`n(x)`$ is the number density of species $`x`$. Note $`y(\mathrm{H}_2)=1/2`$ for fully molecular gas.  We assume that the dust-to-gas ratio is proportional to the mean metallicity, or in other words, that a fixed fraction of the heavy elements in the interstellar medium condenses into dust grains. We adopt the Pollack et al. (1994) model of grains in molecular clouds (see §2.2.2). For the case of local interstellar clouds ($`Z1Z_{\mathrm{}}`$), the mass fraction of grain we used is $`0.934\times 10^2`$ in the lowest temperature regime, and the gas-phase elemental abundances are $`y_\mathrm{C}=0.927\times 10^4`$, and $`y_\mathrm{O}=3.568\times 10^4`$, respectively, which correspond to 46% of oxygen and 72% of carbon depleted onto grains. These values are reduced proportionally for lower metallicity cases. We normalize the mean metallicity (and therefore also the dust-to-gas ratio and gas-phase metallicity in our model) relative to the local interstellar values $`Z_{\mathrm{local}}`$ and denote relative metallicity by $`zZ/Z_{\mathrm{local}}`$.

We calculate the time evolution of the central density, temperature, and chemical composition of the collapsing core.

First, we assume that the dynamics is described by the free-fall relation,

$$\frac{d\rho }{dt}=\frac{\rho }{t_{\mathrm{ff}}},$$
(3)
where $`\rho `$ is the density in the central region and the free-fall time is

$$t_{\mathrm{ff}}\sqrt{\frac{3\pi }{32G\rho }}.$$
(4)
The thermal evolution is followed by solving the energy equation

$$\frac{de}{dt}=p\frac{d}{dt}(\frac{1}{\rho })^{(\mathrm{net})},$$
(5)
where

$$e=\frac{1}{\gamma _{\mathrm{ad}}1}\frac{kT}{\mu m_\mathrm{H}}$$
(6)
is the specific internal energy,

$$p=\frac{\rho kT}{\mu m_\mathrm{H}}$$
(7)
is the pressure for an ideal gas, $`\gamma _{\mathrm{ad}}`$ is the adiabatic exponent, $`T`$ is the temperature, $`\mu `$ is the mean molecular weight, $`m_\mathrm{H}`$ is the mass of hydrogen nucleus, and $`^{(\mathrm{net})}`$ is the net energy loss rate per unit mass (see §2.2).

We neglect any effect owing to rotation or magnetic fields for simplicity. In this case, the actual collapse is expected to proceed like the Penston-Larson similarity solution (Penston 1969; Larson 1969). According to this solution, the cloud consists of two parts, that is, the central core region, which has flat density distribution, and the envelope, where the density decreases outward as $`r^2`$. The size of the central flat region is roughly given by the local Jeans length $`\lambda _\mathrm{J}=\pi c_\mathrm{s}/\sqrt{G\rho }`$ in the core. Since we focus on the evolution of the central region, optical depth is estimated by that across one local Jeans length:

$$\tau _\nu =\kappa _\nu \rho \lambda _\mathrm{J}.$$
(8)
This formulation of Jeans length shielding is the same procedure as that of Low & Lynden-Bell (1976) and Silk (1977), although they introduced it as a result of successive fragmentation into Jeans mass clouds, not because of the Penston-Larson collapse.

### 2.2 Cooling/Heating processes

In addition to compressional heating, we include cooling owing to (i) atomic and molecular line radiation, (ii) energy transfer between the gas and the dust grains, which will be emitted as infrared radiation from the grains, (iii) continuous radiation from the gas, and (iv) cooling and heating associated with chemical reactions. Then

$$^{(\mathrm{net})}=_{\mathrm{line}}+_{\mathrm{gr}}+_{\mathrm{cont}}+_{\mathrm{chem}},$$
(9)
where the terms in the right-hand side of the equation correspond to the cooling rates owing to the processes (i)-(iv) above in the same order. Although we treat clouds with metallicity up to $`1Z_{\mathrm{}}`$, we are concerned mainly with the lower metallicity environments that are relevant to the epoch of galaxy formation. In this paper, we do not try to reproduce the present Galactic environment and neglect any external radiation (e.g., UV radiation, cosmic rays, cosmic background radiation, etc.) for simplicity. Then, the considered heating sources are the compressional work and H2 formation. As will be mentioned below, lower metallicity clouds tend to be warmer. Therefore, the compressional heating, which is proportional to the temperature of the cloud, becomes more important if the external heating rate is constant.

#### 2.2.1 Line Cooling

As line cooling agents, we include atomic fine-structure transitions of CI, CII, and OI and molecular rovibrational transitions of H2, CO, OH, and H2O:

$$_{\mathrm{line}}=_{\mathrm{CI}}+_{\mathrm{CII}}+_{\mathrm{OI}}+_{\mathrm{H}_2}+_{\mathrm{CO}}+_{\mathrm{OH}}+_{\mathrm{H}_2\mathrm{O}}.$$
(10)
The line cooling rate of species $`x`$ is given by

$$_x=\frac{1}{\rho }\underset{(ij)}{}n(x,i)A_{ij}ϵ_{ij}h\nu _{ij},$$
(11)
where $`n(x,i)`$ is the population density of species $`x`$ in level $`i`$, $`A_{ij}`$ is the spontaneous transition probability, $`ϵ_{ij}`$ is the escape probability, and $`h\nu _{ij}`$ is the energy difference between levels $`i`$ and $`j`$. The population density $`n(x,i)`$ is obtained from a solution of the equation of statistical equilibrium

$$n(x,i)\underset{ji}{\overset{n}{}}R_{ij}=\underset{ji}{\overset{n}{}}n(x,j)R_{ji},$$
(12)
where $`n`$ is the total number of lines included and

$$R_{ij}=\{\begin{array}{cc}A_{ij}ϵ_{ij}+C_{ij}\hfill & \text{for }i>j\hfill \\ C_{ij}\hfill & \text{for }i<j,\hfill \end{array}$$
(13)
ignoring the external radiation. Here $`C_{ij}`$ is the collisional rate from level $`i`$ to level $`j`$. Following Takahashi, Hollenbach, & Silk (1983), we use the escape probability for the case that the velocity is proportional to the radius in the central region:

$$ϵ_{ij}=\left(\frac{1\mathrm{e}^{\tau _{ij}}}{\tau _{ij}}\right)\mathrm{e}^{\tau _{\mathrm{cont}}}.$$
(14)
The optical depth averaged over the line and the continuum optical depth are given by

$$\tau _{ij}=\frac{A_{ij}c^3}{8\pi \nu _{ij}^3\eta _\mathrm{T}}[n(x,j)g_i/g_jn(x,i)]l_{\mathrm{sh}}/(2\mathrm{\Delta }v_\mathrm{D}),$$
(15)
and

$$\tau _{\mathrm{cont}}=(\kappa _{\mathrm{gr}}+\kappa _{\mathrm{gas}})\rho \lambda _\mathrm{J},$$
(16)
respectively, where $`\eta _\mathrm{T}`$ is the multiplicity factor, $`g_i`$ is the statistical weight of level $`i`$, $`\mathrm{\Delta }v_\mathrm{D}`$ is the velocity dispersion, and $`l_{\mathrm{sh}}`$ is the shielding length (Takahashi et al. 1983):

$$l_{\mathrm{sh}}=\mathrm{min}(\mathrm{\Delta }s_{\mathrm{th}},\lambda _\mathrm{J});\mathrm{\Delta }s_{\mathrm{th}}=2\mathrm{\Delta }v_\mathrm{D}/(\frac{dv}{dr})=6\mathrm{\Delta }v_\mathrm{D}t_{\mathrm{ff}},$$
(17)
where we have used the relation for the homogeneous collapse $`v=r/3t_{\mathrm{ff}}`$ and, in the case of the small velocity gradient, have assumed Jeans length shielding as introduced in §2. As a source of the velocity dispersion, we consider only the thermal motion of atoms and neglect microturbulent motions. Then

$$\mathrm{\Delta }v_\mathrm{D}=\sqrt{\frac{2kT}{\mu _xm_\mathrm{H}}},$$
(18)
where $`\mu _x`$ is the molecular weight of species $`x`$. The continuum Planck mean opacity of dust $`\kappa _{\mathrm{gr}}`$ and that of gas $`\kappa _{\mathrm{gas}}`$ are adopted from Pollack et al.(1994) and Lenzuni, Chernoff, & Salpeter (1991), respectively (see §2.2.2 and 2.2.3).

In general, the level population and the escape probabilities depend on each other. We then need an iterative procedure to find their consistent solution except for the case of H2. For H2 lines, we need not iterate to find the consistent set of the level populations and the escape probabilities, since in our calculations, clouds become opaque to H2 lines only at the density exceeding the critical density, and in such a case, molecules populate following a Boltzmann distribution regardless of the escape probabilities.

The related parameters of transitions are given in Hollenbach & McKee (1989) for CI, CII, and OI, in Hollenbach & McKee (1979;1989) for OH and H2O 222The OH and H2O collision cross sections and H2O multiplicity factor are given in Hollenbach & McKee (1989). We adopted other parameters from Hollenbach & McKee (1979)., and in McKee et al. (1982) for CO. We computed the population of H2 following the procedure of Hollenbach & McKee (1979) using the collision coefficient given in Hollenbach & McKee (1989). For H2, we considered the first three vibrational states ($`v=02`$) with rotational levels up to $`J=20`$ in each vibrational state. We take into account only the ground vibrational level for other molecular lines since cooling rates owing to vibrational transitions are nearly always dominated by grain cooling or H2 cooling (Hollenbach & McKee 1979).

#### 2.2.2 Gas-Grain Heat Transfer

We adopt the Pollack et al. (1994) model of grains in molecular clouds. In their model, silicates, organics, troilite, metallic iron, and water ice constitute the most abundant grain species and the dominant sources of opacity in molecular cloud cores. The Planck mean opacity due to the grains in the case of mean metallicity $`Z=Z_{\mathrm{local}}`$ is

$$\kappa _{\mathrm{gr}}=4.0\times 10^4T_{\mathrm{gr}}^2\mathrm{cm}^2\mathrm{g}^1,$$
(19)
in $`T_{\mathrm{gr}}50\mathrm{K}`$, where $`T_{\mathrm{gr}}`$ is the effective grain temperature (see below). The mass fraction of grain is $`0.934\times 10^2`$ below the water vaporization temperature (100-200 K, depending on the density), and that of each component; water ice, volatile organics, refractory organics, troilite, orthopyroxene, olivine and metallic iron are $`1.19\times 10^3,6.02\times 10^4,3.53\times 10^3,5.69\times 10^4,7.33\times 10^4,2.51\times 10^3`$, and $`2.53\times 10^4`$, respectively (from Table 2 of Pollack et al. 1994). For the grain size distribution

$$n(a)\{\begin{array}{cc}a^{3.5}\hfill & \text{(}0.005\mu \mathrm{m}<a<1\mu \mathrm{m}\text{)}\hfill \\ a^{5.5}\hfill & \text{(}1\mu \mathrm{m}<a<5\mu \mathrm{m}\text{)}\hfill \end{array}$$
(20)
(Pollack et al. 1994) that is based on that derived by Mathis, Rumpl, & Nordsieck (1977), the energy transfer rate from the gas to the dust grains per unit mass is given by

$$_{\mathrm{gr}}=1.1\times 10^5n(\frac{f_{\mathrm{gr}}}{\rho _{\mathrm{gr}}})(\frac{T}{1000\mathrm{K}})^{1/2}[10.8\mathrm{exp}(75\mathrm{K}/T)](TT_{\mathrm{gr}})$$
(21)
(Hollenbach & McKee 1979), where $`T_{\mathrm{gr}}`$ is an effective grain temperature, which is determined by the energy balance for dust grains;

$$_{\mathrm{gr}}=4\sigma T_{\mathrm{gr}}^4\kappa _{\mathrm{gr}}\beta _{\mathrm{cont}},$$
(22)
where $`\kappa _{\mathrm{gr}}`$ is the Planck mean opacity owing to the dust grains and $`f_{\mathrm{gr}}/\rho _{\mathrm{gr}}`$ is the total volume of dust grains per unit mass of gas and $`f_{\mathrm{gr}}/\rho _{\mathrm{gr}}=5.3\times 10^3`$ at the lowest temperatures (i.e., below the water ice vaporization temperature). These quantities are taken from tables in Pollack et al.(1994).

The continuum energy transport rate $`\beta _{\mathrm{cont}}`$ decreases as $`\tau _{\mathrm{cont}}^2`$ owing to radiative diffusion in the optically thick case (e.g., Masunaga et al. 1998). Then,

$$\beta _{\mathrm{cont}}=\mathrm{min}(1,\tau _{\mathrm{cont}}^2).$$
(23)
#### 2.2.3 Continuum of Gas

Using the Planck mean opacity of gas $`\kappa _{\mathrm{gas}}`$, the cooling rate owing to the gas continuum is given by

$$_{\mathrm{cont}}=4\sigma T^4\kappa _{\mathrm{gas}}\beta _{\mathrm{cont}}.$$
(24)
We take the continuum Planck mean opacity for metal-free gas from Lenzuni et al. (1991), which includes all the important continuum processes, specifically, bound-free absorption by $`\mathrm{H}^0`$ and $`\mathrm{H}^{}`$; free-free absorption by $`\mathrm{H}^0`$, $`\mathrm{H}^{}`$, $`\mathrm{H}_2`$, $`\mathrm{H}_2^{}`$, $`\mathrm{H}_2^+`$, $`\mathrm{H}_3^+`$, $`\mathrm{He}^0`$, $`\mathrm{He}^{}`$; photodissociation of $`\mathrm{H}_2`$, and $`\mathrm{H}_2^+`$ by thermal radiation; Rayleigh scattering by $`\mathrm{H}^0`$, $`\mathrm{H}_2`$, $`\mathrm{He}^0`$; Thomson scattering by $`e^{}`$; and collision-induced absorption by $`\mathrm{H}_2`$ due to collisions with $`\mathrm{H}_2`$, $`\mathrm{He}^0`$, and $`\mathrm{H}^0`$. Among the above processes, the most important cooling mechanism is H2 collision-induced continuum, which dominates the cooling at about $`10^{16}\mathrm{cm}^3`$ for low metallicity (i.e., $`<10^6Z_{\mathrm{}}`$) clouds.

#### 2.2.4 Chemical Cooling/Heating

Following Hollenbach & McKee (1979), we assume the heat deposited per a formed molecular hydrogen as $`0.2+4.2(1+n_{\mathrm{cr}}/n)^1\mathrm{eV}`$ for H2 formation on grain surfaces (reaction 23 in Appendix), $`3.53(1+n_{\mathrm{cr}}/n)^1\mathrm{eV}`$ for H2 formation by H- process (reaction 8), $`1.83(1+n_{\mathrm{cr}}/n)^1\mathrm{eV}`$ for H2 formation by H$`{}_{}{}^{+}{}_{2}{}^{}`$ process (reaction 10), and $`4.48(1+n_{\mathrm{cr}}/n)^1\mathrm{eV}`$ for H2 formation by the three-body reactions (reactions 19 and 20), where

$$n_{\mathrm{cr}}=\frac{10^6T^{1/2}}{1.6y(\mathrm{H})\mathrm{exp}[(400/T)^2]+1.4y(\mathrm{H}_2)\mathrm{exp}[12000/(T+1200)]}\mathrm{cm}^3.$$
(25)
Collisional dissociation and ionization absorb the same energy as the binding energy, that is 4.48 eV per H2 dissociation, 13.6 eV per H ionization, 24.6 eV per He ionization, and 54.4 eV per He+ ionization.

### 2.3 Chemical reactions

We solve nonequilibrium chemistry involving the four elements H, He, C and O, that contains the following 45 species: $`\mathrm{H}`$, $`\mathrm{H}_2`$, $`\mathrm{e}^{}`$, $`\mathrm{H}^+`$, $`\mathrm{H}_2^+`$, $`\mathrm{H}_3^+`$, $`\mathrm{H}^{}`$, $`\mathrm{He}`$, $`\mathrm{He}^+`$, $`\mathrm{He}^{++}`$, $`\mathrm{HeH}^+`$, $`\mathrm{C}`$, $`\mathrm{C}_2`$, $`\mathrm{CH}`$, $`\mathrm{CH}_2`$, $`\mathrm{CH}_3`$, $`\mathrm{CH}_4`$, $`\mathrm{C}^+`$, $`\mathrm{C}_2^+`$, $`\mathrm{CH}^+`$, $`\mathrm{CH}_2^+`$, $`\mathrm{CH}_3^+`$, $`\mathrm{CH}_4^+`$, $`\mathrm{CH}_5^+`$, $`\mathrm{O}`$, $`\mathrm{O}_2`$, $`\mathrm{OH}`$, $`\mathrm{CO}`$, $`\mathrm{H}_2\mathrm{O}`$, $`\mathrm{HCO}`$, $`\mathrm{O}_2\mathrm{H}`$, $`\mathrm{CO}_2`$, $`\mathrm{H}_2\mathrm{CO}`$, $`\mathrm{H}_2\mathrm{O}_2`$, $`\mathrm{O}^+`$, $`\mathrm{O}_2^+`$, $`\mathrm{OH}^+`$, $`\mathrm{CO}^+`$, $`\mathrm{H}_2\mathrm{O}^+`$, $`\mathrm{HCO}^+`$, $`\mathrm{O}_2\mathrm{H}^+`$, $`\mathrm{H}_3\mathrm{O}^+`$, $`\mathrm{H}_2\mathrm{CO}^+`$, $`\mathrm{HCO}_2^+`$ and $`\mathrm{H}_3\mathrm{CO}^+`$. Chemical reactions included are listed in Appendix. The reactions of H and He chemistry are mainly taken from Abel et al. (1997) and Galli & Palla (1998). Other important H2 forming processes included are the three-body reactions (Palla et al. 1983) and the reaction on the surfaces of the dust grains (Tielens & Hollenbach 1985). In addition, we supplement H, He, C and O chemical reactions involving the above species from Millar, Farquhar, & Willacy (1997).

## 3 Results

In this section, we present the results obtained by the method described in §2 and discuss our analysis.

Figure 1 displays the evolutionary trajectories of protostellar clouds whose metallicities are $`z=0,10^6,10^4,10^2,1`$. Here we set the initial condition to be $`T=100\mathrm{K},n=1\mathrm{c}\mathrm{m}^3,y(e)=1\times 10^4,y(\mathrm{H}_2)=1\times 10^6`$. All the carbon is assumed to be in the form of CII, while oxygen is OI at the beginning.

The evolutionary trajectories for two other initial conditions are shown in Figure 2, as well as the same curves in Figure 1. We can see from Figure 2 that trajectories of clouds with a fixed composition converge rapidly toward the dense region for any initial conditions (e.g., Hayashi & Nakano 1965; Low & Lynden-Bell 1976). We have also tested the sensitivity to initial chemical compostions for two cases – (i) where the ionization degree $`y(e)=1\times 10^3`$ and (ii) all the carbon is assumed to be CI,– and have found the results similar to those shown in Figure 2. Hereafter, we discuss only on the case presented in Figure 1.

In general, as seen in Figure 1, the temperatures of lower metallicity clouds are higher because of their lower radiative cooling ability (i.e., lower radiative cooling rate for the same temperature and density) as long as the clouds are transparent to continuum and continue to collapse dynamically, owing to the efficient radiative cooling, in other words, during what is known as the “first collapse” stage.

On account of the total lack of metals and grains, the only cooling agent in the temperature range below $`10^4`$K for primordial ($`z=0`$) clouds is rovibrational transitions of molecular hydrogen (e.g., Matsuda et al. 1969). In the primordial gas, H2 is formed mainly by the H- process,

$`\mathrm{H}+e^{}`$ $``$ $`\mathrm{H}^{}+\gamma ;`$ (26)
$`\mathrm{H}+\mathrm{H}^{}`$ $``$ $`\mathrm{H}_2+e^{},`$ (27)
(Saslaw & Zipoy 1967) until the density reaches $`10^8\mathrm{cm}^3`$, where the three-body reactions

$$3\mathrm{H}\mathrm{H}_2+\mathrm{H}$$
(28)
$$2\mathrm{H}+\mathrm{H}_22\mathrm{H}_2$$
(29)
become efficient (Palla et al. 1983). H2 line emission via electric quadrupole transitions is the only efficient cooling agent until the number density exceeds $`10^{14}\mathrm{cm}^3`$, where H2 continuum emission via collision-induced dipole transitions begins to overwhelm. The drop in temperature at $`n10^{16}\mathrm{cm}^3`$ is due to the quadrature dependence of the collision-induced emission coefficient on the density. The cloud becomes optically thick to the continuum at $`10^{16}\mathrm{cm}^3`$. (See Fig. 3 a)

The gravitational contraction of metal-free protostellar clouds has been investigated with Omukai & Nishi (1998) by detailed hydrodynamical calculations assuming spherical symmetry. Our evolutionary trajectory for the primordial cloud agrees well with those of Omukai & Nishi (1998) at the number density $`n10^{10}\mathrm{cm}^3`$. At lower densities, the temperature of ours is higher owing to the influence of Omukai & Nishi’s (1998) initial condition. They started the calculation from a cloud in hydrostatic equilibrium at $`n10^6\mathrm{cm}^3`$. As long as the density is not so much higher than the initial value, the collapse is slower than the free-fall rate, which we have assumed in this paper. As a result, the temperature is lower because of the lower compressional heating.

The evolution of the cloud with $`z=10^6`$ is essentially the same as that of the primordial cloud, although there is minor deviation caused by extra cooling owing to water molecules and dust grains. (See Fig. 3 b) The dust thermal emission dominates the cooling just before the H2 continuum becomes effective. However, it is only temporary and the grains rapidly sublimate at $`n10^{14}\mathrm{cm}^3`$. The most refractory grain composition that contributes to the Planck mean opacity is iron, which sublimate at about 1200K at this density (Table 3 of Pollack et al. 1994). The slight amount of dust grains in the cloud with $`z=10^6`$ does not contribute to H2 formation significantly, as seen in Figure 4.

For the cloud with $`z=10^4`$, H2, which is formed mainly on the grain surfaces instead of the H- process, is the major coolant at low densities ($`<10^4\mathrm{cm}^3`$). Water molecules, which are mainly formed at $`10^6\mathrm{cm}^3`$ by the reaction

$$\mathrm{H}_2+\mathrm{OH}\mathrm{H}_2\mathrm{O}+\mathrm{H},$$
(30)
dominate the cooling in the range $`10^5\mathrm{cm}^3<n<10^{10}\mathrm{cm}^3`$. In present-day molecular cloud cores, usually neither H2O nor OH plays a key role in the thermal balance since the temperature there is lower than the excitation energy of these molecular species (Neufeld, Lepp, & Melnick 1995). However, in lower metal environments, the temperature is higher because of the lower radiative cooling. In that case, H2O can be a principal cooling agent. In the cloud with $`z=10^4`$, grain surface reactions exceed the H- process as a mode of H2 formation, although hydrogen does not become fully molecular until the three-body reaction begins to work. At $`10^{10}\mathrm{cm}^3`$, the grain thermal emission dominates the cooling and causes the temperature to drop to about 100K. The cloud becomes optically thick to the grain emission at about $`10^{13}\mathrm{cm}^3`$. (See Fig. 3 c)

For the cloud with $`z>10^2`$, cooling is dominated by atomic lines and CO lines in lower densities (that is, $`n<10^6\mathrm{cm}^3`$ for $`z=10^2`$, and $`n<10^3\mathrm{cm}^3`$ for $`z=1`$, respectively), and by the dust thermal emission in higher densities. (See Fig. 3 d, e) The temperature of our $`z=1`$ cloud drops far less than 10K, which is appropriate in present-day protostellar clouds. This is merely because we have ignored external radiation, and the results should not be taken seriously as a model of present-day molecular cloud cores.

Clouds with metallicity $`z>10^4`$ become opaque to the thermal radiation of dust grains before the grain vaporization occurs. Once the cloud becomes opaque, it begins to contract adiabatically since the radiative cooling rate drops rapidly. According to some hydrodynamical simulations (e.g, Larson 1969; Masunaga et al. 1998), a transient core in hydrostatic equilibrium (called a “first core”) forms after the cloud has more contractions. The hydrostatic core continues the adiabatic contraction by accreting the envelope matter until the temperature reaches $`2000`$ K, where molecular hydrogen begins to dissociate.

Clouds with lower metallicity (i.e., $`z<10^6`$) become opaque to the H2 continuum instead of the dust thermal emission, as described above. Since the temperature at that time is already near the dissociation value, molecular hydrogen begins to dissociate just after the cloud becomes opaque. As a result, no transient hydrostatic core forms in these cases as has been demonstrated by the hydrodynamical calculations of Omukai & Nishi (1998).

When the ratio of specific heats $`\mathrm{\Gamma }=\mathrm{dlog}p/\mathrm{dlog}\rho `$ exceeds the critical value 4/3, corresponding to a gradient of 1/3 in the density-temperature ($`nT`$) plane (i.e., Figure 1), the bulk motion of the cloud is decelerated and the free-fall assumption is no longer valid. However, the trajectories in the $`nT`$ plane are not altered since the clouds contract almost adiabatically anyway in the opaque stage because of the rapidly declining radiative cooling rate (e.g., Larson 1969; Narita, Nakano, & Hayashi 1970).

In the course of the adiabatic contraction of the transient cores, all the evolutionary trajectories converge to a certain line (Fig. 1, dashed line) in spite of their different composition and histories. We discuss the reason for the convergence in the following. The trajectories in the $`nT`$ plane just prior to when the clouds become opaque are determined by thermal balance between the compressional heating and the radiative cooling. The compressional heating rate (per unit mass) is given by

$$𝒢_{\mathrm{comp}}=p\frac{d}{dt}(\frac{1}{\rho })=\frac{c_\mathrm{s}^2}{t_{\mathrm{ff}}},$$
(31)
where $`p`$ is the pressure and $`c_\mathrm{s}`$ is the sound speed in the central region of the cloud. The radiative cooling at that time is dominated by continuum for all clouds, namely, the dust thermal emission for clouds with relative metallicity $`z>10^4`$, or the H2 continuum emission for those with lower metallicity. The radiative cooling rate is given by

$$_{\mathrm{rad}}=4\kappa _\mathrm{B}\sigma T^4,$$
(32)
where $`\kappa _\mathrm{B}`$ is the Planck mean continuum opacity. In the equation above, we have used the relation that the grain temperature is the same as the gas temperature, $`T_{\mathrm{gr}}=T`$, which is valid in a cloud so dense that it becomes opaque to the continuum. Note that expression (32) holds for any clouds that cool mainly by continuum by choosing an appropriate value of opacity, where the opacity is an increasing function of relative metallicity $`z`$. Equating (31) and (32), for each value of $`\kappa _\mathrm{B}(z)`$, we obtain a trajectory just before the cloud becomes opaque (we call it Tr($`\kappa _\mathrm{B}`$));

$$\mathrm{Tr}(\kappa _\mathrm{B}):\frac{c_\mathrm{s}^2}{t_{\mathrm{ff}}}=4\kappa _\mathrm{B}\sigma T^4.$$
(33)
As the cloud of opacity $`\kappa _\mathrm{B}`$ collapses along the trajectory represented by (33), the central part of the cloud becomes opaque to the continuum at the point

$$\mathrm{P}(\kappa _\mathrm{B}):\tau _\mathrm{J}(\kappa _\mathrm{B},n,T)=\kappa _\mathrm{B}\rho \pi \frac{c_\mathrm{s}}{\sqrt{G\rho }}=1,$$
(34)
where we have used the assumption that the size of the central region is approximately $`\lambda _\mathrm{J}`$. Eliminating $`\kappa _\mathrm{B}`$ in equations (33) and (34), we obtain the locus L of points P for variable $`\kappa _\mathrm{B}`$;

$$\mathrm{L}:T=(\frac{k^3}{12\sigma ^2m_\mathrm{p}})^{1/5}n^{2/5}.$$
(35)
On the other hand, the trajectory of the cloud of opacity $`\kappa _\mathrm{B}`$ in the opaque stage, O($`\kappa _\mathrm{B}`$), is the line that starts from the point P($`\kappa _\mathrm{B}`$) and whose gradient is approximately equal to the value for adiabatic contraction $`\gamma _{\mathrm{ad}}1`$. The adiabatic exponent $`\gamma _{\mathrm{ad}}`$ for molecular hydrogen is 7/5 in the temperature range $`T200`$K, while it is 2/3 in the lowest temperature. The gradient of the trajectory O($`\kappa _\mathrm{B}`$) in the range $`T200`$K is 2/5. This value happens to be the same as the gradient of the line L. Thus, the trajectories in the opaque stage O all coincide with the line L, which can be written only by physical constants. When the clouds become opaque (i.e., at the point P), their specific entropy has the same value for all clouds regardless of their metallicity.

As a cloud climbs up along line L, the temperature reaches the dissociation value ($`2000`$K) at $`n3\times 10^{16}\mathrm{cm}^3`$. The ratio of specific heats are reduced below the critical value 4/3. Then the cloud begins to collapse dynamically again, or in other words, it begins the second collapse (e.g., Larson 1969). The trajectories in the second collapse stage are also unique because the specific entropy has the same value for all clouds regardless of their metallicity. When most hydrogen molecules have been dissociated, the ratio of specific heats rises above the critical value 4/3 again. After the central part of the cloud contracts almost adiabatically to some extent, a hydrostatic core forms. This core is referred to as the stellar or second core in literature. We cannot accurately estimate the size of the stellar core from our simple one-zone treatment. However, hydrodynamical calculations tell us that the number density and mass of the stellar core at its formation time are on the order of $`10^{22}\mathrm{cm}^3`$ and $`10^3M_{\mathrm{}}`$ respectively, both for present-day (Larson 1969) and primordial protostars (Omukai & Nishi 1998). This value is expected to be universal, namely, independent of metallicity since the equations of state, or trajectories in the $`nT`$ plane, are the same in the second collapse stage.

The stellar core, or protostar, grows in mass by accretion of the envelope material. The temperature rises accordingly and eventually the core will become an ordinary star. Consequently, the final mass of stars is determined by the subsequent accretion onto the core, although the size of the protostar at the time of formation is the same. Stahler, Shu, & Taam (1980) argued that a rough estimate of the protostellar mass accretion rate $`\dot{M}`$ can be obtained from the relation

$$\dot{M}c_\mathrm{s}^3/G,$$
(36)
where $`c_\mathrm{s}`$ is the isothermal sound speed in the initial protostellar cloud. Therefore, the mass accretion rate is higher for a lower metallicity cloud because of the higher temperature of the protostellar cloud (Stahler, Palla, & Salpeter 1986). The higher accretion rate and lower opacity (then smaller radiation force) may result in the higher mass of formed stars for a lower metallicity cloud (e.g., Wolfire & Cassinelli 1987). However, to address these issues fully, further investigations are needed.

## 4 Summary

We have investigated the thermal and chemical evolution of collapsing protostellar clouds with different metallicities. The varied range of metallicity spans the local interstellar value ($`1Z_{\mathrm{}}`$) to zero.

The evolution of the clouds is summarized as follows. While the clouds are transparent to continuous radiation, the temperature of clouds with lower metallicity is higher since their radiative cooling rates are lower for the same density and temperature. However, after the clouds become opaque and begin adiabatic contraction, their evolutionary trajectories converge to a line that is determined only by physical constants. Thereafter, the trajectories coincide with each other regardless of their metallicity. Consequently, the physical dimension of the stellar core at the time of formation is the same for clouds with any composition.

We have also discussed analytically the reasons for the convergence.

The author wishes to express his cordial thanks to Humitaka Sato and Naoshi Sugiyama for their continual interest and advice, to Hiroshi Koyama, Ryoichi Nishi, and Toru Tsuribe for useful discussions, and to the referee for improving this manuscript.


text:
# References

A Novel Approach to Spherical Harmonics Expansion for Electron Transport in Semiconductors

S. F. Liotta and A. Majorana

Dipartimento di Matematica - University of Catania -

Viale A. Doria 6, 95125 Italy

Abstract. A set of equations is derived from the Boltzmann kinetic equation describing charge transport in semiconductors. The unknowns of these equations depend on the space-time coordinates and the electron energy. The non-parabolic and parabolic band approximation are treated in detail. In these cases, the set of equations is equivalent to those obtained in the spherical-harmonics expansion. Stationary and homogeneous solutions are explicitly treated. In order to solve numerically the equations, a cut in the energy range is introduced. The modified model maintains the physical characteristics of the original equations. The solution of an asymptotic equation is found and compared to the numerical solutions.
Keywords: Boltzmann equation; Semiconductors; Spherical harmonics expansion
1. Introduction
The motion of electrons in a crystal is governed by complex physical laws so that accurate models are required to achieve correct results. The drift-diffusion model has been widely used in the past, and it has provided a good description of relevant physical mechanisms. In modern devices, whose sizes are in the submicron range, non-equilibrium effects play an important role and are not adequately modeled by the drift-diffusion approach. Since in many processes a more accurate description than the hydrodynamical setting is required, Boltzmann equation or Monte Carlo simulations are employed. A fully kinetic treatment of carrier dynamics guarantees accurate results but requires very expensive numerical procedures in order to solve realistic problems.

To reduce the complexity of the use of the full Boltzmann equation, many authors (see ref.  and references therein) have introduced simpler models, assuming particular forms for the probability density function. The aim is to yield results, which are more accurate than those obtained by hydrodynamical models but less expensive than direct particle simulations or numerical treatment of the Boltzmann equation. In this framework, a well-known model is derived using a spherical-harmonics expansion , , ,  from the Boltzmann equation.

The paper is organized as follows. In Sec. 2, we introduce the Boltzmann Equation. The collisional operator takes into account the interactions between electrons and phonons, i.e. vibrations of the lattice. This is assumed to be in thermal equilibrium. A set of kinetic equations is derived in Sec. 3. These are obtained for a general expression of the microscopic kinetic energy, so that, for example, they can be used also for non-parabolic bands. The equations are of partial differential and difference type. Equations are specialized in Sec. 4 in the parabolic band approximation, and is Sec. 5, where the Kane band model is assumed. In both cases we show that the equations are equivalent to those obtained in the spherical harmonics framework. In Sec. 6 we write the equations in the stationary and homogeneous case. The rest of the paper is devoted to the study of these equations. In Sec. 7 an asymptotic equation, valid for large value of the electric field, is found and analyzed; the explicit solution is found. In Sec. 8 we discuss some difficulties related to the infinite energy range and we propose a suitable solution by introducing a cut in the scattering rate of the collisional operator. The numerical scheme to solve the equations is presented in Sec. 9, and the results are shown in the last section.

Throughout the paper, boldface and lightface symbols denote vectors and scalar quantities respectively.
2. Basic equations
We consider an electron gas, which moves in a lattice, subjects to an external electric field $`𝐄`$. This can be applied or related, by Poisson equation, to the density of the gas. For the electron gas in a semiconductor device the Boltzmann equation ,  writes

(1)
$$\frac{f}{t}+\frac{1}{\mathrm{}}_𝐤\epsilon _𝐱f\frac{\text{e}}{\mathrm{}}𝐄_𝐤f=Q(f),$$

where the unknown function $`f(t,𝐱,𝐤)`$ represents the probability of finding an electron at the position $`𝐱\text{}^3`$, with the wave-vector $`𝐤\text{}^3`$, at time $`t`$. The parameters $`\mathrm{}`$ and e are the Planck constant divided by $`2\pi `$ and the positive electric charge, respectively. The symbol $`_𝐤`$ stands for the gradient with respect to the variables $`𝐤`$ and $`_𝐱`$ that with respect to the space coordinates $`𝐱`$. The particle energy $`\epsilon `$ is an assigned nonnegative continuous function, defined on $`\text{}^3`$, of the wave-vector $`𝐤`$. Here we assume that the low-density approximation holds, so that $`Q`$ is linear in $`f`$. If $`K(𝐤,\stackrel{~}{𝐤})`$ is the sum of the scattering kernels, which describe the nature of the inelastic collisions (for example, electron–polar or non-polar optical phonon scattering), and $`K_0(𝐤,\stackrel{~}{𝐤})`$ is the kernel for the elastic collisions (for example, electron–impurity scattering), then the collisional operator writes , ,

(2) $`Q(f)=(𝗇_q+1){\displaystyle _\text{}^3}K(𝐤,\stackrel{~}{𝐤})\delta (\epsilon (\stackrel{~}{𝐤})\epsilon (𝐤)\mathrm{}\omega )f(t,𝐱,\stackrel{~}{𝐤})𝑑\stackrel{~}{𝐤}`$
$`+𝗇_q{\displaystyle _\text{}^3}K(\stackrel{~}{𝐤},𝐤)\delta (\epsilon (\stackrel{~}{𝐤})\epsilon (𝐤)+\mathrm{}\omega )f(t,𝐱,\stackrel{~}{𝐤})𝑑\stackrel{~}{𝐤}`$
$`+{\displaystyle _\text{}^3}K_0(𝐤,\stackrel{~}{𝐤})\delta (\epsilon (\stackrel{~}{𝐤})\epsilon (𝐤))f(t,𝐱,\stackrel{~}{𝐤})𝑑\stackrel{~}{𝐤}\overline{\nu }(𝐤)f(t,𝐱,𝐤),`$
where

(3) $`\overline{\nu }(𝐤)=𝗇_q{\displaystyle _\text{}^3}K(𝐤,\stackrel{~}{𝐤})\delta (\epsilon (\stackrel{~}{𝐤})\epsilon (𝐤)\mathrm{}\omega )𝑑\stackrel{~}{𝐤}`$
$`+(𝗇_q+1){\displaystyle _\text{}^3}K(\stackrel{~}{𝐤},𝐤)\delta (\epsilon (\stackrel{~}{𝐤})\epsilon (𝐤)+\mathrm{}\omega )𝑑\stackrel{~}{𝐤}+{\displaystyle _\text{}^3}K_0(\stackrel{~}{𝐤},𝐤)\delta (\epsilon (\stackrel{~}{𝐤})\epsilon (𝐤))𝑑\stackrel{~}{𝐤}.`$
Often, in the following we omit to write explicitly that $`f`$ depends also on $`t`$ and $`𝐱`$, whereas we do not operate with these variables.
The constant quantity $`𝗇_q`$ represents the thermal-equilibrium number of phonons and is given by

$$𝗇_q=\frac{1}{\mathrm{exp}\left(\frac{\mathrm{}\omega }{k_BT_L}\right)1}$$

where $`k_B`$ is the Boltzmann constant and $`T_L`$ is the constant lattice temperature.

The symbol $`\delta (\epsilon (\stackrel{~}{𝐤})\epsilon (𝐤)\pm \mathrm{}\omega )`$ means the composition of the real-valued function $`\epsilon (\stackrel{~}{𝐤})\epsilon (𝐤)\pm \mathrm{}\omega `$ and the Dirac distribution $`\delta `$. Since the function $`\epsilon `$ may have many different expressions, it is not possible in general, by using a unique explicit technique, to transform each of the integral operators (S0.Ex2) in equivalent operators without Dirac distributions. Therefore, we shall maintain the original form of (S0.Ex2) in the following.

As usual, due to the isotropy of phase-space, the kernels $`K(𝐤,\stackrel{~}{𝐤})`$ and $`K_0(𝐤,\stackrel{~}{𝐤})`$ are symmetric with respect to the wave-vectors $`𝐤`$ and $`\stackrel{~}{𝐤}`$. Moreover, since the kernels are scalar, they depend on $`𝐤`$ and $`\stackrel{~}{𝐤}`$ only through the scalar quantities $`\epsilon (𝐤)`$, $`\epsilon (\stackrel{~}{𝐤})`$ and $`𝐤\stackrel{~}{𝐤}`$. This implies in particular that the collision frequency $`\nu `$ depends on $`𝐤`$ only through the variable $`\epsilon `$.
3. Energy-kinetic equations
In the framework of the Boltzmann equation for a perfect gas, a classical procedure, due to Grad , to reduce the dimension of the space of the independent coordinates, is the moment method. It consists in expanding the ratio between the probability density $`f`$ and a local maxwellian function. Inserting this expansion in the Boltzmann equation, an infinite sequence of differential equations are obtained. In general no finite set of equations is uncoupled from the rest of the system. To obtain a determined system up to certain order $`m`$ (usually 13) only $`m`$ equation are considered and the coefficients of the expansion of higher order are assumed negligible.

In order to obtain a new set of equations, where the unknown functions depend on the space-time coordinates and the microscopic kinetic energy, we follow the scheme analyzed in . This scheme partially recalls the moment method.

If $`\text{P}(𝐤)`$ is a polynomial and $`u`$ a real variable, we multiply the Boltzmann equation (1) by $`\text{P}(𝐤)\delta (\epsilon (𝐤)u)`$ and then formally integrate with respect to the variable $`𝐤`$ over the whole space $`\text{}^3`$. It is possible to show  that the following equation can be derived:

(4) $`{\displaystyle \frac{\text{ }}{t}}{\displaystyle _\text{}^3}f(t,𝐱,𝐤)\text{P}(𝐤)\delta (\epsilon (𝐤)u)𝑑𝐤`$
$`+_𝐱{\displaystyle _\text{}^3}𝐯(𝐤)f(t,𝐱,𝐤)\text{P}(𝐤)\delta (\epsilon (𝐤)u)𝑑𝐤`$
$`\text{e}𝐄{\displaystyle \frac{\text{ }}{u}}{\displaystyle _\text{}^3}f(t,𝐱,𝐤)\text{P}(𝐤)𝐯(𝐤)\delta (\epsilon (𝐤)u)𝑑𝐤`$
$`+{\displaystyle \frac{\text{e}}{\mathrm{}}}𝐄{\displaystyle _\text{}^3}f(t,𝐱,𝐤)\left[_𝐤\text{P}(𝐤)\right]\delta (\epsilon (𝐤)u)𝑑𝐤={\displaystyle _\text{}^3}Q(f)\text{P}(𝐤)\delta (\epsilon (𝐤)u)𝑑𝐤.`$
Here the molecular velocity is defined by means of the formula

(5)
$$𝐯(𝐤):=\frac{1}{\mathrm{}}_𝐤\epsilon (𝐤).$$

For any choice of the function $`\text{P}(𝐤)`$ a new equation is obtained. As it happens for the moment method, each equation contains, except for particular cases  at least two unknown functions. For example, if we choose $`\text{P}(𝐤)1`$, then the left hand side of eq. (S0.Ex7) contains in general the following unknowns

(6) $`N(t,𝐱,u)`$ $`:=`$ $`{\displaystyle _\text{}^3}f(t,𝐱,𝐤)\delta (\epsilon (𝐤)u)𝑑𝐤,`$
(7) $`𝐕(t,𝐱,u)`$ $`:=`$ $`{\displaystyle _\text{}^3}f(t,𝐱,𝐤)𝐯(𝐤)\delta (\epsilon (𝐤)u)𝑑𝐤.`$
The scalar quantity $`N(t,𝐱,u)`$ is the probability density function to find an electron with energy $`u`$ at time $`t`$ and position $`𝐱`$. The function $`N`$ is non-negative, because the probability density $`f`$ is always non-negative. We note that the presence of the delta-function implies that $`N`$ and $`𝐕`$ vanish for all $`(t,𝐱,u)`$ such that $`uu_0=\mathrm{min}\left\{\epsilon (𝐤):𝐤\text{}^3\right\}`$. If the function $`N`$ is known, then we can obtain all the hydrodynamical scalars. For example, to find the hydrodynamical density $`\rho (t,𝐱)`$, it is sufficient to integrate $`N`$ with respect to all values of energy $`u`$; i.e.

$$\rho (t,𝐱):=_{u_0}^+\mathrm{}N(t,𝐱,u)𝑑u.$$

The physical meaning of the other quantity $`𝐕(t,𝐱,u)`$ can be understood with a similar argument.

In the rest of this article we assume the function $`\epsilon (𝐤)`$ to be even (i.e. $`\epsilon (𝐤)=\epsilon (𝐤)`$) and the kernels $`K`$ and $`K_0`$ to depend only on $`\epsilon (𝐤)`$ and $`\epsilon (\stackrel{~}{𝐤})`$. The first hypothesis is verified for the most common models used in the numerical simulations. The second are introduced in order to simplify the treatment of the equations. Let

$$\text{K}(\epsilon (𝐤),\epsilon (\stackrel{~}{𝐤})):=K(𝐤,\stackrel{~}{𝐤})\text{ and }\text{K}_0(\epsilon (𝐤),\epsilon (\stackrel{~}{𝐤})):=K_0(𝐤,\stackrel{~}{𝐤}).$$

Then eqs. (S0.Ex2)-(S0.Ex3) become

$`Q(f)`$ $`=`$ $`(𝗇_q+1)\text{K}(\epsilon (𝐤),\epsilon (𝐤)+\mathrm{}\omega )N(t,𝐱,\epsilon (𝐤)+\mathrm{}\omega )`$
$`+𝗇_q\text{K}(\epsilon (𝐤)\mathrm{}\omega ,\epsilon (𝐤))N(t,𝐱,\epsilon (𝐤)\mathrm{}\omega )`$
$`+\text{K}_0(\epsilon (𝐤),\epsilon (𝐤))N(t,𝐱,\epsilon (𝐤))\nu (\epsilon (𝐤))f(t,𝐱,𝐤)`$
with

$`\nu (\epsilon (𝐤))`$ $`=`$ $`(𝗇_q+1)\text{K}(\epsilon (𝐤)\mathrm{}\omega ,\epsilon (𝐤))\sigma (\epsilon (𝐤)\mathrm{}\omega )`$
$`+𝗇_q\text{K}(\epsilon (𝐤),\epsilon (𝐤)+\mathrm{}\omega )\sigma (\epsilon (𝐤)+\mathrm{}\omega )+\text{K}_0(\epsilon (𝐤),\epsilon (𝐤))\sigma (\epsilon (𝐤)),`$
where

(10)
$$\sigma (u):=_\text{}^3\delta (\epsilon (𝐤)u)𝑑𝐤$$

is the density of states.
We are interested in deducing two equations in the only variables $`N`$ and $`𝐕`$. By choosing $`\text{P}(𝐤)1`$ and $`\text{P}(𝐤)=𝐯(𝐤)`$ eqs. (S0.Ex7), (S0.Ex11) and (S0.Ex12) give the following equations

(11) $`{\displaystyle \frac{N}{t}}+_𝐱𝐕\text{e}𝐄{\displaystyle \frac{𝐕}{u}}=G(N)`$
(12) $`{\displaystyle \frac{V_i}{t}}+{\displaystyle \underset{j=1}{\overset{3}{}}}\left[{\displaystyle \frac{\mathrm{\Pi }_{ij}}{x_j}}\text{e}E_j{\displaystyle \frac{\mathrm{\Pi }_{ij}}{u}}+{\displaystyle \frac{\text{e}}{\mathrm{}}}E_j\mathrm{\Delta }_{ij}\right]=\nu (u)V_i(i=1,2,3),`$
where

(13) $`G(N):=(𝗇_q+1)\text{K}(u,u+\mathrm{}\omega )\sigma (u)N(t,𝐱,u+\mathrm{}\omega )`$
$`+𝗇_q\text{K}(u,u\mathrm{}\omega )\sigma (u)N(t,𝐱,u\mathrm{}\omega )`$
$`\left[𝗇_q\text{K}(u,u+\mathrm{}\omega )\sigma (u+\mathrm{}\omega )+(𝗇_q+1)\text{K}(u,u\mathrm{}\omega )\sigma (u\mathrm{}\omega )\right]N(t,𝐱,u),`$
(14) $`\nu (u)=𝗇_q\text{K}(u,u+\mathrm{}\omega )\sigma (u+\mathrm{}\omega )`$
$`+(𝗇_q+1)\text{K}(u,u\mathrm{}\omega )\sigma (u\mathrm{}\omega )+\text{K}_0(u,u)\sigma (u),`$
(15) $`\mathrm{\Pi }_{ij}:={\displaystyle _\text{}^3}f(t,𝐱,𝐤)v_i(𝐤)v_j(𝐤)\delta (\epsilon (𝐤)u)𝑑𝐤,`$
(16) $`\mathrm{\Delta }_{ij}:={\displaystyle _\text{}^3}f(t,𝐱,𝐤){\displaystyle \frac{v_i(𝐤)}{k_j}}\delta (\epsilon (𝐤)u)𝑑𝐤.`$
Also the tensor $`\mathrm{\Delta }_{ij}`$ is symmetric, due to eq. (5). We note that the effects of the elastic collisions, through the kernel $`\text{K}_0`$, disappear in the operator $`G(N)`$.
The set of eqs. (11)-(12) contains $`N`$ and $`𝐕`$ but also two tensors $`\mathrm{\Delta }_{ij}`$ and $`\mathrm{\Pi }_{ij}`$. Therefore, except in the case of spatially homogeneous solutions with null electric field, we need to assume some relations (usually called constitutive equations in the thermodynamic framework), which link $`\mathrm{\Delta }_{ij}`$ and $`\mathrm{\Pi }_{ij}`$ with $`N`$ and $`𝐕`$. These relations should also depend on the form of the microscopic kinetic energy $`\epsilon `$.

We feel that at this stage it would be worthwhile to make a simple observation. We note that the linearity of the Boltzmann equation (1) is maintained in eqs. (11)-(12). To maintain this feature, we suggest to assume constitutive relations of the kind

(17) $`\mathrm{\Pi }_{ij}(t,𝐱,u)`$ $`=`$ $`p_{ij}(u)N(t,𝐱,u)+p_i(u)V_j(t,𝐱,u)+p_j(u)V_i(t,𝐱,u)`$
(18) $`\mathrm{\Delta }_{ij}(t,𝐱,u)`$ $`=`$ $`d_{ij}(u)N(t,𝐱,u)+d_i(u)V_j(t,𝐱,u)+d_j(u)V_i(t,𝐱,u),`$
where the symmetric tensors $`p_{ij}`$ and $`d_{ij}`$ and the vectors $`p_i`$ and $`d_i`$ must be determined. In the next sections we propose a simple choice in the cases of the parabolic and non-parabolic band approximation.
4. The parabolic case
The simplest and widely used expression for the microscopic kinetic energy  is

(19)
$$\epsilon (𝐤)=\frac{\mathrm{}^2}{2m^{}}k^2,$$

where $`m^{}`$ is the value of the effective electron mass in the parabolic mass approximation. In this case the tensor $`\mathrm{\Delta }_{ij}`$ reduces to

(20)
$$\mathrm{\Delta }_{ij}=\frac{\mathrm{}}{m^{}}N(t,𝐱,u)\delta _{ij}(i,j=1,2,3),$$

where $`\delta _{ij}`$ is the Kroneker symbol; so that no further assumption on eq. (18) is required. The tensor $`\mathrm{\Pi }_{ij}`$ is determined by assuming that (17) holds with $`p_i(u)0`$ and moreover that the constitutive relation is exact if $`f`$ depends on $`𝐤`$ only through the variable $`\epsilon `$. For $`f(t,𝐱,𝐤)=\stackrel{~}{f}(t,𝐱,\epsilon (𝐤))`$, we have

$`N(t,𝐱,u)`$ $`=`$ $`\stackrel{~}{f}(t,𝐱,u)\sigma (u)`$
$`\mathrm{\Pi }_{ij}(t,𝐱,u)`$ $`=`$ $`\stackrel{~}{f}(t,𝐱,u){\displaystyle _\text{}^3}v_i(𝐤)v_j(𝐤)\delta (\epsilon (𝐤)u)d𝐤(i,j=1,2,3),`$
so that

$$\sigma (u)p_{ij}(u)=_\text{}^3v_i(𝐤)v_j(𝐤)\delta (\epsilon (𝐤)u)𝑑𝐤.$$

Now, it is easy to verify that

$$\sigma (u)=4\sqrt{2}\pi \left(\frac{\sqrt{m^{}}}{\mathrm{}}\right)^3\theta (u)\sqrt{u},$$

where $`\theta `$ is the Heaviside step function. Then, for every $`u0`$,

$$_\text{}^3v_i(𝐤)v_j(𝐤)\delta (\epsilon (𝐤)u)d𝐤=\frac{2}{3}\frac{u}{m^{}}\delta _{ij}\sigma (u)(i,j=1,2,3).$$

Hence, eq. (12) becomes

(21)
$$\frac{𝐕}{t}+\frac{1}{3m^{}}\left[2u_𝐱N2\text{e}𝐄\frac{\text{ }}{u}(uN)+3\text{e}𝐄N\right]=\nu (u)𝐕.$$

In the following the functions $`N`$ and $`𝐕`$ will be the unknowns, which must satisfy eqs. (11) and (21). It is evident that the knowledge of these quantities cannot provide the probability density $`f`$ uniquely. Also, it is very hard to prove that, if $`N`$ and $`𝐕`$, at the initial time, derive from a probability density $`f`$, there exists a non-negative probability density, which gives $`N`$ and $`𝐕`$ for all times.
Actually the situation is similar to that of the moments in the hydrodynamical equations derived by the Grad method. Here usually, not even for the initial data, the existence of a probability density, which gives all the moments occurring in the equations, is considered.

Therefore, we limit to assuming only the simple conditions that $`N`$ is nonnegative and both the variables $`N`$ and $`𝐕`$ vanish as $`u0`$.

We note that these unknowns are simply related to the first two terms of the spherical harmonic expansion of the distribution function

(22)
$$f(t,𝐱,𝐤)=f_0(t,𝐱,\epsilon (𝐤))+𝐤𝐟_1(t,𝐱,\epsilon (𝐤))+\mathrm{}.$$

It is a simple matter to show that

$`f_0(t,𝐱,\epsilon (𝐤))`$ $`=`$ $`{\displaystyle \frac{N(t,𝐱,\epsilon (𝐤))}{\sigma (\epsilon (𝐤))}},`$
$`𝐟_1(t,𝐱,\epsilon (𝐤))`$ $`=`$ $`{\displaystyle \frac{3\mathrm{}}{2\epsilon (𝐤)\sigma (\epsilon (𝐤))}}𝐕(t,𝐱,\epsilon (𝐤)).`$
The corresponding equations for $`f_0`$ and $`𝐟_1`$ (see ref. ) coincide with our eqs. (11) and (21). This is a consequence of the choice on the approximation of the tensor $`\mathrm{\Pi }_{ij}`$. It is evident that different choices give different sets of equations. In this paper we do not investigate this possibility.

Similar equations were proposed by Hänsch in ref.  in the framework of quantum kinetic transport. Also in that case a problem of equating the number of equations and unknowns arose. It was solved by simple physical arguments, but the equations are different from eqs. (11)-(21).
5. The non-parabolic band approximation
For large values of electron energy eq. (19) is not adequate. The following suitable expression was given by Kane :

(23)
$$\epsilon (1+\alpha _k\epsilon )=\frac{\mathrm{}^2}{2m^{}}k^2,$$

where $`\alpha _k`$ is a constant. Now, it is simple to verify that

$$\sigma (u)=4\sqrt{2}\pi \left(\frac{\sqrt{m^{}}}{\mathrm{}}\right)^3\theta (u)\sqrt{u(1+\alpha _ku)}(1+2\alpha _ku).$$

So eq. (11) can be written explicitly. The second equation requires only to determine $`\mathrm{\Pi }_{ij}`$ and $`\mathrm{\Delta }_{ij}`$ according to eqs. (17)-(18). Due to the spherical symmetry of the band we assume $`p_j(u)=d_j(u)=0`$. Analogously to the parabolic case, we require that eqs. (17)-(18) hold if $`f`$ depends on $`𝐤`$ only through $`\epsilon `$. In the same manner we can see that that

$$p_{ij}(u)=g(u)\delta _{ij},d_{ij}(u)=h(u)\delta _{ij}$$

where

$$g(u)=\frac{2}{3}\frac{u(1+\alpha _ku)}{m^{}(1+2\alpha _ku)^2}$$

and

$$h(u)=\frac{1}{m^{}(1+2\alpha _ku)}\frac{4\alpha _k}{3m^{}}\frac{u(1+\alpha _ku)}{(1+2\alpha _ku)^3}.$$

Therefore eq. (12) becomes

(24)
$$\frac{𝐕}{t}+g(u)_𝐱N\text{e}𝐄\frac{\text{ }}{u}\left[g(u)N\right]+\text{e}𝐄h(u)N=\nu (u)𝐕.$$

Also in this case the equations for $`N`$ and $`V`$ are equivalent to those obtained by means of the spherical harmonics expansion . For $`\alpha _k=0`$ we recover eq. (21).
6. Homogeneous solutions with constant electric field
We consider equations (11) and (24), and we look for a solution depending only on $`u`$. Since in this case the electric field must be a constant vector, we choose the reference frame so that only the first component of $`𝐄`$ is different from zero. As a consequence the only significant component of $`𝐕`$ is the first. We limit to assume constant kernels K and $`\text{K}_0`$. This fact holds when a homogeneous, intrinsic silicon at room temperature is considered.

It is useful to introduce dimensionless variables. Let be

$$t_{}:=\left[4\sqrt{2}\pi \left(\frac{\sqrt{m^{}}}{\mathrm{}}\right)^3\sqrt{k_BT_L}\text{K}𝗇_q\right]^1,\mathrm{}_{}:=\sqrt{\frac{k_BT_L}{m^{}}}t_{},u_{}:=k_BT_L,$$

$$w:=\frac{u}{u_{}},n(w):=u_{}\mathrm{}_{}^3N(u),v(w):=u_{}\mathrm{}_{}^2t_{}V_1(u),$$

$$\alpha :=\frac{\mathrm{}\omega }{k_BT_L},a:=\frac{𝗇_q+1}{𝗇_q}=e^\alpha ,\beta =\frac{\text{K}_0}{𝗇_q\text{K}},\stackrel{~}{\alpha }:=\alpha _ku_{},\zeta :=\text{e}E\frac{\mathrm{}_{}}{u_{}}.$$

It is a simple matter that the following equations can be obtained

(25) $`\zeta {\displaystyle \frac{dv}{dw}}=s(w)\left[an(w+\alpha )+n(w\alpha )\right]\left[s(w+\alpha )+as\left(w\alpha \right)\right]n(w)`$
(26) $`\zeta \left[{\displaystyle \frac{d\text{ }}{dw}}(r(w)n)+q(w)n\right]=\left[s(w+\alpha )+as\left(w\alpha \right)+\beta s(w)\right]v(w),`$
where

(27) $`s(w)`$ $`=`$ $`\theta (w)\sqrt{w(1+\stackrel{~}{\alpha }w)}(1+2\stackrel{~}{\alpha }w)`$
(28) $`r(w)`$ $`=`$ $`{\displaystyle \frac{2}{3}}{\displaystyle \frac{w(1+\stackrel{~}{\alpha }w)}{(1+2\stackrel{~}{\alpha }w)^2}}`$
(29) $`q(w)`$ $`=`$ $`{\displaystyle \frac{1}{1+2\stackrel{~}{\alpha }w}}{\displaystyle \frac{4\stackrel{~}{\alpha }}{3}}{\displaystyle \frac{w(1+\stackrel{~}{\alpha }w)}{(1+2\stackrel{~}{\alpha }w)^3}}.`$
As previously mentioned, we recall that we assume the conditions

$$n(0)=v(0)=0.$$

We note that in this case eq. (26) gives $`v(0)=0`$ as a consequence of the boundary condition $`n(0)=0`$.

We look for solutions of eqs. (25)-(26) satisfying the following conditions

(30) $`n(0)=0,\underset{w+\mathrm{}}{lim}v(w)=0`$
(31) $`n(w)0\text{for every }w0\text{ and }{\displaystyle _{\text{}}}n(w)𝑑w>0.`$
Since eqs. (25)-(26) are linear and homogeneous, if a solution $`(n(w),v(w))`$, satisfying the above conditions exists, then, for every $`c>0`$, also $`(cn(w),cv(w))`$ is solution.

It is important to note that the conservation of particle number holds for the system (25)-(26). In fact, from eq. (25) it is easy to verify that

(32)
$$_0^+\mathrm{}s(w)\left[an(w+\alpha )+n(w\alpha )\right]\left(s(w+\alpha )+as(w\alpha )\right)n(w)dw=0,$$

for any nonnegative function $`n`$ such that $`n(w)s(w)`$ is integrable and $`n(0)=0`$. This property forces to choose the boundary on $`v`$ in (30), being $`v(0)=0`$.
7. Asymptotic equations
A simple approximate solution of the previous problem can be obtained in the parabolic band case. We introduce a new variable $`\psi `$ defined by

$$n(w)=\sqrt{w}\psi (w).$$

Then, taking into account that now $`\stackrel{~}{\alpha }=0`$, eqs. (25)-(26) become

(33) $`\zeta {\displaystyle \frac{dv}{dw}}=\sqrt{w}\left[a\sqrt{w+\alpha }\psi (w+\alpha )+\sqrt{(w\alpha )_+}\psi ((w\alpha )_+)\right]`$
$`\text{ }\sqrt{w}\left(\sqrt{w+\alpha }+a\sqrt{(w\alpha )_+}\right)\psi (w)`$
(34) $`{\displaystyle \frac{2}{3}}\zeta \sqrt{w^3}{\displaystyle \frac{d\psi }{dw}}=\left(\sqrt{w+\alpha }+a\sqrt{(w\alpha )_+}+\beta \sqrt{w}\right)v(w),`$
where $`(z)_+=\mathrm{max}\{z,0\}`$.
We look for asymptotic equations of system (S0.Ex32)-(34) for large values of the energy $`w`$. To this scope we expand the coefficients of the equations up to the zeroth order (for example $`\sqrt{w+\alpha }\sqrt{w}`$). We obtain a new set of equations

(35) $`\zeta {\displaystyle \frac{dv_A}{dw}}`$ $`=`$ $`w\left\{\left[a\psi _A(w+\alpha )+\psi _A(w\alpha )\right](a+1)\psi _A(w)\right\}`$
(36) $`v_A(w)`$ $`=`$ $`{\displaystyle \frac{2\zeta }{3(1+a+\beta )}}w{\displaystyle \frac{d\psi _A}{dw}},`$
where the subscript $`A`$ indicates the new unknowns. Inserting $`v_A(w)`$ in the first equation and neglecting a coefficient proportional to $`w^1`$, a single equation in $`\psi _A`$ is obtained

(37)
$$\frac{2\zeta ^2}{3(a+1+\beta )}\frac{d^2\psi _A}{dw^2}+a\psi _A(w+\alpha )+\psi _A(w\alpha )(a+1)\psi _A(w)=0.$$

This is a linear difference-differential equation. It is easy to see that $`\mathrm{exp}(\lambda w)`$ is a solution of eq. (37) if and only if $`\lambda `$ satisfies the transcendent equation

(38)
$$\frac{2\zeta ^2}{3(a+1+\beta )}\lambda ^2+ae^{\lambda \alpha }+e^{\lambda \alpha }(a+1)=0.$$

In Appendix A, we prove that this equation admits only two solutions: $`\lambda =0`$ and $`\overline{\lambda }(1,0)`$. We give also a simple approximation for $`\overline{\lambda }`$.
By means of this solution of eq. (37), we get the original quantities

(39)
$$n_A(w)=c\sqrt{w}e^{\overline{\lambda }w},v_A(w)=\frac{2\zeta w}{3(1+a+\beta )}\overline{\lambda }ce^{\overline{\lambda }w},$$

where $`c`$ is a positive constant, in order to satisfy the condition (31). Despite this solution is obtained for large values of $`w`$ it often agrees with the numerical results (which we present in next section) in the whole range of $`w`$.
In the limit case $`\zeta =0`$ (no electric field) $`\overline{\lambda }=1`$ and the solution of (37) gives the correct function $`N`$ obtained from the equilibrium maxwellian distribution function .

For high electric fields approximate solutions are known (see ,  for example). These are obtained by means of an asymptotic expansion of the original differential-difference equations, which are transformed in simple ordinary differential equations and are explicitly solved. Our solution, obtained by the asymptotic equation, is different with respect to the previous ones, because we do not approximate the difference terms $`f_0(ϵ\pm \mathrm{}\omega )`$ ($`\omega `$ is the constant optical phonon frequency), appearing in the equations, using Taylor formula. In fact the asymptotic equation remains of differential-difference type.

We are interested in calculating the hydrodynamical velocity of electrons

$$v_h:=\frac{{\displaystyle _{u_0}^+\mathrm{}}V_1(u)𝑑u}{{\displaystyle _{u_0}^+\mathrm{}}N(u)𝑑u}=\sqrt{\frac{k_BT_L}{m^{}}}\frac{{\displaystyle _0^+\mathrm{}}v(w)𝑑w}{{\displaystyle _0^+\mathrm{}}n(w)𝑑w}.$$

By a simple calculation, it is possible to verify that, using (39), the following value is obtained

$$\left|v_{A_h}\right|=\sqrt{\frac{k_BT_L}{m^{}}}\frac{4}{3(1+a+\beta )\sqrt{\pi }}\zeta \sqrt{|\overline{\lambda }|}.$$

If we use the approximate value (see Appendix A) of $`\overline{\lambda }`$, we obtain the saturation velocity

$$\underset{\zeta +\mathrm{}}{lim}\left|v_{A_h}\right|=\sqrt{\frac{k_BT_L}{m^{}}}\sqrt{\frac{8\alpha (a1)}{3\pi (1+a+\beta )}}.$$

8. Finite energy model
The difficulties to solve numerically eqs. (25)-(26) are given by the deviating arguments in the right side of equations and by the infinite interval for the variable $`w`$. The first suggests to use a difference scheme to discretize the equations with the constraint that the step $`\mathrm{\Delta }w`$ be such that $`{\displaystyle \frac{\alpha }{\mathrm{\Delta }w}}`$ is an integer. The second difficulty is solved by considering a finite interval for the variable $`w`$. This truncation must be made carefully. To explain this aspect of the problem, we recall that we are looking for a solution which describes the no-runaway phenomena. It is a particular configuration of the gas, where the applied electric field is balanced by the effect of the collisions with the background. In other words, the energy that an electron receives from the external electric field is transferred to the phonons so that a spatial homogeneous and static configuration is maintained.

In order to find this solution also in the case of a finite range for the energy $`w`$, we introduce a cut in the kernels K and $`\text{K}_0`$. Let us choose a positive integer $`\overline{N}`$. Put $`M=\overline{N}\mathrm{}\omega `$, we consider the interval $`[0,M]`$ and the function $`H(u)`$ which coincides with the Heaviside step function $`\theta (u)`$ except for $`u=0`$ where $`H(u)=0`$. Then the new kernels are defined by

$`\text{K}^M(\epsilon (𝐤),\epsilon (\stackrel{~}{𝐤}))`$ $`=`$ $`\text{K}(\epsilon (𝐤),\epsilon (\stackrel{~}{𝐤}))H(M\mathrm{max}\{\epsilon (𝐤),\epsilon (\stackrel{~}{𝐤})\}),`$
$`\text{K}_0^M(\epsilon (𝐤),\epsilon (\stackrel{~}{𝐤}))`$ $`=`$ $`\text{K}_0(\epsilon (𝐤),\epsilon (\stackrel{~}{𝐤}))H(M\mathrm{max}\{\epsilon (𝐤),\epsilon (\stackrel{~}{𝐤})\}).`$
These kernels imply that a particle has zero probability of collision in the following two cases:

* the particle before the collision has an energy less than $`M`$, but after it should have an energy equal or greater than $`M`$;
* the particle before the collision has an energy equal or greater than $`M`$.

Using the new kernels eqs. (25)-(26) are substituted by the following

(40) $`\zeta {\displaystyle \frac{dv}{dw}}=A(n)(w)`$
(41) $`\zeta \left[{\displaystyle \frac{d\text{ }}{dw}}(r(w)n)+q(w)n\right]=B(v)(w),`$
where
$`A(n)(w):=s(w)\left[aH(\alpha \overline{N}w\alpha )n(w+\alpha )+H(\alpha \overline{N}w)n(w\alpha )\right]`$
$`\left[H(\alpha \overline{N}w\alpha )s(w+\alpha )+aH(\alpha \overline{N}w)s(w\alpha )\right]n(w)`$
$`B(v)(w):=\left\{H(\alpha \overline{N}w\alpha )s(w+\alpha )+H(\alpha \overline{N}w)\left[as(w\alpha )+\beta s(w)\right]\right\}v(w).`$
The corresponding boundary conditions becomes

(42) $`n(0)=0,v(\alpha \overline{N})=0`$
(43) $`n(w)0\text{for every }w0\text{ and }{\displaystyle _0^{\alpha \overline{N}}}n(w)𝑑w>0.`$
The boundary condition on $`v`$ is the natural consequence of the condition (30). The choice of the above boundary conditions and the cut are suggested by a simple physical picture. No-runaway phenomena may be considered as the limit for $`t+\mathrm{}`$ of a suitable spatial homogeneous solution of eqs. (11)(21). For time-depending solutions, the new kernels and the condition $`v(\alpha \overline{N})=0`$ imply that a particle with initial energy less than $`M`$ will have, for all time, energy less than $`M`$. For the same reason particles having energy greater than $`M`$ will have, for all time, energy greater than $`M`$. Then, the total number of the particles having energy less than $`M`$ will be constant in time. Different choices can violate the conservation of the particle number for electrons having energy less than $`M`$. In these cases the number of particles having energy less than $`M`$ could tend to zero or to infinity.
9. Numerical solutions
In order to discretize eqs. (40)-(41), we fix the step $`\mathrm{\Delta }w`$ and consider in the interval $`[0,\alpha \overline{N}]`$ the points

$$w_i=i\times \mathrm{\Delta }w(i=0,1,2,\mathrm{}.,\text{N})$$

where $`\text{N}={\displaystyle \frac{\alpha \overline{N}}{\mathrm{\Delta }w}}`$. If we integrate eqs. (40)-(41) in the interval $`[w_i,w_{i+1}]`$ exactly where it is possible and using trapezoidal rule otherwise, we obtain the linear algebraic system

(44) $`{\displaystyle \frac{2\zeta }{\mathrm{\Delta }w}}(v_{i+1}v_i)+A_{i+1}+A_i=0(i=0,1,\mathrm{},\text{N}1),`$
(45) $`v_0=0`$
(46) $`{\displaystyle \frac{2\zeta }{\mathrm{\Delta }w}}(n_{i+1}r_{i+1}n_ir_i)\zeta (q_{i+1}n_{i+1}+q_in_i)B_{i+1}B_i=0`$
$`(i=1,2,\mathrm{},\text{N}1),`$
where, for any function $`\phi (w)`$, $`\phi _i=\phi (w_i)`$. The corresponding boundary conditions are

(47) $`n_0=0,v_\text{N}=0`$
(48) $`n_i0i=0,1,\mathrm{}\text{N}\text{and }{\displaystyle \frac{\mathrm{\Delta }w}{2}}{\displaystyle \underset{i=0}{\overset{\text{N}1}{}}}\left(n_{i+1}+n_i\right)>0.`$
Equations (S0.Ex39) are obtained considering only the interval $`[w_i,w_{i+1}]`$ for $`i1`$, because in the first point $`w_0`$ we have written the exact equation $`v_0=0`$ (i.e. eq. (45)), which derives from eq. (41) with the boundary condition $`n(0)=0`$. It is simple to deduce from (44) the following equations

(49)
$$v_j=v_\text{N}+\frac{\mathrm{\Delta }w}{2\zeta }\underset{i=j}{\overset{\text{N}1}{}}(A_i+A_{i+1})(j=0,1,\mathrm{},\text{N}1).$$

Moreover, as a consequence of the appropriate choice of the cut, the conservation of the particle number holds (see ref. ), being

(50)
$$\underset{i=0}{\overset{\text{N}1}{}}(A_i+A_{i+1})=0.$$

Eq. (50) is the corresponding of eq. (32) introducing the cut in the kernels and performing the integral using the trapezoidal rule.
Now, eq. (49) gives $`v_0=v_\text{N}`$ for $`j=0`$. Due to the boundary condition $`v_\text{N}=0`$, there follows that (45) is a consequence of eqs. (49), or equivalently eq. (49) for $`j=0`$ coincides with eq. (45). For numerical convenience, we prefer to eliminate the case $`j=0`$ in (49). Now, we must add the condition (48). To select an unique solution, we choose the equation $`n_1=1`$. This, with the conditions $`n_i0`$ for $`i=2,3,\mathrm{}\text{N}`$ guarantees (48). This makes the number of equations equal to the number of the unknowns. For the sake of clarity, we rewrite the complete linear system

(51) $`v_\text{N}=0,v_0=0,n_0=0,n_1=1,`$
(52) $`v_j=v_\text{N}+{\displaystyle \frac{\mathrm{\Delta }w}{2\zeta }}{\displaystyle \underset{i=j}{\overset{\text{N}1}{}}}(A_i+A_{i+1})(j=1,2,\mathrm{},\text{N}1)`$
(53) $`{\displaystyle \frac{2\zeta }{\mathrm{\Delta }w}}\left[n_{i+1}r_{i+1}n_ir_i\right]\zeta (q_{i+1}n_{i+1}+q_in_i)B_{i+1}B_i=0`$
$`(i=1,2,\mathrm{},\text{N}1).`$
10. Numerical results and conclusion
We are interested into solve eqs. (40)-(41) in the case of a silicon bulk device. The appropriate values for the parameters are given in the following table.

| $`m^{}=0.32m_e`$ | $`T_L=300`$ K | $`\mathrm{}\omega =0.063`$ eV |
| --- | --- | --- |
| $`\text{K}={\displaystyle \frac{\left(D_tK\right)^2}{8\pi ^2\rho \omega }}`$ | $`D_tK=11.4`$ eV $`\stackrel{}{\text{A}}`$<sup>-1</sup> | $`\rho =2330`$ Kg m<sup>-3</sup> |
| $`\text{K}_0={\displaystyle \frac{k_BT_L}{4\pi ^2\mathrm{}v_0^2\rho }}\mathrm{\Xi }_d^2`$ | $`\mathrm{\Xi }_d=9`$ eV | $`v_0=9040`$ m sec<sup>-1</sup>. |
| $`\alpha _k=0.5eV^1`$ |  |  |

Here, $`m_e`$ denotes the electron rest mass. Using these parameters, we get $`\alpha 2.437`$ and $`\beta 5.986`$. The values for the electric field are $`|𝐄|=10^3`$, $`10^4`$ and $`10^5`$ Vcm<sup>-1</sup>. Except for low electric field there is a great difference between P (parabolic) and NP (non-parabolic) band. In fig. 1 we plot the velocity versus the electric field, obtained using the asymptotic solutions (39) and values yielded by numerical integration of eqs. (40)-(41) both in P and NP cases. In the first case the value for the saturation velocity obtained is $`1.27\times 10^5m/sec`$. The second model does not show a limit value as in the Monte Carlo simulations . For the largest value of $`|𝐄|`$ the velocity is $`0.95\times 10^5m/sec`$ The experimental value is $`1\times 10^5m/sec`$. In order to compare numerical and asymptotic solutions, we have chosen the condition

(54)
$$_0^{\alpha \overline{N}}n_A(w)𝑑w=_0^{\alpha \overline{N}}n(w)𝑑w$$

where $`n_A(w)`$ indicates the asymptotic solution. The integral in (54) of the function $`n_A`$ is performed analytically and the second numerically.
In such way the constant $`c`$ in (39) is chosen. The meaning of eq. (54) is evident and coherent with the criterion of the choice of the cut. We have made many numerical experiments, varying the values of the electric field, the step $`\mathrm{\Delta }w`$ and the cut $`\overline{N}`$. In fig. 2 the quantity $`N(u)`$, both for P and NP approximation, is shown for a low value of the electric field. Nevertheless a small range of the energy $`04\mathrm{}\omega `$ is sufficient to contain almost all the electrons, the numerical results for P case agree with the asymptotic solution. This is due to the fact that for a small values of the electric field $`\overline{\lambda }1`$. The quantity $`V_1(u)`$ is shown in fig. 3. Here the asymptotic solution is given by the formula

$$v_A(w)=\frac{2\zeta \sqrt{w^3}c\overline{\lambda }e^{\overline{\lambda }w}}{3\left(\sqrt{w+\alpha }+a\sqrt{w\alpha }+\beta \sqrt{w}\right)}$$

This choice is not coherent with the expansion, but it allows a better agreement with the numerical solutions of P case and it puts in evidence a discontinuity in the derivative at the point $`w=\alpha u=\mathrm{}\omega `$, due to the term $`\sqrt{w\alpha }`$. This square root is also present in the complete set of equations (25)-(26). Then, such discontinuity is also expected in the solutions. A deeper analysis could be necessary to prove mathematically this fact. Figs. 4 and 5 shown an intermediate case where the agreement is poor.
When the electric field is high (figs. 6 and 7) very small variations are observable. This agreement was reasonable for large values of $`u`$. In the other cases this is due to the number of particles having low energy which is enough small with respect to that of the whole gas.
The last figures show a quite surprising result in the parabolic case. We have used the same value of the electric field $`E`$ as in fig. 6 and 7. Here, a small value of $`\overline{N}`$ is used. The density $`N`$ agrees with the asymptotic solutions. The quantity $`V_1`$ subject to the boundary condition $`V_1(\overline{N}\alpha )=0`$, follows the asymptotic solutions almost to the end of the interval.
Acknowledgments
We acknowledge partial support from Italian Consiglio Nazionale delle Ricerche (Prog. N. 96.03855.CT01) and TMR project n. ERBFMRCT970157 “Asymptotic Methods in Kinetic Theory”.
Appendix A
The function

$$\eta (\lambda )=\frac{2\zeta ^2}{3(a+1+\beta )}\lambda ^2+ae^{\lambda \alpha }+e^{\lambda \alpha }(a+1)$$

for $`\zeta 0`$ is convex, because $`\eta ^{\prime \prime }(\lambda )0`$ for every real $`\lambda `$.
Since $`\eta (0)=0`$ and $`\eta ^{}(0)=\alpha (a1)>0`$, then there follows that the minimum is negative and is reached at a point $`\lambda _m<0`$.
From $`\eta (1)=2\zeta ^2\left(3(a+1+\beta )\right)^1`$ (recall that $`a=e^\alpha `$) we get $`1<\lambda _m<0`$. Then, there exists only $`\overline{\lambda }0`$ such that $`\eta (\overline{\lambda })=0`$. Moreover $`\overline{\lambda }[1,0[`$. An approximation of this value is obtained expanding $`e^{\pm \lambda \alpha }`$ around $`\lambda =0`$ up the term of second order. A second order algebraic equations follows

$$\left[\frac{2\zeta ^2}{3(a+1+\beta )}+\frac{\alpha ^2}{2}(a+1)\right]\lambda ^2+\alpha (a1)\lambda =0.$$

Then

$$\overline{\lambda }\frac{6(a+1+\beta )\alpha (a1)}{4\zeta ^2+3(a+1+\beta )\alpha ^2(a+1)}.$$

For $`\zeta 1`$ we get

$$\overline{\lambda }\frac{3\alpha (a1)(a+1+\beta )}{2\zeta ^2}.$$