text:
# Cosmic Rays below 𝑍=30 in a diffusion model: new constraints on propagation parameters.

## 1 Introduction

Understanding the composition and spectral features of cosmic rays has always been an astrophysical challenge. On one hand, the observational data have long been scarce and suffered from large uncertainties. On the other hand, the theoretical predictions to which these data should be compared to have also suffered from several drawbacks. Composition and spectra arise from the nuclear interaction of an initial distribution of energetic particles with interstellar matter (spallations) and their electromagnetic interactions with galactic magnetic fields (acceleration and diffusive reacceleration). First, the nuclear cross sections to be used were not very well known until recently. Second, our knowledge of the galactic magnetic field is far from complete. Cosmic rays are sensitive to its scale inhomogeneities (diffusion) which are not well observed. They are also sensitive to the presence of plasma shock–waves (acceleration in localized sources and diffusive reacceleration). Third, composition and spectra are altered as the cosmic rays enter the solar magnetic field, so that some more modelling has to be done in order to infer interstellar spectra from observations.

However, despite these sources of uncertainty, some gross features of the cosmic ray properties are well established. First, the study of secondary–to–primary ratio shows that a cosmic ray has crossed an average of about $`x=9\text{g}\text{cm}^2`$ of interstellar matter between its acceleration and its detection. Second, the isotopic ratio of radioactive species shows that it took about $`\tau =20\text{Myr}`$ between the same events. As cosmic rays have a velocity close to $`c`$, we can infer the average density of the medium to be $`n=x/(\tau cm_H)0.3\text{cm}^3`$. As the disk density is about $`1\text{cm}^3`$, the cosmic rays must spend a fraction of time in an empty region, called the diffusion halo. Cosmic rays are produced and destroyed in the galactic disk, they diffuse in a larger zone whose characteristics are not well known, and eventually they can escape from this zone. To have a good modelling of cosmic ray propagation, we should know the geometry and size of the diffusion zone, the characteristics of the galactic magnetic field, and the sources.

In this paper, we use the existing data on nuclei, gathered by balloon–borne and space experiments during the last thirty years, to put constraints on the parameters describing the propagation of cosmic ray nuclei. One consequence is that many of the large uncertainties affecting the calculation of the antiproton flux could be strongly reduced (F. Donato & al., in preparation la donato ). This is of utmost relevance for the study of exotic sources of antiprotons and antideuterons.

## 2 The diffusion model

It has been recognized for a long time that the relevant physical propagation model to be used is the diffusion model (Berezinskii & al. 1990; Ginzburg & Syrovatskii 1964), though the so–called leaky box model has been widely preferred for decades because of its simplicity.

The steady–state differential density $`N^j(E,\stackrel{}{r})`$ of the nucleus $`j`$ as a function of energy $`E`$ and position $`\stackrel{}{r}`$ in the Galaxy, is given by (see for instance Berezinskii & al. 1990)

$`.(K^jN^jV_cN^j){\displaystyle \frac{}{E}}\left({\displaystyle \frac{.V_c}{3}}E_k\left({\displaystyle \frac{2m+E_k}{m+E_k}}\right)N^j\right)+{\displaystyle \frac{}{E}}(b^jN^j){\displaystyle \frac{1}{2}}{\displaystyle \frac{^2}{E^2}}(d^jN^j)+\stackrel{~}{\mathrm{\Gamma }}^jN^j=q^j+{\displaystyle \underset{m_k>m_j}{}}\stackrel{~}{\mathrm{\Gamma }}^{kj}N^k`$ (1)
The first terms represent diffusion ($`K^j`$ is the diffusion coefficient) and convection ($`V_c`$ is the convection velocity). The divergence of this velocity, expressed in the next term gives rise to an energy loss term connected with the adiabatic expansion of cosmic rays. Further, we have to take into account ionization and coulombian losses (the energy loss coefficients are specified below), plus a reacceleration term in first order derivative (all included in $`b^j`$) and finally a second order derivative in $`E`$ for the associated second order term in reacceleration ($`d^j`$ is the energy diffusion coefficient). These stand for the continuous losses. The last term of the l.h.s. takes care of the disappearance of the nucleus $`j`$ ($`\stackrel{~}{\mathrm{\Gamma }}^j`$ for short) due to its collisions with interstellar matter (ism). In the r.h.s., the source term $`q^j`$ takes into account the primary production and acceleration of nuclei described by an injection spectrum (for the sake of clarity, we have not written down the terms describing the contribution of radioactive species). Finally, the last term is for the secondary $`j`$ sources, namely spallation contribution111We use here the so–called straight–ahead approximation which implies the conservation of kinetic energy per nucleon during a spallation process
$$_0^{\mathrm{}}n_Hv^{}N^k(T^{})\sigma ^{kj}(T,T^{})𝑑T^{}=_0^{\mathrm{}}n_Hv^{}N^k(T^{})\sigma ^{kj}(T)\delta (TT^{})𝑑T^{}=n_HvN^k(T)\sigma ^{kj}(T)$$
(2)   $`\stackrel{~}{\mathrm{\Gamma }}^{kj}`$ from all other heavier nuclei. We use here the notation $`E`$ for total energy and $`E_k`$ for kinetic energy.

One needs to solve a complete triangular–like set of coupled equations since only heaviest nuclei contribute to a given nucleus. Quantities in this equation are functions of spatial coordinates (not time, steady–state being assumed) and of energy. Unless stated otherwise, the word energy will stand for kinetic energy per nucleon ($`T`$), since this is the appropriate parameter to be used, as it is conserved in spallation reactions (see footnote 2).

The other very popular model is the leaky box, in which all the quantities are spatially averaged (so that convection has no meaning), and the diffusion term is then replaced by an escape term

$$.(K^jN^j)\overline{N^j}/\tau _{esc}$$

which has the meaning of a residence time $`\tau _{esc}`$ (Myr) in the confinement volume. Note that in this case, the mean time $`\tau `$ spent in the box and the mean column density $`x`$ (in units of g$`\text{cm}^2`$) of crossed matter are treated on the same footing since the model is homogeneous. This is not true for multi–zone models and especially for the diffusion model where cosmic rays spend most of their time in the diffusion halo which has a zero matter density. It may seem very surprising that this simplified model is able to reproduce many observations on stable nuclei. Actually, it can be shown that leaky box models are often equivalent to diffusion models, as far as stable species are concerned (for a generic discussion on the “leakage–lifetime” approximation in cosmic ray diffusion, see Jones 70). This is not obvious from the equations, but it can be seen more readily if we solve them in a formalism called the weighted slab technique. Basically, it consists in writing a general solution under the form (see for example Berezinskii et al. 1990; Ginzburg & Syrovatskii 1964)

$$N^j(𝐫)=_0^{\mathrm{}}\stackrel{~}{N}^j(x)G(𝐫,x)𝑑x$$
(3)
Thus, (1) decouples into two independent equations. The first one involves $`G(𝐫,x)`$. It depends only on the geometry of the problem and on the chosen diffusion scheme (or the leaky box parameters), but not on the particular species $`j`$. The other one involves $`\stackrel{~}{N}^j(x)`$ and contains the nuclear physics aspects of the propagation. Different models (diffusion or leaky box) correspond to different $`G(𝐫,x)`$, but the equations on $`\stackrel{~}{N}^j(x)`$ are the same. $`G`$ is the path length distribution and is interpreted as the distribution of probability that a nucleus $`j`$ goes through a column density $`x`$ before reaching Earth. The global solution is then obtained by a convolution of the surviving quantity $`\stackrel{~}{N}^j(x)`$ under nuclear processes after $`x\text{g}\text{cm}^2`$, with the weighted probability associated.

In a leaky box model, it can be shown that the normalized path length distribution function is given by

$$G(x)=\frac{1}{\lambda _{\text{esc}}}\mathrm{exp}\left(\frac{x}{\lambda _{\text{esc}}}\right)$$
(4)
where the average quantity of matter crossed by a cosmic ray is given by

$$x=_0^{\mathrm{}}xG(x)𝑑x=\lambda _{\text{esc}}$$
(5)
In the diffusion model, it has been shown (Owens 1976) that for a wide class of geometries, the function $`G`$ is given by an infinite series of exponentials involving the diffusion coefficient. Moreover, the physical parameters are such that for sufficiently large grammage $`x`$, only the first exponential contributes so that we recover a leaky box model (Jones 1970; Ginzburg, Khazan & Ptuskin 1980; Schlickeiser & Lerche 1985; Lerche & Schlickeiser 1985). For the model used in this paper, it can be shown that $`\lambda _{\text{esc}}=\overline{m}nvhL/K`$ in the limit $`LR`$. However, this is no longer true for lower $`x`$, and then the leaky box model fails. This may explain, in conjunction with poorly known cross sections for heavy nuclei, why for many years cosmic ray physicists had to use modifications (Lezniak & Webber 1979; Garcia–Munoz & al. 1987; Webber & al. 1998b) of the leaky box path length distribution at small grammage to reconcile observations of B/C with those of sub–Fe/Fe. Even worst, this approximate equivalence between leaky box and diffusion model breaks down for radioactive species (Prishchep & Ptuskin 1975).

This motivates us to choose the diffusion model and solve the corresponding equation (1). From a more general point of view, note also that the weighted slab technique mentionned above should be used carefully. Actually, it was recently realized (Ptuskin, Jones & Ormes 1996) that pre–1996 estimations of cosmic ray parameters with this technique were not exact. First, at fixed energy per nucleon, the rigidity and therefore the propagation features depend on the nuclear species, whatever the energy. This prevents the above–mentionned decoupling. Second, energy losses that strongly depend on $`A`$ and $`Z`$ also take place at low energy. Ptuskin et al. 1996 presented a modification of this formalism which is correct in both the highly and non–relativistic regimes. These fundamental findings were confirmed by a more complete numerical treatment (Stephens & Streitmatter 1998) and give rise to differences of as much as 20%.

For the nuclear part, we adopted the direct approach of the shower technique. This means that the flux is first evaluated for the heavier primary cosmic ray, for which the diffusion equation is simple and does not couple to any other species. Then, the flux of the next nucleus is computed, with a spallation term depending only on the heavier nucleus, whose flux is known from the previous step. This procedure is repeated for all the nuclei down to the lightest one. We started at $`Z=30`$, the heavier non negligeable primary found in cosmic radiation (Binns & al. 1989), since $`(Z>30)10^3`$ Fe. For all light nuclei, we checked that it is sufficient to start from $`Z=16`$ (S). For the detail of solutions, the reader is referred to Appendix A. The next section is devoted to the inputs of the model.

## 3 The parameters of the model

### 3.1 Geometry of the Galaxy

The geometry of the problem used here is a classical cylindrical box (see for example Webber et al. 1992) whose radial extension is $`R`$, with a disk of thickness $`2h`$ and a halo of half–height $`L`$. Only $`L`$ is not fixed (see $`\mathrm{\S }`$5) and $`R=20\text{kpc}`$, $`h=100\text{pc}`$. Sources and interactions with matter are confined to the thin disk and diffusion which occurs throughout disc and halo with the same strength is independent of space coordinates. The solar system is located in the galactic disc ($`z=0`$) and at a centrogalactic distance $`R_{}=8\text{kpc}`$ (Stanek & Garnavich 1998; Alves 2000).

### 3.2 Cross sections

Cross sections play a crucial role in propagating the cosmic rays throughout the Galaxy. We have to distinguish bewtween the total inelastic cross section and the spallation cross section.

#### 3.2.1 Total inelastic cross section

The total inelastic cross section, which is actually a reaction cross section, is defined by

$$\sigma _{inel}^{tot}=\sigma ^{tot}\sigma _{elastic}^{tot}$$
(6)
Various empirical modifications of the original (Bradt & Peters 1950) geometrical approach have been used to parameterize this total cross section. Letaw, Silberberg & Tsao (1983) and Silberberg & Tsao (1990) produced the basic energy dependent equations further refined in Sihver & al. (1993) and in Wellish & Axen (1996). Recently Tripathi, Cucinotta & Wilson (1997a, 1997b, 1999) proposed a universal parameterization which is valid for all nucleus–nucleus reactions, including neutron induced reaction. All of these approaches have been compared in Silberberg, Tsao & Barghouty (1998), where the preference were given to Tripathi et al. (1997a, 1997b, 1999). For this reason, we adopted the latter parameterization throughout this study. The precision can be estimated to be about $`5\%`$.

#### 3.2.2 Spallation cross section

The fragmentation cross sections are somewhat more complicated. Basically, we can separate three distinct approaches : (i) microscopical description (Bondorf & al. 1995; Ramsey & al. 1998) (ii) semi–empirical formulae (Silberberg et al. 1998; Tsao, Silberberg & Barghouty, 1998; Tsao et al. 1999) and (iii) empirical formulae (Webber, Kish & Schrier 1990d). The first one takes into account the fundamental nuclear physics but is very computer time consuming and has not reached a precision as good as the two other approaches. The second method takes into account more phenomenological inputs, and is better designed for the unmeasured fragmentation cross sections which correspond mainly to the ultra heavy ($`Z30`$) nuclei. The third one breaks down in this region, but can be used for $`3Z30`$ because it is adapted to fit measurements. Consequently, we used the code of Webber et al. (1990d) available on the web 222http://spdsch.phys.lsu.edu/SPDSCH\_Pages/Software\_Pages/Cross\_Section\_Calcs/CrossSectionCalcs.html updated with new parameters (table V of Webber & al. 1998a). This parameterization takes advantage of the fact that about $`98\%`$ of all the reaction involved in cosmic ray propagation have now at least one point of energy measured. This set of formulae gives a precision of about $`10\%`$ but can not be further refined (Webber et al. 1998a). The code is extended for spallation on He with the parameterization given in Ferrando et al. (1988)333The parameterization given in this paper should be modified at low energy and for large $`\mathrm{\Delta }Z`$ (P. Ferrando 2001, private communication)..

To go further than this treatment, one should keep in mind the following points. First, the process of single or double nucleon removal is an important channel for the production of numerous nuclei. It has been parameterized in detail in Norbury & Townsend (1993) and in Norbury & Mueller (1994), although in the following we will consider Webber’s approach sufficient for our purpose. Second, most of the parameterizations make use of the straight–ahead approximation (see footnote 2). Tsao & al. (1995) relaxed this condition and find a 5% effect around the $`1\text{GeV}`$ peak of B/C (note that this effect is negligable for primary to primary and secondary to secondary ratios). Finally, it was found that all the measurements are not always consistent and systematic effects are likely to be present (Taddeucci & al. 1997; Vonach & al. 1997; Zeitlin & al. 1997; Flesh & al. 1999; Korejwo & al. 1999).

To summarize, we consider that Webber’s code is very well suited for modelling spallation cross sections, thanks to the many experiments developed these last years (Webber, Kish & Schrier 1990a, 1990b, 1990c; Webber et al. 1998b, 1998c; Chen et al. 1997). The overall precision can be roughly estimated to be better than 10% on H, and better than 20% on He.

### 3.3 Nuclear parameters

Apart from the stable species, we can distinguish two categories of nuclei to be propagated, depending on whether they live long enough to be observed or they decay quickly into other nuclei.

#### 3.3.1 Ghost nuclei

Ghost nuclei are defined as those having a half–lifetime larger than $`1\text{ms}`$ (so that they are detected in cross section measurements) and less than few hundred years (so that we can consider that disintegration is fast on a propagation timescale). Several reaction channels may yield a given nucleus. This means that for this nucleus, we have to add to the direct production the contributions of all the short lived parent nuclei

$$\sigma _{kj}^{effective}=\sigma _{kj}+\underset{nucleiX}{\overset{\text{ghost}}{}}\sigma _{kX}r(Xj)$$
(7)
where $`r(Xj)`$ is the branching ratio of the corresponding channel444To give an example, such a reaction reads for $`{}_{4}{}^{9}Be`$
$$\sigma _{k_4^9Be}^{effective}=\sigma _{i_4^9Be}+49.2\%\sigma _{k_3^9Li}+4.1\%\sigma _{k_3^{11}Li}$$
(8)  . A table of all the disintegration chains can be found in Letaw & al. 1984. To see the effect of the recent discovery of new intermediate short–lived nuclei and of the changes in branching ratio, we reconstructed from a recent compilation of nuclear properties (Audi et al. 1997) all the ghost nuclei and compared them with the Letaw table. We found that the net effect is negligible. This is because the discovered radioelements are very far from the valley of stability and their production is suppressed (as it can be observed experimentally). Such reaction chains should accordingly be taken as a fixed input for the future.

#### 3.3.2 Radioactive nuclei

Radioactive nuclei can be separated into three groups : unstable under $`\beta `$ decay, unstable under electronic capture and unstable under both reactions. All the physical effects related to these nuclei will be analysed in detail in a separate paper (F. Donato & al., in preparation la donato2 ).

### 3.4 Diffusion coefficient

Physically, diffusion arises because charged particles interact with the galactic magnetic field inhomogeneities. This is an energy–dependent process because higher energy particles are sensitive to larger spatial scales. As a result, the diffusion coefficient is related to the power spectrum of these inhomogeneities, which is poorly known. Several analytical forms of this energy dependence have been assumed in the literature. In particular, leaky box models give a purely phenomenological form

$$\{\begin{array}{cc}\lambda _{\text{esc}}=\lambda _0\beta \hfill & \text{if }<_c\text{in GV}\hfill \\ \lambda _{\text{esc}}=\lambda _0\beta \left(/_c\right)^\delta \hfill & \text{if }>_c\text{in GV}\hfill \end{array}$$

where $`=p/Z`$ stands for the particle rigidity. A more fundamental origin may be found for the previous expression, as one can show in a one–dimensional diffusion model with convection (Jones & al. 2001). All the effects that led to that phenomenological form are explicitely taken into account in our treatment and we may therefore consider a more natural form for $`K(E)`$, inferred from magnetohydrodynamics considerations (Ptuskin et al. 1997):

$$K(E)=K_0\beta \times ^\delta $$
(9)
where the normalisation $`K_0`$ is expressed in kpc$`{}_{}{}^{2}\text{Myr}_{}^{1}`$. Nevertheless, we have also tested leaky box inspired forms, but these never give good fits to the data (see next section).

### 3.5 Sources

A charged particle is called a cosmic ray when it has been accelerated to an energy greater than $`100\text{MeV/n}`$ as far as our study is concerned. Several physical processes may be responsible for this energy gain; they all involve an interaction of the particle with a plasma shockwave front. This is responsible for the power–law energy dependence (see below). The term source refers to the place where the charged particles have been promoted to cosmic rays. We present here the various hypotheses about their spatial and spectral distributions and their compositions.

#### 3.5.1 Spatial distribution

First, as explained above, we assumed that the Galaxy has a cylindrical symmetry and that the sources are located in the galactic disc. So we only have to specify a radial distribution. We have used the Case & Bhattacharya (1996, 1998) distributions, but the results presented below are not sensitive to this particular choice.

#### 3.5.2 Spectral distribution

The source spectral distributions can be split in two terms

$$q^j(E)=q_0^jQ^j(E)$$
(10)
where $`q_0^j`$ represent the composition and $`Q^j(E)`$ the energy dependence. Several spectra have been used in the past (Engelmann & al. 1985),

$`Q^j(E)`$ $``$ $`p^{\alpha _j}`$
$`Q^j(E)`$ $``$ $`\beta E_{kin}^{\alpha _j}`$ (11)
$`Q^j(E)`$ $``$ $`E_{tot}^{\alpha _j}`$
where $`p`$ is the momentum and where the spectral index $`\alpha _j`$ may depend on the species $`j`$ considered. These forms are equivalent at high energies but we focused on the first one only (equivalent to a power–law in rigidity) as the others do not yield good fits to the data. The sum $`\delta +\alpha _j`$ gives the overall power–law index of each species at high energy ($`1\text{TeV}/n`$): this is because the measured flux at these energy is proportional to $`Q^j(E)/K(E)^{\delta +\alpha _j}`$. In other words, at sufficiently high energy, energetic redistribution and spallation reactions do not affect propagation. Consequently, measured spectral indexes enter as a fixed input, and we took them in Wiebel–Sooth, Biermann & Meyer (1998). For consistency, we required that whatever $`\delta `$, the parameters $`\alpha _j`$ are adjusted so as to reproduce the above–mentionned values. This ensures that the spectra we compute agree with data at high energy.

#### 3.5.3 Composition

The composition $`q_0^j`$ is usually inferred from a comparison between leaky box calculations and data. Alternatively, a more fundamental approach can be used. It may be taken either as the solar system composition convoluted with the fip (First Ionisation Potential), or convoluted with the volatility. As pointed out in Webber (1997), the two alternatives are not so different. In this work, the isotopic composition of each element was taken from Anders & Grevesse (1989) and Grevesse & Sauval (1998), and the $`q_0`$ dependence of the fip has been taken from Binns et al. (1989). Then the primary cosmic ray elemental composition has been adjusted to fit the heao–3 data at $`10.6\text{GeV/n}`$, while keeping constant the relative isotopic abundances.

### 3.6 Galactic wind, energy loss, reacceleration

Once we have chosen a generic model, we have to include effects which affect the spectra at low energies. A first step is the inclusion of energetic losses. Then a most complete treatment must take into account reacceleration and convection.

#### 3.6.1 Energy losses

There are two types of energy losses which are relevant for nuclei: ionization losses in the ism neutral matter ($`90`$% H and $`10`$% He), and Coulomb energy losses in a completely ionized plasma, dominated by scattering off the thermal electrons. We will use $`<n_e>0.033\text{cm}^3`$, and $`T_e10^4\text{K}`$ (Nordgren 1992). Complete formulae for these two effects are compiled, for example, in Strong & Moskalenko (1998) or in Mannheim & Schlickeiser (1994). There is also a third type of loss related to the existence of a galactic wind. We recall that in equation (1), there is a term proportional to $`.V_c`$. As we choose a constant wind velocity in the $`z`$ direction (see farther), we could conclude at a first glance that this term vanishes. Actually, this is the case everywhere in the Galaxy except at $`z=0`$ where a discontinuity occurs, due to the opposite sign of the wind velocity above and below the galactic plane. One gets a term that can be expressed in the same form as ionization and coulombian losses with an effective term

$$\frac{dE}{dt}_{Adiab}=E_k\left(\frac{2m+E_k}{m+E_k}\right)\frac{V_c}{3h}$$
(12)
$`E_k`$ stands for the total kinetic energy and it should not be confused with the kinetic energy per nucleus frequently used in this paper. Finally, because of these effects the propagation equation must be solved numerically.

#### 3.6.2 Reacceleration

As the fractional energy changes in a single collision are small, we can treat Fermi acceleration (Fermi 1949; 1954) using the Fokker–Planck formalism (Blandford & Eichler 1987). In one dimension

$`{\displaystyle \frac{f}{t}}+v{\displaystyle \frac{f}{x}}`$ $`=`$ $`{\displaystyle \frac{}{p}}\left[{\displaystyle \frac{\mathrm{\Delta }p}{\mathrm{\Delta }t}}f+{\displaystyle \frac{1}{2}}{\displaystyle \frac{}{p}}\left({\displaystyle \frac{(\mathrm{\Delta }p)^2}{\mathrm{\Delta }t}}f\right)\right]`$ (13)
where $`\mathrm{\Delta }p/\mathrm{\Delta }t`$ and $`(\mathrm{\Delta }p)^2/\mathrm{\Delta }t`$ are expressed in terms of the probability for changing the momentum $`p`$ to $`p+\mathrm{\Delta }p`$ during time $`\mathrm{\Delta }t`$. A simplification occurs when the recoil of the scatterer can be ignored. In general, the principle of detailed balance ensures that the probability for a gain is equal to that for a loss. This implies

$$\frac{\mathrm{\Delta }p}{\mathrm{\Delta }t}=\frac{1}{2}\frac{}{p}\frac{(\mathrm{\Delta }p)^2}{\mathrm{\Delta }t}$$
(14)
which reduces the equation into

$`{\displaystyle \frac{f}{t}}+v{\displaystyle \frac{f}{x}}={\displaystyle \frac{}{p}}K_{pp}{\displaystyle \frac{}{p}}\left(f\right)`$ (15)
where

$$K_{pp}\frac{1}{2}\frac{(\mathrm{\Delta }p)^2}{\mathrm{\Delta }t}$$
(16)
From another point of view, we can start with a collisionless Boltzmann equation and expand it up to second order in perturbed quantities (magnetic irregularities). It gives the full diffusion equation in quasi–linear theory (Schlickeiser 1986; Kulsrud & Pearce 1969). Assuming a constant flow of magnetic irregularities (corresponding to a constant galactic wind) in the $`z`$ direction, we find

$$\frac{f}{t}=V_cf+(Kf)+\frac{1}{p^2}\frac{}{p}\left[K_{pp}p^2\frac{f}{p}\right]+Q(r,p,t)$$
(17)
which is the propagation equation we used in Appendix A (up to several nuclear terms, cf eq. (A1)). In this theory, the momentum diffusion coefficient $`K_{pp}`$ is naturally connected to the space diffusion coefficient $`K`$ by the relation

$$K_{pp}=\frac{V_{a}^{}{}_{}{}^{2}}{9K}p^2$$
(18)
where $`V_a`$ is the Alfvén speed of the scattering centers.

We emphasize that expressions (15) and (17) are just simplifications of the full Fokker–Planck equations. As pointed out by Ostrowski & Siemieniec-Oziȩblo (1997), if this simplification is natural in the kinetic description of a gas of scattering particles, this is no longer the case for particles scattered off external ‘heavy’ scattering centers. Nevertheless, this approximation which arises in a number of astrophysical cases is sufficient for our purpose. We thus have to treat a simple energy diffusion coefficient, which we evaluated in the no–recoil hard sphere scattering centers approximation as

$$K_{EE}\frac{1}{2}\frac{(\mathrm{\Delta }E)^2}{\mathrm{\Delta }t}=\frac{2}{3}V_{a}^{}{}_{}{}^{2}\frac{\sigma _Kn_K}{c}E^2\beta ^3$$
(19)
where $`\sigma _Kn_K`$ describes the rate of collisions with scatterers which is related to the mean free path by $`\lambda _K=1/\sigma _Kn_K`$. As $`K(E)=(1/3)\lambda _Kv`$, the previous relation gives the final coefficient to be used in our calculation

$$K_{EE}=\frac{2}{9}V_{a}^{}{}_{}{}^{2}\frac{E^2\beta ^4}{K(E)}$$
(20)
Our approch is basically similar to that of Heinbach & Simon (1995). The meaning of the Alfvénic speed derived in their analysis is not very clear, since the numerical factor in front of expression (20) may vary according to the hypotheses made for the scattering process. The assumption that reacceleration occurs only in the thin disc is supported by recent complete magnetohydrodynamics simulations (Ptuskin et al. 1997) describing the evolution of the hot gas, cosmic rays and magnetic field.

#### 3.6.3 Integration into equation

We define

$$b_{loss}^j(E)=\frac{dE}{dt}_{Ion}+\frac{dE}{dt}_{Coul}+\frac{dE}{dt}_{Adiab}$$
(21)
and with $`K_{EE}^j(E)`$ given by (20), one obtains (see Appendix)

$$A_i^jN_i^j(0)=\overline{𝒬^j}2h\frac{}{E}\left\{b_{loss}^j(E)N_i^j(0)K_{EE}^j(E)\frac{}{E}N_i^j(0)\right\}$$
(22)
This is a second order equation which is solved for each order $`i`$ of the Bessel decomposition. $`A_i^j`$ is given by (A10) and $`\overline{𝒬^j}`$ (A5) is a generic source term (see Appendix A for details). Relation (22) has been developed on our energy array so as to be transformed into a matrix equation that is directly inverted. Two boundary conditions need to be implemented. The effects of energy losses and diffusive reacceleration have been assumed to be negligible at the high energy tip ($`100\text{GeV/n}`$) of the spectrum. We have furthermore imposed a vanishing curvature ($`^2N_i^j/E^2=0`$) at the lowest energy point.

#### 3.6.4 Galactic wind

It has been recognized for a long time that a thin disc configuration would be disrupted by cosmic ray pressure (Parker 1965, 1966; Ko, Dougherty & McKenzie 1991). It can be stabilized by the presence of a halo, but further considerations imply that this halo would not be static either (presence of a convective wind). Consequences of a wind has been firstly investigated by Ipavitch (1975) and since then it has been observed in other galaxies (Lerche & Schlickeiser 1982a, 1982b; Reich & Reich 1988; Werner 1988). In our own, its effects on nuclei have been investigated in various models (Jones 1979; Kóta & Owens 1980; Freedman & al. 1980; Bloemen & al. 1993), and in this work we adopted a very simple and tractable form for $`V_c`$ throughout the diffusive volume ($`dV_c/dz=0`$), following Webber, Lee & Gupta (1992). A constant wind with free escape boundary is a reasonable approximation to a more sophisticated variable–wind model (see also conclusions from magnetohydrodynamics simulations (Ptuskin et al. 1997)).

#### 3.6.5 Comparison with other works

Our semi–analytical two–dimensional model lies, in some sense, between the complete numerical resolution of Strong & Moskalenko (1998) and the recent study of Jones et al. (2001). As compared to the first one, we have a simplified description of the matter distribution, for which the full analysis of the parameter space does not require too much cpu time. The second approach can be considered as the limiting case for our two–dimension model when the radial extension is much larger than the halo size ($`RL`$). Note that a similar analysis was performed by Letaw, Silberberg & Tsao (1993) in the framework of a leaky box model or Seo & Ptuskin (1994) in a one–dimensional diffusion model but their less systematic approach does not allow a full exploration of parameter space.

## 4 Data Analysis

For the aim of our analysis, we can consider different classes of flux ratios: primary–to–primary (e.g. C/O), secondary–to–primary (e.g. B/C or sub–Fe/Fe), secondary–to–secondary (e.g. Li/B or Be/B), ratio of isotopes either stable (e.g. $`{}_{}{}^{10}\mathrm{B}/^{11}\mathrm{B}`$) or unstable (e.g. $`{}_{}{}^{10}\mathrm{Be}/^9\mathrm{Be}`$). Each of these may be an indicator of some dominant physical phenomenon and be particularly sensitive to the corresponding diffusion parameters. The ratio of two primaries is practically insensitive to changes in all the parameters, since they have the same origin and undergo the same physical processes. So, it can be a very useful tool to fix their source abundances. We will return to this subject. Similar conclusions, even if less strong, may be drawn for the ratio of two isotopes of the same species, such as $`{}_{}{}^{10}\mathrm{B}/^{11}\mathrm{B}`$. Indeed, at very low energy values this quantity is slightly affected by changes in the injection spectra, but the effect is too weak to constrain free parameters.

The most sensitive quantity is B/C, as B is purely secondary and its main progenitor C is primary. The shape of this ratio is seriously modified by changes in the diffusion coefficient, in the height of the diffusion halo, in the convective and alfvénic velocities. Moreover, it is also the quantity measured with the best accuracy, so that it is ideal to test models. Indeed, being the ratio of two nuclei having similar Z, it is less sensitive to systematic errors than single fluxes or other ratios of nuclei with more distant Z. For the same reasons, the sub–Fe/Fe may also be useful. Unfortunately, since existing data are still affected by sizeable experimental errors, we can only use them to cross–check the validity of B/C but not to further constrain the parameters under scrutiny. Another particularly interesting quantity is the ratio $`{}_{}{}^{10}\mathrm{Be}/^9\mathrm{Be}`$. Since $`{}_{}{}^{10}\mathrm{Be}`$ is a radioactive element, it is very sensitive to the processes which can occur in the halo. Therefore, we could use this ratio in particular to bracket the size of the halo (F. Donato & al., in preparation).

### 4.1 Solar modulation

We have estimated the effect of the solar wind on the energies and intensities of cosmic–rays following the prescriptions of the force field approximation (Perko 1987). The modification in the total interstellar energies of a nucleus with charge $`Z`$ and atomic number $`A`$, corresponds to the shift $`Z\varphi `$ (we set the electronic charge equal to unity):

$$E^{\mathrm{TOA}}/A=E^{\mathrm{IS}}/A\left|Z\right|\varphi /A.$$
(23)
Here $`E^{\mathrm{TOA}}`$ and $`E^{\mathrm{IS}}`$ correspond to the top–of–atmosphere (modulated) and interstellar total energy, respectively. The solar modulation parameter $`\varphi `$ has the dimension of a rigidity (or an electric potential), and its value varies according the 11–years solar cycle, being greater for a period of maximal solar activity. Often people refer to the equivalent quantity $`\mathrm{\Phi }=\left|Z\right|\varphi /A\frac{1}{2}\varphi `$. Once the momenta at the Earth $`p^{\mathrm{TOA}}`$ and at the boundaries of the heliosphere $`p^{\mathrm{IS}}`$ are determined, the interstellar flux of the considered nucleus is related to the toa flux according to the simple rule:

$$\frac{\mathrm{\Phi }^{\mathrm{TOA}}\left(E^{\mathrm{TOA}}\right)}{\mathrm{\Phi }^{\mathrm{IS}}\left(E^{\mathrm{IS}}\right)}=\left\{\frac{p^{\mathrm{TOA}}}{p^{\mathrm{IS}}}\right\}^2$$
(24)
The determination of the modulation parameter $`\mathrm{\Phi }`$ is totally phenomenological and suffers from some uncertainty. As explained in the following, we will deal with data taken around period of minimal solar activity, for which we fixed $`\mathrm{\Phi }=250\text{MV}`$.

### 4.2 The dataset we used

In order to test our diffusion model we have widely employed the data taken by the experiment which was launched on board the nasa heao–3 satellite (Engelmann et al. 1990). The relative abundances of elements, with charge from 4 to 28 and for energy ranging from 0.6 to 35 GeV/n, have been measured with unprecedent accuracy. Data have been taken in 1979–80, around a solar minimal activity. For the case of B/C, heao–3 quoted errors are 2–3$`\%`$. To perform the analysis, we also included data from balloons (Dwyer & Meyer 1987) and from the isee–3 experiment (Krombel & Wiedenbeck 1988), even if the relevant error bars are wider. The first one collected data in 1973–75 for energies spanning from around 1.7 to 7 GeV/n. The second experiment was operating during 1979–81 on board a spacecraft, in the energy range 100–200 MeV/n. In some of the figures presented below, we also plot – for purely illustrative goals – the data point from imp–8 (Garcia–Munoz et al. 1987) and the voyager experiments (Lukasiak, McDonald & Webber 1999) (we did not include Ulysse data point (DuVernois & al. 1996), since it corresponds to a period of maximal solar activity). They were not included in the $`\chi ^2`$ analysis since they are the most sensitive to solar modulation. Nevertheless, we have checked that the best $`\chi ^2`$ values were not significantly modified when these points were added.

As regards the sub–Fe(Sc+V+Ti)/Fe ratio, we used data from heao–3 (Engelmann et al. 1990) and from balloons (Dwyer & Meyer 1987). In both cases the error bars, around 10$`\%`$, are significantly larger than for B/C.

## 5 Results

We varied the relevant parameters $`K_0`$, $`L`$, $`V_c`$, $`V_a`$, and $`\delta `$ of the diffusion model described above. The constraints are much simpler to express with the combinations $`L`$, $`K_0/L`$, $`V_c`$ and $`V_a/\sqrt{K_0}`$ (this last expression appears naturally in $`K_{EE}(E)`$ describing reacceleration). As for the primaries, we have adjusted the source abundance for nitrogen which is not a pure secondary. We calculated the $`\chi ^2`$ over the 26 experimental points listed in 4.2, for each possible combination obtained varying the free parameters in the whole parameter space. In the following, we will focus on B/C and we cross–check the compatibility of the obtained parameters when applied to other fluxes.

### 5.1 Results of B/C and sub–Fe/Fe analysis for $`\delta =0.6`$

In a first step, the power $`\delta `$ of rigidity in the expression for the diffusion coefficient (see Eq. (9)) has been fixed to 0.6. This will be the reference value for the bulk of our analysis. We only show models giving a $`\chi ^2`$ less than 40. The best models have $`\chi ^228`$. Some of these configurations have a very small halo size in which case the condition $`hL`$ is no longer valid. Thus we also required $`L>1\text{kpc}`$ in the whole analysis.

The first result is that only the source spectrum $`q^j(E)=q_0^jp^{\alpha _j}`$ gives a good fit to the data. It will be used all the way long. For each $`L`$ and $`K_0/L`$, we vary $`V_a/\sqrt{K_0}`$ and $`V_c`$ and the best $`\chi ^2`$ are shown in Fig. 1. All the good models are located in a narrow strip in the plane $`L`$$`K_0/L`$, showing a strong correlation between these two parameters. However, there is still a degeneracy. It should be lifted with further specific analysis of radioactive nuclei (F. Donato & al., in preparation). Conversely, for each $`V_a/\sqrt{K_0}`$ and $`V_c`$ we varied $`L`$ and $`K_0/L`$ and the best values for the $`\chi ^2`$ are depicted in Fig. 2. We see that they are gathered in a narrow range, namely: $`8.5V_c12`$ km sec-1 and $`410V_a/\sqrt{K_0}530\text{km}\text{s}^1/\text{Mpc}\text{Myr}^{1/2}`$. We do not find any model having a good $`\chi ^2`$ without convection ($`V_c=0`$) or without reacceleration ($`V_a=0`$). We emphasize the fact that these four parameters are strongly correlated so that, when one of them is determined (for example $`L`$ with radioactive nuclei), the allowed ranges for the others will be narrower than what could naively appear from the above two figures.

The observed sub–Fe/Fe flux ratio does not give very powerful constraints. Nevertheless, all the models which give a good fit to B/C data are consistent with sub–Fe/Fe data. Figures 3 and 4 illustrate these results with the particular choice $`L=9.5\text{kpc}`$, $`K_0/L=0.00345\text{kpc}\text{Myr}^1`$, $`V_c=10.5\text{km}\text{s}^1`$ and $`V_a/\sqrt{K_0}=470\text{km}\text{s}^1/\text{Mpc}\text{Myr}^{1/2}`$, giving a reduced $`\chi _r^21.3`$. Note that all the configurations lying in the allowed ranges of Fig. 1 and 2 will provide very similar curves.

Up to now we have dealt with flux ratios, which partially hide the spectral features of the fluxes themselves. So we have also checked that the above configurations were consistent with the oxygen flux as measured by heao–3. In Fig. 5 the oxygen flux has been calculated for the same configuration as in Fig. 3 and 4. We can notice that the agreement with heao–3 data is not so satisfactory. The parameter which mainly prevents us from obtaining a good fit is the source spectrum. We recall that we set the sum $`\delta +\alpha _j`$ to the power–low index as determined in Wiebel–Sooth et al. (1998). Indeed, this number for the oxygen species is equal to $`2.68`$ as derived from heao–3 and many other measurements. If instead we fix $`\delta +\alpha _{Oxygen}`$ to be $`2.8`$ (as in Swordy & al. (1993)), the fit improves (see Fig. 5).

### 5.2 Variation of the diffusion coefficient spectral index

We also tested different values for $`\delta `$, and we find that correspondingly the diffusion parameters giving a good fit to B/C change. As an example, for a fixed value of the halo thickness $`L=3\text{kpc}`$, we find that $`\delta `$ is allowed to vary between $`0.5`$ and $`0.84`$. This is displayed in Fig. 6. In the whole parameter space, the range of $`\delta `$ extends from approximately $`0.45`$ to $`0.85`$, as one can see in Fig. 7. In particular the value $`\delta =0.33`$ corresponding to a Kolmogorov–like turbulence spectrum is strongly disfavoured ($`\chi ^2>100`$). For intermediate values of $`\delta `$, good models are obtained for the full range in $`L`$ (as in Fig. 1). For low values of $`\delta `$, models with a small halo size $`L`$ are excluded; in particular for $`\delta <0.45`$, there is no good model with $`L<15\text{kpc}`$. Finally, for high values of $`\delta `$, models with a large halo $`L`$ are excluded, and for $`\delta >0.85`$, there is no good model with $`L>1\text{kpc}`$. As can be seen in Fig. 7, each value of $`\delta `$ gives a different contour plot in the $`K_0/L`$$`L`$ plane. It appears that they all can be superimposed to a single curve by a rescaling $`K_0/LK_0/L\times f(\delta )`$, where $`f`$ is a function of $`\delta `$ only. For the contours displayed in Fig. 7, it takes the values $`f(0.46)=0.51`$, $`f(0.5)=0.62`$, $`f(0.6)1`$, $`f(0.7)=1.54`$ and $`f(0.85)=2.78`$.

In the $`V_a/\sqrt{K_0}`$$`V_c`$ plane, the values of $`V_c`$ are shifted downward as $`\delta `$ is decreased but the allowed range of $`V_a/\sqrt{K_0}`$ does not significantly move. Nevertheless, the allowed values for $`V_c`$ never reach $`0`$, so that no–wind models can be excluded.

## 6 Discussion and conclusion

We obtain good quantitative constraints on the diffusion parameters from B/C data. In particular, there is a very strong correlation between $`L`$, $`K_0/L`$, $`V_c`$, $`V_a/\sqrt{K_0}`$ and $`\delta `$. For $`\delta =0.6`$, we find that $`8.5\text{km}\text{s}^1<V_c<12\text{km}\text{s}^1`$ and $`410<V_a/\sqrt{K_0}<530`$ (where $`V_a`$ is expressed in $`\text{km}\text{s}^1`$ and $`K`$ in $`\text{kpc}^2\text{Myr}^1`$). Furthermore, we show that the power law index for the diffusion coefficient is restricted to the interval $`[0.45,0.85]`$, the best $`\chi ^2`$ being $`25.5`$ for $`\delta =0.70`$. For any $`\delta `$ in this interval, the good parameters in the $`K_0/L`$$`L`$ and $`V_a/\sqrt{K_0}`$$`V_c`$ planes can be straightforwardly deduced from the corresponding values for $`\delta =0.6`$ by a simple scaling law. We exclude any model without a convective velocity or without reacceleration for any combination of the three other diffusion parameters.

Our conclusions could get more stringent by new measurements in the whole energy spectrum for all nuclei. We emphasize that all our results were obtained using the best data, which are rather scarce and more than 20 year–old; new data are thus strongly needed. The ams experiment on board the International Space Station will have in principle the ability to provide some of these data.

The next steps of this analysis will be to study radioactive species. In particular, we expect the recent smili data (giving the spectral distribution of Be isotopes over a large energy range) to provide a new insight on cosmic ray propagation, and thus to constraint further the diffusion parameters. We are also investigating the standard antiproton signal, using the results described in this paper.

## Acknowledgments

We thank Aimé Soutoul for having pointed out a mistake in an early version of this paper. D.M. is particularly grateful for his remarks. F.D. gratefully acknowledges a fellowship by the Istituto Nazionale di Fisica Nucleare. We also would like to thank the French Programme National de Cosmologie for its financial support.

## Appendix A Solution for the diffusion model

Sources and interactions with matter are confined to the thin disk and diffusion which occurs throughout disc and halo with the same strength is independent of space coordinates. Sources follow a universal form $`Q(E)q(r)`$ with an associated normalization $`q_0^j`$ for each nucleus (see $`\mathrm{\S }`$ 3.5.2). Furthermore, we consider here only constant wind $`V_c`$ in the $`z`$ direction so that adiabatic losses (third term in (1)) vanish.

### A.1 High energy limit

The dependence in energy $`E`$ is implicit, and if we add radioactive contributions localized in both disc and halo

$`\left(_{diff}\mathrm{\Gamma }_{rad}^j\right)N^j(r,z)+{\displaystyle \underset{k=1}{\overset{j1}{}}}\mathrm{\Gamma }_{rad}^{kj}N^k+2h\delta (z)\left(q_0^jQ(E)q(r)+{\displaystyle \underset{k=1}{\overset{j1}{}}}\stackrel{~}{\mathrm{\Gamma }}^{kj}N^k(r,0)\stackrel{~}{\mathrm{\Gamma }}^jN^j(r,0)\right)=0`$ (A1)
$`_{diff}=V_c{\displaystyle \frac{}{z}}+K\left({\displaystyle \frac{^2}{z^2}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{}{r}}(r{\displaystyle \frac{}{r}})\right)`$ (A2)
In our geometry a solution can be found in term of Bessel functions of zeroth order taking advantage of their properties (Jackson 1975). One can expand all the quantities over the orthogonal set of Bessel functions $`\{J_0(\zeta _ix)\}^{i=1\mathrm{}\mathrm{}}`$,

$`N^j(r,z)={\displaystyle \underset{i=1}{\overset{\mathrm{}}{}}}N_i^j(z)J_0(\zeta _i{\displaystyle \frac{r}{R}})`$ (A3)
$`q(r)={\displaystyle \underset{i=1}{\overset{\mathrm{}}{}}}\widehat{q_i}J_0(\zeta _i{\displaystyle \frac{r}{R}}),\widehat{q_i}={\displaystyle \frac{1}{\pi R^2J_1^2(\zeta _i)}}{\displaystyle \frac{_0^1\rho q(\rho )J_0(\zeta _i\rho )𝑑\rho }{_0^1\rho q(\rho )𝑑\rho }},\rho r/R`$ (A4)
Such an expression respects automatically one of the boundary condition $`N(r=R,z)=0`$. We will use a compact notation to describe the most general form for a source term

$$\overline{𝒬^j}=q_0^jQ(E)\widehat{q_i}+\underset{k}{\overset{m_k>m_j}{}}\stackrel{~}{\mathrm{\Gamma }}^{kj}N_i^k(0)$$
(A5)
##### Stable progenitors $`𝑵^𝒌`$

As equations have to be valid order by order (second term of (A1) vanishes)

$$\left[\frac{d^2}{dz^2}\frac{V_c}{K}\frac{d}{dz}(\frac{\zeta _i^2}{R^2}+\frac{\mathrm{\Gamma }_{rad}^{N^j}}{K})\right]N_i^j(z)=\left(\frac{\overline{𝒬^j}}{K}+\frac{2h\stackrel{~}{\mathrm{\Gamma }}^j}{K}N_i^j(0)\right)\delta (z)$$
(A6)
A fast procedure to solve (A6) is given by the following three steps:

1. Solution in the halo (no right–hand side in eq. (A6)) with the boundary condition $`N_i^j(z=\pm L)=0`$
2. Integration $`(lim_{h0}_h^{+h}\mathrm{}𝑑z)`$ of equation (A6) through the thin disc555 In terms of distribution (quoted in braces), defining $`\sigma _0`$ and $`\sigma _1`$ as the discontinuities of $`0^{th}`$ et $`1^{st}`$ order, remember that
$$\frac{^2}{z^2}\{(z)\}=\left\{\frac{^2(z)}{z^2}\right\}+\sigma _1\delta (z)+\sigma _0\frac{\delta (z)}{z}$$
Imposing the continuity of the vertical cosmic ray current across the plane $`z=0`$, we thus have
$`\sigma _1lim_{ϵ0}\left[dN_i^j(z)/dz\right]_ϵ^{+ϵ}=2N_i^j(0)\frac{V_c}{K}`$ and $`\sigma _0=0`$. , which gives

$$2N_{}^{}{}_{i}{}^{j}(z)|_{z=0}2N_i^j(0)\frac{V_c}{K}2hN_i^j(0)\frac{\stackrel{~}{\mathrm{\Gamma }}^j}{K}+\overline{𝒬^j}=0$$
(A7)
3. Put the halo solution in equation (A7) to ensure continuity beetwen the two zones.

We finally obtain the solutions for stable progenitors in relativistic regime:

$`N^j(r,z)=\mathrm{exp}\left({\displaystyle \frac{V_cz}{2K}}\right){\displaystyle \underset{i=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\overline{𝒬^j}}{A_i^j}}{\displaystyle \frac{\mathrm{sinh}\left[\frac{S_i^j(Lz)}{2}\right]}{\mathrm{sinh}\left[\frac{S_i^jL}{2}\right]}}J_0(\zeta _i{\displaystyle \frac{r}{R}})`$ (A8)
$`\overline{𝒬^j}q_0^jQ(E)\widehat{q_i}+_k^{m_k>m_j}\stackrel{~}{\mathrm{\Gamma }}^{kj}N_i^k(0)`$ (A9)
$`S_i^j(\frac{V_c^2}{K^2}+4\frac{\zeta _i^2}{R^2}+4\frac{\mathrm{\Gamma }_{rad}^{N^j}}{K})^{1/2}A_i^j2h\stackrel{~}{\mathrm{\Gamma }}_{N^j}^{tot}+V_c+KS_i^j\mathrm{coth}(\frac{S_i^jL}{2})`$ (A10)
For a primary $`\overline{𝒬^j}=q_0^jQ(E)\widehat{q_i}`$, and for a pure secondary $`\overline{𝒬^j}=_k^{m_k>m_j}\stackrel{~}{\mathrm{\Gamma }}^{kj}N_i^k(0)`$. Note that solutions given in Webber et al. (1992) for secondary takes advantage of the primary form of $`N_i^k(0)`$. Since we are here interested in a shower–like (see $`\mathrm{\S }`$ 3.2.2) resolution, the form given here is more adapted.

##### $`𝜷`$ decay contribution from $`𝑵^𝒌`$

For all the nuclei treated here, $`N^j`$ never has more than one unstable contribution, so that the sum over $`k`$ for $`N_{rad}^k`$ reduces to one term in equation (A1). Resolution is complicated by the localisation of this source in the whole halo. Focalising on this specific term, neglecting for a while primary source and classical spallative secondary contribution $`2h\delta (z)_k\stackrel{~}{\mathrm{\Gamma }}^{kj}N^k(r,0)`$, one obtains (following the same procedure as described in the previous section666The contribution of these radioactive nuclei may be unimportant in some cases, but we should take it into account as it is the dominant process for some others. In the simple example of 10Be$`^{10}`$B, neglecting this channel would give an error of about $`10\%`$ on the B flux, whereas considering that this term is only located in the disc would give an error of about $`3\%`$ compared to the rigourous treatment given above. Notice finally that at fixed energy per nucleon, the rigidity depends on the nuclear species at stake. The diffusion coefficient $`K^j`$ of the child nucleus $`j`$ is therefore different from its progenitor’s one $`K^k`$. The difference $`K^jK^k`$ tends to vanish for the heaviest nuclei.)

$`N_{\mathrm{\Gamma }_{rad}^k}^j(r,z)=`$ $`{\displaystyle \underset{i=0}{\overset{\mathrm{}}{}}}J_0(\zeta _i{\displaystyle \frac{r}{R}})\times {\displaystyle \frac{\mathrm{\Gamma }_{rad}^{kj}}{K^j(a_i^2a^2)}}{\displaystyle \frac{N_i^k(0)}{\mathrm{sinh}\left(\frac{S_i^kL}{2}\right)}}`$ (A11)
$`\times `$ $`\left\{\begin{array}{cc}\alpha a_i\mathrm{cosh}[\frac{S_i^j(Lz)}{2}]e^{(\frac{V_cz}{2K^j})}+\left(a\mathrm{sinh}[\frac{S_i^k(Lz)}{2}]+a_i\mathrm{cosh}[\frac{S_i^k(Lz)}{2}]\right)e^{(\frac{V_cz}{2K^k})}\hfill & \\ +\frac{e^{(\frac{V_cz}{2K^j})}}{A_i^j}\frac{\mathrm{sinh}\left[\frac{S_i^j(Lz)}{2}\right]}{\mathrm{sinh}\left[\frac{S_i^jL}{2}\right]}\left[\begin{array}{ccc}\alpha a_i\mathrm{cosh}\left[\frac{S_i^jL}{2}\right]\left(V_c+2h\stackrel{~}{\mathrm{\Gamma }}^j+K^jS_i^j\mathrm{tanh}\left[\frac{S_i^jL}{2}\right]\right)\hfill & & \\ +\mathrm{sinh}\left[\frac{S_i^kL}{2}\right]\left[a\left(\frac{V_cK^j}{K^k}2V_c2h\stackrel{~}{\mathrm{\Gamma }}^j\right)a_iK^jS_i^k\right]\hfill & & \\ +\mathrm{cosh}\left[\frac{S_i^kL}{2}\right]\left[a_i\left(\frac{V_cK^j}{K^k}2V_c2h\stackrel{~}{\mathrm{\Gamma }}^j\right)aK^jS_i^k\right]\hfill & & \end{array}\right]\hfill & \end{array}\right\}`$ (A17)
where $`S_i^j`$ and $`A_i^j`$ have already been defined in (A10) and

$`\alpha \mathrm{exp}\left(\frac{V_cL}{2}\left(\frac{1}{K^k}\frac{1}{K^j}\right)\right),a\frac{V_c^2}{2K^k}\left(\frac{1}{K^k}\frac{1}{K^j}\right)+\frac{\mathrm{\Gamma }_{rad}^{kj}}{K^k}\frac{\mathrm{\Gamma }_{rad}^j}{K^j},a_i\frac{S_i^kV_c}{2}\left(\frac{1}{K^k}\frac{1}{K^j}\right)`$ (A18)
### A.2 Full solution

Modifying relation (A1) so as to take into account energy losses and diffusive reacceleration is straightforward since those processes take place only in the disc and not in the halo. Following step by step the three point procedure previously described, one gets the differential equation777For the $`\beta `$–unstable progenitor, the derivation is a bit more tedious but leads to the same final compact result.

$$A_i^jN_i^j(0)=\overline{𝒬^j}2h\frac{}{E}\left\{b_{loss}^j(E)N_i^j(0)K_{EE}^j(E)\frac{}{E}N_i^j(0)\right\}$$
(A19)
where $`b_{loss}^j(E)`$ and $`K_{EE}^j`$ respectively correspond to energy losses (described in section 3.6.1 and defined by (21)) and energy diffusion (reacceleration term given by (20)).


format: text

text:
# 1 Introduction

## 1 Introduction

The CKM matrix element $`V_{ub}`$ plays an important role in the determination of the unitarity triangle. The cleanest method to obtain the absolute value $`|V_{ub}|`$ is through the measurement of semileptonic $`bu`$ transitions, which will eventually give the length of one of the sides of the unitarity triangle.

The possibilities of determining $`V_{ub}`$ from semileptonic decays have been studied in detail in the BaBar Workshop . Here the conclusion was reached that the determination of $`V_{ub}`$ will be performed by a mixture of different methods.

One of the theoretically cleanest possibilities is to use inclusive semileptonic decays. Placing a cut on the hadronic invariant mass of the final state can in principle eliminate the charm contribution, which otherwise would be overwhelming. This cut on the hadronic invariant mass can be implemented at the asymmetric $`B`$ factories, making this method experimentally feasible.

From the theoretical side the decay rate, including cuts on the lepton energy as well as on the hadronic invariant mass of the final state, can be computed systematically within the framework of the $`1/m_b`$ expansion, except for certain regions of phase space, where the $`1/m_b`$ expansion has to be replaced by an expansion in twist. To describe these regions of phases space one has to introduce a so-called “shape function”, which in principle introduces a large hadronic uncertainty.

Another method that has been proposed  avoids the twist expansion and relies only on a standard $`1/m_Q`$ expansion. This method needs a measurement of the lepton invariant mass spectrum, in which the regions of phase space where the shape function plays a role are kinematically suppressed. This method will also allow a clean determination of $`V_{ub}`$.

An alternative approach has been put forth in  where the factorization into soft, jet and hard sub-processes  has been employed to relate the radiative $`b`$ decays to the semileptonic ones in a way which explicitly reduces the impact of the shape function uncertainties. Then a prediction of the ratio $`|V_{ub}/V_{ts}|`$ can be obtained with a good accuracy in a model independent way. This approach is similar to the one in .

Exclusive decays will open a completely different window on $`V_{ub}`$; however, in these decays a certain model dependence seems to be unavoidable, unless lattice data become reasonably precise.

From this variety of methods to determine $`|V_{ub}|/|V_{cb}|`$, we shall expand in this paper on the one in which a cut is applied on the hadronic invariant mass in semileptonic $`B`$ decays to filter out the $`bu`$ transitions. The advantage is that the hadronic invariant mass may be easier to measure, however, this method involves the shape function and potentially has larger theoretical uncertainties than the inclusive method using the leptonic invariant mass.The method based on the hadronic invariant mass spectrum has already been discussed in , where the main focus was on the perturbative contributions of order $`\alpha _s`$ and $`\beta _0\alpha _s^2`$.

The purpose of this paper is to consider this approach in detail and to try to estimate the uncertainties, including perturbative as well as non-perturbative contributions and, in addition, a cut on the lepton energy. It turns out that the main uncertainties originate from the heavy quark mass $`m_b`$ (or, equivalently, from $`\overline{\mathrm{\Lambda }}=M_Bm_b`$ where $`M_B`$ is the $`B`$-meson mass) and the strong coupling $`\alpha _s`$, while the uncertainties introduced by the shape function are small.

The next section deals with the kinematics. Then, in Section 3, we give the radiative corrections for the partonic process $`bu\mathrm{}\overline{\nu }_{\mathrm{}}`$ to order $`\alpha _s`$. In section 4 we include the leading-twist non-perturbative effects; this requires the introduction of the light-cone distribution function for the heavy quark, for which we use a simple parametrization. We combine the perturbative and non-perturbative corrections and study the uncertainties in the determination of $`|V_{ub}|/|V_{cb}|`$ in section 5.

## 2 Kinematical Relations and Definitions

We shall first define the kinematic variables for the partonic process $`bu\mathrm{}\overline{\nu _{\mathrm{}}}`$. Although we also use results for the semileptonic $`bc`$ decay rate, the latter is considered in the partonic framework solely. In fact, we only need the total rate with a single lower cut on the electron energy for $`bc`$; we thus refer the reader to  for further discussion of the kinematics of this process. The initial state $`b`$ quark has a momentum $`p_b=m_bv`$, where $`v`$ is the velocity of the $`B`$ meson. With the momentum transfer to the leptons $`q=k+k^{}`$ the variable $`p=qm_bv`$ is the partonic momentum of the final state. Writing $`E_{\mathrm{}}=vk`$ for the energy of the lepton, $`E_p=vp`$ for the partonic energy of the final state and $`p^2`$ for the partonic invariant mass, we define the following re-scaled variables

$$x=\frac{2E_{\mathrm{}}}{m_b},x_p=\frac{2E_p}{m_b},y=\frac{q^2}{m_b^2},z=\frac{p^2}{m_b^2}.$$
(1)
One of these variables is superfluous when calculating the triple differential decay rate and can be substituted by using the relation

$$x_p+yz=1.$$
(2)
Using these definitions, we compute the triple differential partonic rate to order $`\alpha _s`$ and order $`1/m_b`$, and split it into the tree-level term $`\mathrm{\Gamma }^0`$ and the $`𝒪(\alpha _s)`$ correction $`\mathrm{\Gamma }^1`$:

$$\text{d}\mathrm{\Gamma }^{pert}=\text{d}\mathrm{\Gamma }^{parton,0}+\frac{2\alpha _s}{3\pi }\text{d}\mathrm{\Gamma }^{pert,1}.$$
(3)
The radiative corrections to $`𝒪(\alpha _s)`$ have recently been calculated by ; we have checked their result and find full agreement with ours. The formulae up to terms $`𝒪(\alpha _s)`$ are quite tedious and can be found in . However, we present them in the Appendix in a form suitable for our numerical evaluation.

We shall now discuss the hadronic kinematics. The momentum of the initial $`B`$ meson is $`M_Bv=(m_b+\overline{\mathrm{\Lambda }})v`$, where we have used the relation between the $`B`$ meson and the $`b`$-quark mass

$$M_B=m_b+\overline{\mathrm{\Lambda }}$$
(4)
which holds to leading order in the $`1/m_b`$ expansion. Consequently, the hadronic mass of the final state is $`M_X^2=(M_Bvq)^2=(p+\overline{\mathrm{\Lambda }}v)^2=p^2+2E_p\overline{\mathrm{\Lambda }}+\overline{\mathrm{\Lambda }}^2`$ and hence involves both the partonic invariant mass and the partonic energy. Thus we have

$`{\displaystyle \frac{\text{d}\mathrm{\Gamma }^{parton}}{\text{d}M_X^2\text{d}E_l}}(m_b)`$ $`=`$ $`{\displaystyle \underset{0}{\overset{1}{}}}\text{d}x{\displaystyle \underset{1x}{\overset{2x}{}}}\text{d}x_p{\displaystyle \underset{z_{min}}{\overset{z_{max}}{}}}\text{d}z{\displaystyle \frac{\text{d}^3\mathrm{\Gamma }^{pert}}{\text{d}x\text{d}x_p\text{d}z}}\times `$ (5)
$`\delta (M_X^2m_b^2zx_p\overline{\mathrm{\Lambda }}m_b\overline{\mathrm{\Lambda }}^2)\delta \left(E_l{\displaystyle \frac{m_bx}{2}}\right).`$
where the kinematic limits of the $`z`$ integration depend on the cuts as well as on $`x_p`$. They are

$`z_{min}=\{\begin{array}{cc}0\hfill & \text{for }x_p1,\hfill \\ x_p1\hfill & \text{for }x_p>1,\hfill \end{array}`$ (8)
and

$$z_{max}=(1x)(x_p+x1).$$
(9)
## 3 Perturbative corrections

The analysis we perform requires the knowledge of the triple differential partonic rate to order $`\alpha _s`$, as can already be seen from Eq. (5) as well as from the final formula, Eq. (14). The Born approximation is proportional to the delta function of argument $`z`$,

$$\frac{\text{d}^3\mathrm{\Gamma }^{pert,0}}{\text{d}x\text{d}x_p\text{d}z}=12\mathrm{\Gamma }_0(2xx_p)(x+x_p1)\delta (z),$$
(10)
where

$$\mathrm{\Gamma }_0=\frac{G_F^2m_b^5}{192\pi ^3}.$$
(11)
As is well known, the virtual one-loop correction contains an infrared divergence that cancels with the real gluon emission. We shall discuss in the next section how to include the leading-twist non-perturbative effects by “smearing” over some small window in the hadronic invariant mass. While choosing an excessively small value for this window would still yield a meaningless result, a region of size $`\overline{\mathrm{\Lambda }}m_b`$ is believed to exist, where this smearing provides a realistic approximation. Instead of keeping track of the virtual and real parts, one can now use an appropriately integrated distribution, which is subjected to the smearing procedure. The tree-level term is rather simple to implement numerically but, on the other hand, we have found it very useful to perform one integration of the one-loop correction analytically. On inspection of Eq. (14) in conjunction with Eq. (5), it is clear that the following integral is of use:

$$F(x,x_p,z)=_{z_{min}}^z\text{d}z^{}\frac{\text{d}\mathrm{\Gamma }^{pert,1}}{\mathrm{\Gamma }_0\text{d}x\text{d}x_p\text{d}z^{}}.$$
(12)
This function is a lengthy and tedious expression, and it is given in the Appendix.

## 4 Leading-twist non-perturbative corrections

It has been shown that the leading-twist non-perturbative corrections can be implemented at tree level by redefining the heavy-quark mass and a subsequent convolution with a so-called shape function . This convolution corresponds to an integration over the light-cone variable $`k_+`$, namely

$$\text{d}\mathrm{\Gamma }^{hadron}=\underset{m_b}{\overset{\overline{\mathrm{\Lambda }}}{}}\text{d}\mathrm{\Gamma }^{parton}(m_b+k_+)f(k_+)\text{d}k_+$$
(13)
in $`\text{d}\mathrm{\Gamma }^{HQET}`$. Although this formula is quite suggestive , it has been shown recently that it contains,in fact, spurious contributions of sub-leading twist . Furthermore, once radiative corrections are implemented, this simple convolution formula will probably no longer hold . However, there is no fully consistent way yet to include radiative corrections into this convolution, and hence we proceed in a naive way as suggested in .

This naive way is to actually use the convolution formula (13) also beyond tree level, which is at least as consistent as using the ACCMM model  beyond tree level, which is common practice. In fact, the connection between the ACCMM model and the shape function formalism has been pointed out in .

It is convenient to change the variable of integration according to

$$\text{d}\mathrm{\Gamma }^{hadron}=\underset{0}{\overset{M_B}{}}\text{d}\mathrm{\Gamma }^{parton}(m^{})f(m^{})\text{d}m^{}$$
(14)
with $`m^{}=m_b+k_+`$.

The shape function is a nonperturbative function which has to be determined either from experiment or by some model. A few relations are known for the moments of the shape function $`k_+^n`$:

$$k_+^0=1,k_+^1=0\text{and}k_+^2=\frac{1}{3}\lambda _1$$
(15)
For our study we shall use an ansatz for the shape function. Taking a simple three parameter function

$$f(x)=N\text{e}^{cx}(1x)^a,x=1\frac{M_Bm^{}}{\overline{\mathrm{\Lambda }}},$$
(16)
we can express the parameters of this ansatz in terms of the HQET parameters (15)

$$N=\frac{c^c}{\overline{\mathrm{\Lambda }}\text{e}^c\mathrm{\Gamma }(c)},a=c1,c=\frac{3\overline{\mathrm{\Lambda }}^2}{\lambda _1}.$$
(17)
We shall discuss the dependence of our results on this ansatz in the next section.

## 5 The measurement of $`|V_{ub}|`$

For the measurement of $`|V_{ub}|`$ we propose to study semileptonic $`B`$ decays with certain cuts. The first cut is on the lepton energy $`E_l`$, which is mainly given by the experimental limitations on detecting electrons with small momenta. The second cut is on the hadronic invariant mass $`M_X`$ of the final state, which serves to suppress charm. We define the semileptonic rates including cuts as

$$\mathrm{\Gamma }(M_{cut}^2,E_{cut})=\underset{0}{\overset{M_{cut}^2}{}}\text{d}M_X^2\underset{E_{cut}}{\overset{\frac{M_B^2M_X^2}{2M_B}}{}}\text{d}E_l\frac{\text{d}\mathrm{\Gamma }(BX_ue\overline{\nu _e})}{\text{d}M_X^2\text{d}E_l},$$
(18)
Clearly the region $`0<M_X<M_D`$ is dominated by $`bu`$ transitions, and we thus use in this region the expressions for $`bu`$ decays only.

We shall normalize everything to the rate with no cut on the hadronic invariant mass, and thus obtain, for the ratio $`r`$:

$$\frac{|V_{ub}|^2}{|V_{cb}|^2}r=\frac{\mathrm{\Gamma }(M_{cut}^2,E_{cut})}{\mathrm{\Gamma }(M_B,E_{cut})},$$
(19)
where in the denominator we can safely take into account the $`bc`$ channel only, the $`bu`$ contribution being only about 1%.

For the charmless decay rate entering the numerator, we include the results discussed above. The perturbative corrections are taken into account to one loop , while the non-perturbative ones are included to leading twist. The $`bc`$ decay rate in the denominator is evaluated including $`𝒪(\alpha _s)`$ corrections.

Unlike the approach advertised in  and , the shape of the hadronic invariant mass spectrum depends on the shape function, for which we use the parametrization (16). However, as it will turn out, the precise form of the shape function is irrelevant for the ratio (19), the main sources of uncertainties are the value of the strong coupling and that of the quark masses. The latter may be replaced by $`\overline{\mathrm{\Lambda }}`$, we thus use $`m_b=\overline{M}_B\overline{\mathrm{\Lambda }}`$ and $`m_c=\overline{M}_D\overline{\mathrm{\Lambda }}`$, where $`\overline{M}`$ denotes the spin-averaged meson mass.

The dependence on the shape function can be studied by checking the sensitivity of our result with respect to higher moments of this function. Since we have not included any $`1/m_b^2`$ terms in our calculation the ratio (19) should be reasonably insensitive already to the second moment of the shape function, which is given in terms of the kinetic energy operator $`\lambda _1`$ in (15). This forces us to make $`M_{cut}`$ as large as possible.

We have examined the ratio defined in Eq. (19), allowing the variations $`0.4\text{ GeV}\overline{\mathrm{\Lambda }}0.75\text{ GeV}`$, $`0.2\alpha _s0.3`$ and $`0.6\text{ GeV}^2\lambda _10.1\text{ GeV}^2`$.

The main uncertainty is induced by $`\alpha _s`$ and $`\overline{\mathrm{\Lambda }}`$, where $`\overline{\mathrm{\Lambda }}`$ is equivalent to the heavy quark mass. It has been argued that these two quantities are correlated; the size of the radiative corrections depends on the particular choice of the mass.

Using the pole mass scheme, it has been shown that the radiative corrections are large and the perturbation series converges very slowly. Treating both the pole mass and the perturbative contributions independently, the uncertainties of the quantities would simply add, leaving us with a large (and certainly overestimated) uncertainty. In this case we get (see the lighter band in Fig. 1):

$$1.15r1.96\text{at }M_{cut}^2=4\text{ GeV}^2$$
(20)
where we have used $`m_b^{pole}=(4.75\pm 0.15)`$ GeV and $`\alpha _s`$ between 0.2 and 0.3. In this way we obtain a theoretical uncertainty of $`25\%`$ in the determination of the ratio.

Another option is to switch to a short-distance mass definition, such as $`m_b^{\overline{MS}}`$ by replacing

$$m_b^{pole}=m_b^{\overline{MS}}\left(1+\frac{4}{3\pi }\alpha _s\right),$$
(21)
which reduces the size of the coefficients of the perturbation series, and the pertrubative uncertainties thus become smaller. Using a recent value for $`M_b^{\overline{MS}}`$

$$m_b^{\overline{MS}}=(4.25\pm 0.08)\text{GeV},$$
(22)
we arrive at an estimate for $`r`$ with a smaller uncertainty (the darker band in Fig. 1)

$$1.28r1.85\text{at }M_{cut}^2=4\text{ GeV}^2.$$
(23)
In order to display the dependence on the input parameters, we choose the “average” values of the three parameters:

$$\overline{\mathrm{\Lambda }}^{aver}=0.55\text{ GeV},\alpha _s^{aver}=0.25,\lambda _1^{aver}=0.35\text{ GeV}^2,$$
(24)
and obtain up to linear terms in the variations

$$r=1.580.86\frac{\mathrm{\Delta }\overline{\mathrm{\Lambda }}}{\overline{\mathrm{\Lambda }}^{aver}}0.26\frac{\mathrm{\Delta }\alpha _s}{\alpha _s^{aver}}+0.06\frac{\mathrm{\Delta }\lambda _1}{\lambda _1^{aver}}.$$
(25)
This explicitly shows that the dependence on $`\lambda _1`$, which is the second moment of the shape function, is weak. The dependence on even higher moments is expected to be further suppressed by inverse powers of the heavy quark mass.

In conclusion, our best estimate for $`r`$ is

$$r=1.57\pm 0.3,$$
(26)
which corresponds to a theoretical uncertainty in the determination of $`|V_{ub}/V_{cb}|`$ of about ten percent

$$\left(\frac{\mathrm{\Delta }|V_{ub}/V_{cb}|}{|V_{ub}/V_{cb}|}\right)_{theor}10\%$$
(27)
## 6 Conclusions

We have performed a detailed analysis of one of the possibilities to obtain $`V_{ub}`$ from inclusive semileptonic $`B`$ decays by placing a cut on the hadronic invariant mass to get rid of the charm background. This method has been criticized since it depends on the shape function, which describes the endpoint of the hadronic invariant mass spectrum. This function is not very well known and hence it has to be modelled, which will introduce some systematic uncertainty. However, integrating over the window in hadronic invariant masses relevant to $`bu`$ transitions, we have shown that the dependence on the shape function is much smaller than the uncertainties induced by the quark mass and by the truncation of the perturbative series.

The the ratio between the semileptonic rates including a cut on the hadronic invariant mass and the semileptonic rate without a cut yields $`|V_{ub}/V_{cb}|`$ up to a quantity $`r`$, which we have computed in leading twist approximation and to order $`\alpha _s`$. Based on our calculations the uncertainty in this quantity is $`20\%`$ leaving us with a $`10\%`$ theoretical uncertainty in the determination of $`V_{ub}`$. This method is thus one of the cleanest possible to obtain $`V_{ub}`$ at the ongoing $`B`$-factory experiments.

## Acknowledgements

TM and BP acknowledge support from the DFG Graduiertenkolleg “Elementarteilchenphysik an Beschleunigern”and from the DFG Forschergruppe “Quantenfeldtheorie, Computeralgebra und Monte Carlo Simulationen”. TM also acknowledges support from the BMBF. This work is partly supported by the KBN grants 2P03B05418 and 5P03B09320 and by the European Commission 5th Framework contract HPRN-CT-2000-00149. PU would like to thank the Polish–French Collaboration within IN2P3 through Annecy.

## Appendix

We present here the convoluted spectra in some more detail, showing in particular a convenient way of calculating the contribution from the terms proportional to $`\alpha _s`$. The convolution of the partonic spectra results in the leading-twist approximated rate of $`BX_ul\overline{\nu }`$ decay. The resulting distribution can be written in the following form:

$`{\displaystyle \frac{\text{d}\mathrm{\Gamma }}{\text{d}M_X^2\text{d}E_l}}`$ $`=`$ $`{\displaystyle _0^1}\text{d}x{\displaystyle _{1x}^{2x}}\text{d}x_p{\displaystyle _{z_{min}}^{z_{max}}}\text{d}z{\displaystyle _{\mathrm{}}^{M_B}}\text{d}m^{}f(m^{}){\displaystyle \frac{\text{d}\mathrm{\Gamma }^{parton}}{\text{d}x\text{d}x_p\text{d}z}}\times `$ (A.1)
$`\delta (M_X^2m^2zx_p\overline{\mathrm{\Lambda }}^{}m^{}\overline{\mathrm{\Lambda }}^2)\delta (E_l{\displaystyle \frac{m^{}x}{2}}).`$
The partonic rate can be split into the tree level term and the $`𝒪(\alpha _s)`$ correction,

$$\text{d}\mathrm{\Gamma }^{parton}=\text{d}\mathrm{\Gamma }^{parton,0}+\frac{2\alpha _s}{3\pi }\text{d}\mathrm{\Gamma }^{parton,1},$$
(A.2)
and the convoluted rate can then be divided up accordingly. Using the delta functions, the integrations involved in these rates can be simplified, yielding

$`{\displaystyle _{E_{cut}}^{M_B/2}}\text{d}E_l{\displaystyle _0^{M_{cut}^2}}\text{d}M_X^2{\displaystyle \frac{\text{d}\mathrm{\Gamma }^{(0)}(BX_ue\overline{\nu _e})}{\text{d}M_X^2\text{d}E_l}}=`$
$`{\displaystyle _{E_{cut}}^{M_B/2}}\text{d}E_l{\displaystyle _{2E_l/M_B}^1}\text{d}x{\displaystyle _{y_{min}}^x}\text{d}y{\displaystyle \frac{2}{x}}f\left({\displaystyle \frac{2E_l}{x}}\right){\displaystyle \frac{\text{d}\mathrm{\Gamma }^{parton,0}}{\text{d}x\text{d}y}},`$ (A.3)
$$y_{min}=\mathrm{max}\{0,\frac{1M_{cut}^2(M_B2E_l/x)^2}{(2E_l/x)(M_B2E_l/x)}\}$$
(A.4)
for the tree-level contribution, while the $`\alpha _s`$ correction contributes

$`{\displaystyle _{E_{cut}}^{M_B/2}}\text{d}E_l{\displaystyle _0^{M_{cut}^2}}\text{d}M_X^2{\displaystyle \frac{\text{d}\mathrm{\Gamma }^{(1)}(BX_ue\overline{\nu _e})}{\text{d}M_X^2\text{d}E_l}}=`$
$`{\displaystyle _{E_{cut}}^{M_B/2}}\text{d}E_l{\displaystyle _{2E_l/M_B}^1}\text{d}x{\displaystyle _{1x}^{2x}}\text{d}x_p{\displaystyle _{z_{min}}^{z_{max}}}\text{d}z{\displaystyle \frac{2}{x}}f\left({\displaystyle \frac{2E_l}{x}}\right){\displaystyle \frac{\text{d}\mathrm{\Gamma }^{parton,1}}{\text{d}x\text{d}x_p\text{d}z}},`$ (A.5)
where

$$z_{min}=\{\begin{array}{cc}0,& \text{f}orx_p1,\\ x_p1& \text{f}orx_p>1,\end{array}$$
(A.6)
and

$$z_{max}=\mathrm{min}\{\frac{M_{cut}^2x_p(M_B2E_l/x)(M_B2E_l/x)^2}{2E_l/x},(1x)(x_p+x1)\}.$$
(A.7)
With these formulae, it is easy to numerically convolute the Born approximated rate. However, it is useful to eliminate one integral from the convolution of the $`\alpha _s`$ term. To this end, the integral over $`z`$ in Eq. (Appendix) has been performed analytically, so that we use the prime function $`F(x,x_p,z)`$ defined as

$$F(x,x_p,z)=_{z_{min}}^z\text{d}z^{}\frac{\text{d}\mathrm{\Gamma }^{parton,1}}{\mathrm{\Gamma }_0\text{d}x\text{d}x_p\text{d}z^{}},$$
(A.8)
where $`\mathrm{\Gamma }_0`$ is defined in $`(\text{11})`$. Then,

$$F=\{\begin{array}{cc}F_1(z)F_1(0)+F_2(z)+F_3,\hfill & x_p<1,\hfill \\ F_1(z)F_1(z_{min})+F_2(z)F_2(z_{min})\hfill & x_p>1,\hfill \end{array}$$

where

$$z_{min}=x_p1,$$
(A.9)
and

$`F_1`$ $`=`$ $`{\displaystyle \frac{1}{8}}\left[{\displaystyle \frac{1}{t}}\mathrm{ln}{\displaystyle \frac{1+t}{1t}}[x_p^24(x1+x_p/2)^2/t^2]+8(x1+x_p/2)^2/t^2\right]`$ (A.10)
$`\times (3639x_p+12x_p^23x_p^3/4)+{\displaystyle \frac{1}{t}}\mathrm{ln}{\displaystyle \frac{1+t}{1t}}(x1+x_p/2)[4878x_p`$
$`+45x_p^221x_p^3/2+3x_p^4/4+(x1+x_p/2)(4215x_p+27x_p^2`$
$`21/8x_p^3)]+(A_1t+A_2)\mathrm{ln}{\displaystyle \frac{1+t}{1t}}+A_3\mathrm{log}(1t)+A_4,`$
$`F_2`$ $`=`$ $`12[2\mathrm{ln}^2z(8\mathrm{ln}x_p7)\mathrm{ln}z]v(\overline{x}x_p),`$ (A.11)
$`F_3`$ $`=`$ $`(x_p\overline{x})\{24\mathrm{ln}x_p(1+5v4v\mathrm{ln}x_p)`$ (A.12)
$`+v[16\pi ^26048\text{Li}_2(1x_p)]\}.`$
In the above formulae,

$`A_1`$ $`=`$ $`(x_p\overline{x})(x_pz\overline{x}/2+21x_p\overline{x}/26x_pv12x_p^2\overline{x}+3x_p^2v+7x_p^3\overline{x}/4)`$ (A.13)
$`+(2x_pzv9x_pz/29/10x_pz^2x_p^2zv+7x_p^2z+12x_p^2v+15x_p^2/2`$
$`21x_p^3z/20+17x_p^3v/2105/4x_p^35/4x_p^4v+53/4x_p^4`$
$`219/160x_p^5),`$
$`A_2`$ $`=`$ $`(x_p\overline{x})(36x_p(1x_px)+6x_p^2v+4x_p^3\overline{x}+12\overline{x}96v\mathrm{ln}2+48v)`$ (A.14)
$`24x_p+48x_p^2+16x_p^3v60x_p^32x_p^4v+26x_p^412/5x_p^5,`$
$`A_3`$ $`=`$ $`2A_2+96(x_p\overline{x})(1+\overline{x}x_p)\mathrm{ln}(1+t),`$ (A.15)
$`A_4`$ $`=`$ $`(x_p\overline{x})(5x_pz\overline{x}24z\overline{x}+12zv)+32x_pzv81x_pz21/5x_pz^2`$ (A.16)
$`4x_p^2zv+40x_p^2z81/20x_p^3z+36zv+12z6z^2v+24z^2,`$
and

$$\overline{x}=1x,v=2xx_p,t=\sqrt{14z/x_p^2}.$$
(A.17)


text:
# REFERENCES

Pattern Competition in the Photorefractive Semiconductors

Yuo-Hsien Shiau1,∗ and Chin-Kun Hu2,3,+

1 National Center for Theoretical Sciences (Physics Division) 101, Section 2 Kuang Fu Road,

Hsinchu, Taiwan 300, ROC

2 Institute of Physics, Academia Sinica, Taipei 11529, Taiwan

3Department of Physics, National Dong Hwa University, Hualien 97401, Taiwan

## Abstract

We analytically study the photorefractive Gunn effect in n-GaAs subjected to two external laser beams which form a moving interference pattern (MIP) in the semiconductor. When the intensity of the spatially independent part of the MIP, denoted by $`I_0`$, is small, the system has a periodic domain train (PDT), consistent with the results of linear stability analysis. When $`I_0`$ is large, the space-charge field induced by the MIP will compete with the PDT and result in complex dynamics, including driven chaos via quasiperiodic route.

PACS number(s): 72.20.Ht, 42.70.Nq, 05.45.-a, 05.45.Yv

Keywords: A. Semiconductors; D. Pattern Competition; E. Nonlinear Analysis

Spatiotemporal behavior of Gunn domains  in semiconductors is an interesting and important current research problem . Gunn diodes are not only good systems for studying nonlinear dynamics, but also potentially important for various applications . In 1996-1997 Segev et al.  and Subačius et al.  showed, respectively, that optical waves can excite multiple Gunn-domain formation with spatial periodicity, called periodic domain trains (PDT), in deep-impurity doped GaAs and semi-insulating GaAs; both groups used the well known Kroemer’s criterion  to determine the number of Gunn domains. In 1998, Bonilla et al.  numerically studied the same dynamical system as Segev et al.; however, their findings are mostly inconsistent with the predictions of Segev et al.. Bonilla et al.  found that the number of Gunn domains in the system is not consistent with the prediction of Segev et al. and in many cases the system shows complex or chaotic behavior. Very recently, Shiau et al. used a linear stability analysis (LSA) to study a shallow-impurity doped GaAs with length $`L`$ subjected to two external laser beams which form a moving interference pattern (MIP) of intensity $`I(x,t)=I_0[1+m\text{cos}(Kx+\mathrm{\Omega }t)]`$ and predicted PDT in the Gunn diode . Here $`\mathrm{\Omega }`$ is the frequency difference of the two optical waves, $`K=2\pi /\mathrm{\Lambda }`$ is the interference wave number, $`\mathrm{\Lambda }`$ is the grating period and is much smaller than $`L`$, $`m`$ is the modulation depth, and $`I_0`$ is the average intensity. With consideration of time-delay feedback, the diode will emit high-dimensional chaotic microwaves, which is potentially useful for secure microwave communications, memory devices, applications involving photorefractive effects, etc. .

The purpose of this Communication is to clarify different findings of , , and  using a nonlinear analytic method to study the system of  as $`L\mathrm{}`$. In this approach, the intensity of the spatially independent part of the MIP, $`I_0`$, plays an important role in determining the simple dynamics (i.e., PDT) or complex dynamics of Gunn domains. When $`I_0`$ is small enough, the system shows simple dynamics consistent with the result of the LSA, and Kroemer’s criterion is still valid. When $`I_0`$ is large, the space-charge field induced by the MIP will compete with the PDT and result in complex dynamics. This picture is also useful for understanding the results of other groups  discussed above.

We consider two optical waves incident on a $`n`$-GaAs with a biased voltage $`V`$. The photon energy is just above the bandgap of GaAs, i.e., 1.42 eV, so that electron-hole pairs can be generated by optical excitation. The generation-recombination processes include complete thermal ionization of the donor, generation of electron-hole pairs by the MIP with rate $`g`$, and recombination of electron-hole pairs with rate $`\gamma `$. Therefore we can write the dynamical equations for the electric field $`E`$, the free electron density $`n`$, and the free hole density $`p`$ as:

$$\frac{E}{x}=\frac{e}{ϵ}(nN_D^{}p),$$
(1)
$$\frac{n}{t}=gI(x,t)\gamma np\frac{}{x}\left[nv(E)D_n\frac{n}{x}\right],$$
(2)
$$\frac{p}{t}=gI(x,t)\gamma np+\frac{}{x}\left[p\mu _pE+D_p\frac{p}{x}\right],$$
(3)
$$V=_0^LE(x,t)𝑑x.$$
(4)
Eq. (1) is for the Gauss law. Eq. (2) and Eq. (3) represent the continuity equations for the electrons and holes, respectively. The circuit condition is given at Eq. (4). Here $`v(E)`$ displays an N-shaped negative differential mobility (NDM). The homogeneous solution of Eqs. (1)-(4) is $`E_0=V/L`$, $`n_0=(N_D^{}+\sqrt{N_D^2+4gI_0/\gamma })/2`$, and $`p_0=(N_D^{}+\sqrt{N_D^2+4gI_0/\gamma })/2`$. When $`I_0`$ is very small, $`n_0N_D^{}+gI_0/N_D^{}\gamma `$, $`p_0gI_0/N_D^{}\gamma `$, and a natural expansion parameter for high order corrections is $`\overline{ϵ}=p_0/N_D^{}=gI_0/N_{D}^{}{}_{}{}^{2}\gamma `$. We can rearrange the model equations by the following rescaling process for the dynamical variables and parameters in Eqs. (1)-(4): $`EE,nn,p\overline{ϵ}p,I\overline{ϵ}I.`$ The small parameter $`\overline{ϵ}`$ can distinguish the large variables (i.e., $`E`$ and $`n`$) from small variable $`p`$ in our model. Rescaling variables in Eqs. (1)-(4), we have  $`E/x=e(nN_D^{})/ϵ+O(\overline{ϵ})`$, $`n/t=[nv(E)D_nn/x]/x+O(\overline{ϵ})`$. These results and Eq. (3) implies

$$\frac{E}{t}=\frac{1}{ϵ}eN_D^{}v(E)v(E)\frac{E}{x}+D_n\frac{^2E}{x^2}+\frac{1}{ϵ}J_{tot}(t),$$
(5)
$$\frac{p}{t}=gI(x,t)\gamma (N_D^{}+\frac{ϵ}{e}\frac{E}{x})p+\frac{}{x}\left(p\mu _pE+D_p\frac{p}{x}\right),$$
(6)
where $`J_{tot}(t)`$ is time-dependent current density. The physical meaning of Eq. (5) is quite clear, which represents the Gunn-domain formation and with drift velocity $`v(E_0)`$ . Then, Eq. (5) confirms our LSA results in Ref. . In the following we want to show the detailed structure of PDT.

In order to get the solution of Eq. (5), an ansatz of the electric field, we set $`E(x,t)=E_0+[E_s(x,t)e^{i\overline{K}(xv_0t)}+E_s^{}(x,t)e^{i\overline{K}(xv_0t)}]/2`$, where $`v_0=v(E_0)`$. $`E_s(x,t)`$ is a new variable corresponding to the amplitude of electric-field domains and $`E_s^{}(x,t)`$ is its complex conjugate; $`\overline{K}`$ ($`=2\pi /L`$) is a bulk property of semiconductor. Since $`\overline{K}`$ is small, the spatiotemporal behavior of the dynamical system within the linear dimension $`\mathrm{\Lambda }`$ is dominated by the behavior of $`E_s(x,t)`$. Actually, $`\overline{K}`$ is an irrelevent parameter which is verified by the numerical works of . When $`E_0|E_s|`$, we expand $`v(E)`$ as: $`v(E)v_0+v_0^{(1)}(EE_0)+v_0^{(2)}(EE_0)^2/2+v_0^{(3)}(EE_0)^3/6`$, where $`v_0^{(j)}=d^jv(E)/dE^j|_{E=E_0}`$. Equation (5) and above expansions for $`E(x,t)`$ and $`v(E)`$ imply

$`{\displaystyle \frac{E_s}{t}}`$ $`=`$ $`({\displaystyle \frac{e}{ϵ}}N_D^{}v_0^{(1)}D_n\overline{K}^2)E_s+(i2D_n\overline{K}v_0){\displaystyle \frac{E_s}{x}}{\displaystyle \frac{1}{8}}v_0^{(2)}\left(E_s^2{\displaystyle \frac{E_s^{}}{x}}+2|E_s|^2{\displaystyle \frac{E_s}{x}}\right)`$ (8)
$`+D_n{\displaystyle \frac{^2E_s}{x^2}}({\displaystyle \frac{1}{8}}{\displaystyle \frac{e}{ϵ}}N_D^{}v_0^{(3)}+i{\displaystyle \frac{1}{8}}v_0^{(2)}\overline{K})|E_s|^2E_s.`$
With $`\overline{K}`$ 0 (as length $`L\mathrm{}`$) and $`E_s\mathrm{\Psi }(t)e^{iKx}`$, Eq. (8) implies

$$\frac{d\mathrm{\Psi }}{dt}=C\mathrm{\Psi }\frac{1}{8}\left(\frac{e}{ϵ}N_D^{}v_0^{(3)}+iKv_0^{(2)}\right)|\mathrm{\Psi }|^2\mathrm{\Psi }.$$
(9)
Here $`C=(eN_D^{}v_0^{(1)}/ϵ+D_nK^2+iKv_0)`$. With $`\mathrm{\Psi }(t)=R(t)`$exp$`i\mathrm{\Theta }(t)`$, Eq. (9) implies $`R(t)=(a/b)^{1/2}[r_0^2e^{2at}/(r_0^2(1e^{2at})+a)]^{1/2}`$, $`\mathrm{\Theta }(t)=Kv_0t\frac{1}{8}Kv_0^{(2)}^tR^2(t^{})𝑑t^{}`$, where $`a=eN_D^{}v_0^{(1)}/ϵ+D_nK^2`$, $`b=eN_D^{}v_0^{(3)}/8ϵ`$, and $`r_0`$ is related to initial value of $`R(t)`$. When $`a>0`$ (i.e., the operating point is within the regime of positive differential mobility PDM), $`b>0`$, and $`t\mathrm{}`$, $`R(t)`$ and $`\mathrm{\Theta }(t)`$ tend to $`0`$ and $`Kv_0t`$, respectively. These results mean that the origin on the complex plane of $`\mathrm{\Psi }(t)`$ is a stable focus, which also implies electric-field profile in material is homogeneous, i.e., $`E(x,t)=E_0`$. When $`a<0`$ (i.e., the NDM regime), $`b>0`$, and $`t\mathrm{}`$, $`R(t)=(a/b)^{1/2}`$ and $`\mathrm{\Theta }(t)=K(v_0+av_0^{(2)}/8b)t`$. This solution describes a circular trajectory with counterclockwise angular frequency $`K(v_0+av_0^{(2)}/8b)`$, radius $`(a/b)^{1/2}`$, and transient response time $`|a|^1`$. The transition from stable focus to stable limit cycle is a supercritical Hopf bifurcation. Since $`\overline{K}`$ 0, $`E_s\mathrm{\Psi }(t)e^{iKx}`$ and $`\mathrm{\Psi }(t)=R(t)`$exp$`i\mathrm{\Theta }(t)`$ imply for $`a<0`$, $`b>0`$, and $`t\mathrm{}`$:

$$E(x,t)=E_0+(a/b)^{1/2}\mathrm{cos}[K(xv_0tav_0^{(2)}t/8b)],$$
(10)
which means that the PDT has a spatial period $`\mathrm{\Lambda }`$ and a traveling velocity $`v_0+av_0^{(2)}/8b`$. Since $`v_0`$ ($`10^7`$ cm/s) is much larger than $`|av_0^{(2)}/8b|`$, so we neglect the term $`av_0^{(2)}/8b`$ below.

In general, $`E_s`$ may be written as

$$E_s=\underset{q=1}{\overset{\mathrm{}}{}}\mathrm{\Psi }_q(t)e^{iqKx}=\underset{q=1}{\overset{\mathrm{}}{}}\left(\frac{a_q}{b}\right)^{1/2}\left[\frac{r_0^2e^{2a_qt}}{r_0^2(1e^{2a_qt})+a_q}\right]^{1/2}e^{iqK(xv_0t)}.$$
(11)
where $`a_q=eN_D^{}v_0^{(1)}/ϵ+D_nK^2q^2`$. However, Eq. (11) needs to be corrected when we consider the coupling between spatial modes. This is very important for self-organizing systems. For example, consider two spatial modes $`\mathrm{\Psi }_1e^{ikx}`$ and $`\mathrm{\Psi }_2e^{i2kx}`$ with $`a_2>0>a_1`$ and $`b>0`$. From equations for $`R(t)`$ and $`\mathrm{\Theta }(t)`$ , we know that $`|\mathrm{\Psi }_1|(a_1/b)^{1/2}`$ and $`|\mathrm{\Psi }_2|0`$ when $`t\mathrm{}`$. Then $`|\mathrm{\Psi }_2|/|\mathrm{\Psi }_1|`$ will become zero. This result implies that only unstable spatial modes (i.e., $`a_q<0`$) have contributions for pattern-forming systems. However, if we combine these two modes, i.e., $`E(x,t)=E_0+_{q=1}^2[\mathrm{\Psi }_q(t)e^{iqKx}+\mathrm{\Psi }_q^{}(t)e^{iqKx}]/2`$, and substitute it into Eq. (5), we have

$$\frac{d\mathrm{\Psi }_1}{dt}=(a_1+iKv_0)\mathrm{\Psi }_1\underset{m_1+m_2+\mathrm{}+m_s=1}{}\frac{\beta _{s,1}}{s!}\mathrm{\Psi }_{m_1}\mathrm{\Psi }_{m_2}\mathrm{}\mathrm{\Psi }_{m_s},$$
(12)
$$\frac{d\mathrm{\Psi }_2}{dt}=(a_2+i2Kv_0)\mathrm{\Psi }_2\underset{m_1+m_2+\mathrm{}+m_s=2}{}\frac{\beta _{s,2}}{s!}\mathrm{\Psi }_{m_1}\mathrm{\Psi }_{m_2}\mathrm{}\mathrm{\Psi }_{m_s},$$
(13)
where $`m_s=\pm 1`$ or $`\pm 2`$, $`\beta _{s,m}=eN_D^{}v_0^{(s)}/ϵ+imsv_0^{(s1)}K`$, and $`\mathrm{\Psi }_q^{}=\mathrm{\Psi }_q`$ . Based on the slaving principle , we know that the stable mode $`\mathrm{\Psi }_2`$ will be determined by the unstable mode $`\mathrm{\Psi }_1`$, and that $`\mathrm{\Psi }_2`$, and hence the ratio $`|\mathrm{\Psi }_2|/|\mathrm{\Psi }_1|r_{2,1}(K,N_D^{})`$, will be nonzero. For simplicity, in the following we assume that $`r_{q+1,q}`$ ($`|\mathrm{\Psi }_{q+1}|/|\mathrm{\Psi }_q|`$) is independent of $`q`$ and denote the number by $`r`$, which can be determined by experiments (see below). From $`E(x,t)=E_0+[E_s(x,t)+E_s^{}(x,t)]/2`$ and Eq. (11), for large $`t`$ we get

$$E(x,t)=E_0+\left(\frac{a_1}{b}\right)^{1/2}\frac{r+\mathrm{cos}[K(xv_0t)]}{1+r^22r\mathrm{cos}[K(xv_0t)]}.$$
(14)
Therefore, the maximum and minimum fields are equal to $`E_{max}=E_0+(a_1/b)^{1/2}(1r)^1`$ and $`E_{min}=E_0(a_1/b)^{1/2}(1+r)^1`$, respectively. According to Pockel’s effect, the difference in refractive index between $`E_{max}`$ and $`E_{min}`$ is $`\overline{n}^3\gamma _{41}(a_1/b)^{1/2}(1r^2)^1`$, where $`\overline{n}`$ is refractive index of pure material and $`\gamma _{41}`$ is the electro-optic coefficient ($`\overline{n}=3.3`$ and $`\gamma _{41}=1.43\times 10^{10}`$ cm/V for GaAs when the wavelength of the incident light is 1 $`\mu `$m). The profile of Eq. (14), shown in Fig. 1, is strongly $`r`$ dependent. If $`r`$ is near 0, $`E(x,t)`$ is approxmately equal to $`\mathrm{cos}[K(xv_0t)]`$. If $`r`$ is near 1, Eq. (14) shows an obvious shock-wave structure. In some senses, $`r`$ is related to the product $`N_D^{}\mathrm{\Lambda }`$ in Kroemer’s criterion for domain formation. To see this, we first consider $`r1`$. In this case, $`(|\mathrm{\Psi }_q||\mathrm{\Psi }_{q+1}|)/|\mathrm{\Psi }_q|`$, and hence $`(|a_2||a_1|)/|a_1|`$, is much smaller than 1 and there are many excited unstable Fourier modes. Thus, $`eN_D^{}v_0^{(1)}/ϵ`$ will be much larger than $`D_nK^2`$ and it follows from Eq. (11) that Kroemer’s criterion for sharp domain formation is $`a_1\mathrm{\Lambda }/v_0eN_D^{}v_0^{(1)}\mathrm{\Lambda }/(v_0ϵ)1`$, i.e. $`N_D^{}\mathrm{\Lambda }v_0ϵe^1/v_0^{(1)}`$ ($`10^{12}`$ cm-2). On the other hand, in the case $`r0`$, $`eN_D^{}v_0^{(1)}/ϵ`$ is of the same order as $`D_nK^2`$ and Kroemer’s criterion for the appearance of sinusoidal electric-field domain is $`a_1\mathrm{\Lambda }/v_0eN_D^{}v_0^{(1)}\mathrm{\Lambda }/(v_0ϵ)1`$, i.e. $`N_D^{}\mathrm{\Lambda }v_0ϵe^1/v_0^{(1)}`$ ($`10^{12}`$ cm-2). Therefore, Kroemer’s criterion appears in Eq. (14) via the representation of $`r`$ parameter. To determine the relation between $`r`$ and $`N_D^{}\mathrm{\Lambda }`$, we can measure $`(E_{max}E_0)/(E_0E_{min})`$ ($`=(1+r)/(1r)`$) at different values of $`N_D^{}\mathrm{\Lambda }`$. Then we will find that $`r`$ is a function of $`N_D^{}\mathrm{\Lambda }`$. In other words, $`r`$ can be determined by experiments.

In order to get hole distribution, we assume $`p(x,t)=p_0+[p_s(x,t)e^{i\overline{K}(x+\mu _pE_0t)}+p_s^{}(x,t)e^{i\overline{K}(x+\mu _pE_0t)}]/2`$, and substitute such $`p(x,t)`$ into Eq. (6) to find

$$\frac{p_s}{t}=(D_p\overline{K}^2+\gamma N_D^{})p_s+D_p\frac{^2p_s}{x^2}+(2i\overline{K}D_p+\mu _pE_0)\frac{p_s}{x}+S(x,t).$$
(15)
Here $`S(x,t)=2gI_0e^{i\overline{K}(x+\mu _pE_0t)}m\mathrm{cos}(Kx+\mathrm{\Omega }t)+p_0(\mu _p\gamma ϵ/e)(\frac{E_s}{x}+i\overline{K}E_s)e^{i\overline{K}(v_0+\mu _pE_0)t}.`$ Equation (15) is a linear partial differential equation with external driving function $`S(x,t)`$ including the MIP and multiple Gunn domains. Thus $`p_s`$ can be considered as a linear combination of a homogeneous solution $`p_s^h(x,t)`$ and an inhomogeneous solution $`p_s^i(x,t)`$. Since the characteristic length of $`p_s`$ is equal to $`\mathrm{\Lambda }`$ ($`L`$), we may treat $`\overline{K}`$ as zero in Eq. (15) and write $`p_s^h/t=\gamma N_D^{}p_s^h+D_p^2p_s^h/x^2+\mu _pE_0p_s^h/x`$,

$$\frac{p_s^i}{t}=\gamma N_D^{}p_s^i+D_p\frac{^2p_s^i}{x^2}+\mu _pE_0\frac{p_s^i}{x}+2gI_0m\mathrm{cos}(Kx+\mathrm{\Omega }t)+p_0(\mu _p\gamma \frac{ϵ}{e})\frac{E_s}{x}.$$
(16)
It is easy to find solution of $`p_s^h`$ and show that $`p_s^h=0`$ as $`t\mathrm{}`$. When $`E_s=_{q=1}^{\mathrm{}}r^{q1}(a_1/b)^{1/2}`$exp$`[iqK(xv_0t)]`$ is substituted into Eq. (16), we can find solution for $`p_s^i`$. Therefore, the hole distribution $`p(x,t)=p_0+Re(p_s^i)`$ is

$`p(x,t)={\displaystyle \frac{gI_0}{\gamma N_D^{}}}+{\displaystyle \frac{2gI_0m}{A}}[(D_pK^2+\gamma N_D^{})\mathrm{sin}(Kx+\mathrm{\Omega }t)(K\mu _pE_0\mathrm{\Omega })\mathrm{cos}(Kx+\mathrm{\Omega }t)]`$ (17)
$`D{\displaystyle \underset{q=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{qr^{q1}\{q(v_0+\mu _pE_0)\mathrm{cos}[qK(xv_0t)]+(\gamma N_D^{}+q^2K^2D_p)\mathrm{sin}[qK(xv_0t)]\}}{(\gamma N_D^{}+q^2K^2D_p)^2+q^2K^2(v_0+\mu _pE_0)^2}}.`$ (18)
Here $`A=(D_pK^2+\gamma N_D^{})^2+(K\mu _pE_0\mathrm{\Omega })^2`$ and $`D=gI_0K(\mu _p\gamma ϵ/e)(a_1/b)^{1/2}/\gamma N_D^{}`$. The pattern profiles of Eqs. (14) and (18) correspond to simple dynamics.

When $`I_0`$ is large and the higher order corrections in Eqs. (1)-(4) need to be included, then we have

$$E/xe(nN_D^{})/ϵ=\overline{ϵ}e/ϵp^{(0)},$$
(19)
$$\frac{n}{t}+\frac{}{x}\left[nv(E)D_n\frac{n}{x}\right]=\overline{ϵ}\left[gI(x,t)\gamma (N_D^{}+\frac{ϵ}{e}\frac{E^{(0)}}{x})p^{(0)}\right],$$
(20)
where $`p^{(0)}`$ and $`E^{(0)}`$ are solutions of Eqs. (18) and (14), respectively. These results imply that up to order $`O(\overline{ϵ})`$ we have

$$\frac{E}{t}+\frac{1}{ϵ}eN_D^{}v(E)+v(E)\frac{E}{x}D_n\frac{^2E}{x^2}\frac{1}{ϵ}J_{tot}(t)=\overline{S}(x,t),$$
(21)
where $`\overline{S}(x,t)=\overline{ϵ}e\left[p^{(0)}v(E^{(0)})\mu _pE^{(0)}p^{(0)}(D_p+D_n)p^{(0)}/x\right]/ϵ`$ and represents the electron-hole coupling. To understand the underlying physics of Eq. (21), we follow the same procedure as the case of simple dynamics, to get

$$\frac{d\mathrm{\Psi }}{dt}=C\mathrm{\Psi }\frac{1}{8}\left(\frac{e}{ϵ}N_D^{}v_0^{(3)}+iKv_0^{(2)}\right)|\mathrm{\Psi }|^2\mathrm{\Psi }+_0^\mathrm{\Lambda }\overline{S}(x,t)e^{iKx}𝑑x,$$
(22)
which is a generalization of Eq. (9) and represents a self-sustained oscillator (i.e., PDT) under the external driving force $`\overline{ϵ}F(t)/ϵ`$ ($`_0^\mathrm{\Lambda }\overline{S}(x,t)e^{iKx}𝑑x`$) and

$`F`$ $`(t)=\left[iE_0\mu _p{\displaystyle \frac{\pi }{K}}iv_0{\displaystyle \frac{\pi }{K}}(D_p+D_n)\pi \right]G_1e^{i\mathrm{\Omega }t}\left[E_0\mu _p{\displaystyle \frac{\pi }{K}}v_0{\displaystyle \frac{\pi }{K}}+i(D_p+D_n)\pi \right]G_2e^{i\mathrm{\Omega }t}`$ (23)
$``$ $`\left[E_0\mu _p{\displaystyle \frac{\pi }{K}}v_0{\displaystyle \frac{\pi }{K}}+i(D_p+D_n)\pi \right]G_3e^{iKv_0t}+\left[iE_0\mu _p{\displaystyle \frac{\pi }{K}}iv_0{\displaystyle \frac{\pi }{K}}(D_p+D_n)\pi \right]G_4e^{iKv_0t}`$ (24)
$`+`$ $`{\displaystyle \frac{\pi }{K}}\left({\displaystyle \frac{a_1}{b}}\right)^{1/2}\left(v_0^{(1)}{\displaystyle \frac{a_1}{8b}}v_0^{(3)}\mu _p\right)G_0e^{iKv_0t},`$ (25)
where $`G_0=gI_0/\gamma N_D^{},G_1=2gI_0m(D_pK^2+\gamma N_D^{})/A,G_2=2gI_0m(K\mu _pE_0\mathrm{\Omega })/A`$, $`G_3=gI_0K^2(\mu _p\gamma ϵ/e)(a_1/b)^{1/2}(v_0+\mu _pE_0)/B`$, $`G_4=gI_0K(\mu _p\gamma ϵ/e)(a_1/b)^{1/2}(\gamma N_D^{}+K^2D_p)/B`$, and $`B=\gamma N_D^{}(\gamma N_D^{}+K^2D_p)^2+K^2(v_0+\mu _pE_0)^2`$. The intrinsic oscillating frequency in Eq. (22) is $`v_0/\mathrm{\Lambda }`$ and the external driving frequencies include $`v_0/\mathrm{\Lambda }`$ and $`\mathrm{\Omega }/2\pi `$. If the operating point is in the PDM regime (i.e., $`a_1>0`$), the last three terms in Eq. (25) will be zero. This means that there is no Gunn-domain formation in the semiconductor, but the external optical waves, i.e., the first two terms on the right-hand side of Eq. (25), will trigger a space-charge field (SCF) $`e^{i(Kx+\mathrm{\Omega }t)}`$, which is well known in photorefractive material . On the other hand, if the operating point is in the NDM regime, in addition to the SCF, the PDT will be generated by the intrinsic electrical instability with the oscillating frequency $`Kv_0`$. The relative strength of the SCF and PDT depends on the strength of spatially independent part of the MIP, i.e. $`I_0`$. When $`I_0`$ is small, SCF may be neglected and the behavior of the system is dominated by the PDT. When $`I_0`$ is large, the SCF will compete with the PDT and the system will show complex dynamics. For example, these two competing frequencies $`\mathrm{\Omega }/2\pi `$ and $`v_0/\mathrm{\Lambda }`$ can lead to quasiperiodic behavior.

The picture presented above is useful for understanding the results of other groups. The system numerically studied by Subačius et al.  may be represented by Eqs.(1)-(4) of this Communication with $`N_D^{}=0`$ and $`\mathrm{\Omega }=0`$. When the system is illuminated by a pulse of light interference field with intensity $`I(x,t)=I_0f(t)[1+m\mathrm{cos}(Kx)]`$ and $`m<<1`$ ($`f(t)`$ is the temporal shape of a laser pulse), the temporal part $`I_0f(t)`$ generates background electrons in the conduction band and the very small spatially dependent part of the light makes a small modulation on the background. The situation is similar to the system of the present manuscript with very small $`I_0`$. Fig. 1 in the present Communication is very similar as the numerical results of Subačius et al.. On the other hand, Segev et al.  and Bonilla et al.  considered deep-impurity doped GaAs under the externally optical waves, i.e., $`I(x,t)=I_0[1+m\mathrm{cos}(Kx+\mathrm{\Omega }t)]`$. In their case the modulation depth $`m`$ also plays a crucial role. For example, Segev et al. considered $`m`$ as a small parameter. Then the dynamical variables can be analytically expanded to Fourier series based on the form of $`I(x,t)`$. Therefore they can get multiple Gunn-domain order by order. Actually the physical idea of Segev et al. is the same as Subačius et al.. The spatial period of $`I(x,t)`$ determines the spacing of Gunn domains and $`I_0`$ (or $`I_0f(t)`$) generates background electrons. However the work of Segev et al. didn’t discuss the stability of multiple Gunn-domain, which will lead to the travelling domains with the same moving velocity (i.e., $`\mathrm{\Omega }/K`$) of $`I(x,t)`$. Furthermore, according to the prediction of Segev et al. the moving velocity of electric-field profile in PDM or NDM regime has no any difference. This unsuitable prediction can be addressed and verified by Eq. (20). As for the numerical work of Bonilla et al., they showed complicated spatiotemporal behaviors of Gunn domains, e.g., quasiperiodic route to chaos, are mostly observed in their simulations. In other words, PDT didn’t appear in their system. Therefore they concluded that Kroemer’s criterion can not give the correct numbers of Gunn domains. We think that $`m=0.1`$ used in their study is large enough to generate driven chaos due to the competition between SCF and PDT. Therefore both considerations of Kroemer’s criterion and pattern competition may give correct physics in light-triggered Gunn-domain systems.

In conclusion, using a nonlinear analytic method we confirm the LSA results of  and find that the competition between the SCF and the PDT could be a dominant factor to determine the spatiotemporal pattern in semiconductor. We hope that the formulation of the present Communication can be further verified by numerical calculations or experiments.

We thank Jonathan Dushoff for a critical reading of the paper. This work was supported in part by the National Science Council of the Republic of China (Taiwan) under Contract No. NSC 89-2112-M-001-084.

Figure Caption

FIG. 1. The profile of $`E(x,t)`$ represented by Eq. (14) for several values of $`r`$, where $`E_0=10`$ kV/cm, $`(a_1/b)^{1/2}=9.8`$ kV/cm, $`\mathrm{\Lambda }=15`$ $`\mu `$m, and $`L=45`$ $`\mu `$m.