text:
# 1 Introduction

## 1 Introduction

Nonequilibrium phenomena in real systems are usually described by sets of nonlinear differential or integro-differential equations in partial derivatives. The solutions to such equations are in many cases nonuniform in space exhibiting the formation of various spatial structures. It often happens that a given set of equations possesses several solutions corresponding to different spatial patterns. And, moreover, it is often impossible to distinguish between these solutions being based on stability analysis since all these solutions can be stable. If an ensemble of solutions can be parametrized by a multiparameter $`\beta `$ from a manifold $``$, so that all solutions corresponding to any $`\beta `$ are stable, then the manifold $``$ is called the stability balloon .

In the case of such a nonuniqueness of solutions one may wander which of the multiplicity of allowed states would actually be found for a given experimental protocol. This is how the problem of pattern selection arises, when the considered equations have many solutions for given external conditions, whereas observed patterns may constitute a more restricted set. And there no doubts are numerous systems and phenomena in real life with such a complicated behaviour .

Not only the solutions to equations can be nonunique but real phenomena may also display a variety of patterns for the same given conditions. This means that the multiplicity of solutions is not just an artifact but a feature expressing the intrinsic complexity of natural phenomena. If real patterns show multiplicity, is there any ordering between them, such that one solution is preferred over the other? In equilibrium thermodynamics this is a familiar concept, with the free energy providing the ordering principle. But in nonequilibrium systems, there is no such a general organizing principle to apply . For some particular cases several heuristic recipes have been suggested. For example, one state may be preferred over the other if it has a larger basin of attraction for typical initial conditions. Or one may think that the fastest growing mode dominates the evolution. Or selection is made assuming that the pattern with the fastest spatial decay rate is preferable. The limits of such heuristic arguments are discussed in detail in Cross and Hehenberg  where one can find extensive citations to literature on the problem.

In the present paper a general approach is developed for treating the problem of pattern selection. The approach, which is formulated in Sec. 2, is based on defining the probability distribution of patterns. In Sections 4 to 5, this approach is applied to the problem of turbulent photon filamentation. An introduction to the latter problem is given in Section 3. The obtained results are in very good agreement with experiment. The final Section 7 contains conclusions.

## 2 Probability Distribution of Patterns

Let us consider a system of evolution equations, which displays the multiplicity of solutions describing different spatial structures. Assume that these solutions can be parametrized by a multiparameter $`\beta `$. All admissible values of $`\beta `$ form a manifold $`=\{\beta \}`$. In many cases, differential equations in partial derivatives can be reduced to a $`d`$-dimensional system of ordinary differential equations, with a dimension $`d`$ that may equal infinity . For the simplicity of notation, we shall keep in mind this possibility of working with a $`d`$-dimensional dynamical system. A generalization of this case will be given at the end of the present Section.

The state of a $`d`$-dimensional dynamical system is the set

$$y(t)=\{y_i(t)=y_i(\beta ,t)|i=1,2,\mathrm{},d\}$$
(1)
of solutions to the system of differential equations, which can always be presented in the normal form

$$\frac{d}{dt}y(t)=v(y,t),$$
(2)
where $`v=\{v_i|i=1,2,\mathrm{},d\}`$ is a velocity field. By assumption, the multiplicity of solutions is parametrized by a multiparameter $`\beta `$, but in intermediate calculations we shall not, for the compactness of notation, label functions with $`\beta `$, restoring this dependence in final formulas.

Our aim is to classify the states (1) labelled by $`\beta `$ by defining a probability measure on the manifold $``$. To introduce the probability distribution $`p(\beta ,t)`$ of patterns labelled by a multiparameter $`\beta `$, we may resort to the ideas of statistical mechanics , where a probability $`p`$ can be connected with entropy $`S`$ by the relation $`pe^S`$. In the nonequilibrium case, the entropy is a function of time, $`S(t)`$. Since it is not the entropy itself but rather its change that is measurable, it is convenient to count entropy from its initial value $`S(0)`$, thus, considering the entropy variation

$$\mathrm{\Delta }S(t)=S(t)S(0),$$
(3)
which is a kind of relative entropy . The entropy may be defined as the logarithm of an elementary phase volume , which for the nonequilibrium case can be expressed as

$$S(t)=\mathrm{ln}|\delta \mathrm{\Gamma }(t)|,$$
(4)
with the elementary phase volume

$$\delta \mathrm{\Gamma }(t)\underset{i}{}\delta y_i(t).$$
(5)
Hence the entropy variation (3) is

$$\mathrm{\Delta }S(t)=\mathrm{ln}\left|\frac{\delta \mathrm{\Gamma }(t)}{\delta \mathrm{\Gamma }(0)}\right|.$$
(6)
The probability distribution $`pe^{\mathrm{\Delta }S}`$, being normalized with respect to $`\beta `$, takes the form

$$p(\beta ,t)=\frac{e^{\mathrm{\Delta }S(\beta ,t)}}{Z(t)},Z(t)e^{\mathrm{\Delta }S(\beta ,t)}𝑑\beta ,$$
(7)
where the integration over $`\beta `$ runs through the manifold $``$. If the latter manifold is not continuous, the integration is to be replaced by summation.

The elementary phase volume (5) can be presented as

$$\delta \mathrm{\Gamma }(t)=\underset{i}{}\underset{j}{}M_{ij}(t)\delta y_j(0),$$
(8)
where $`M_{ij}(t)`$ are the elements of the multiplier matrix  defined through the variational derivatives

$$M_{ij}(t)\frac{\delta y_i(t)}{\delta y_j(0)},M_{ij}(0)=\delta _{ij}.$$
(9)
Then the entropy variation (6) writes

$$\mathrm{\Delta }S(t)=\underset{i}{}\mathrm{ln}|M_{ii}(t)|.$$
(10)
The multiplier matrix $`\widehat{M}(t)=[M_{ij}(t)]`$ satisfies the equation

$$\frac{d}{dt}\widehat{M}(t)=\widehat{J}(y,t)\widehat{M}(t),$$
(11)
which follows from the variation of the evolution equation (2) and where $`\widehat{J}(y,t)=[J_{ij}(y,t)]`$ is the Jacobian matrix with the elements

$$J_{ij}(y,t)=\frac{\delta v_i(y,t)}{\delta y_j(t)}.$$
(12)
The initial condition for Eq. (11), according to the definition (9), is $`\widehat{M}(0)=[\delta _{ij}]`$. Combining Eqs. (7) and (10), we get the probability distribution of patterns,

$$p(\beta ,t)=\frac{1}{Z(t)}\underset{i}{}\frac{1}{|M_{ii}(\beta ,t)|},$$
(13)
with the normalization factor

$$Z(t)=\underset{i}{}\frac{1}{|M_{ii}(\beta ,t)|}d\beta .$$

This is the probability distribution of solutions classified with a multiparameter $`\beta `$. Since each solution, by definition, represents a particular pattern, expression (13) is the probability distribution of patterns. This expression naturally connects the notion of probability and the notion of stability. Really, the multipliers are smaller by modulus for more stable solutions and, consequently, for more probable patterns.

To calculate the probability distribution of patterns (13), we need to know the multipliers (9) which are defined by the equation (11). Using the latter equation, the pattern distribution (13) can be transformed as follows. Introduce the matrix $`\widehat{L}(t)=[L_{ij}(t)]`$ with the elements

$$L_{ij}(t)\mathrm{ln}|M_{ij}(t)|.$$
(14)
Then the entropy variation (10) writes

$$\mathrm{\Delta }S(t)=\mathrm{Tr}\widehat{L}(t).$$
(15)
The trace of a matrix does not depend on the matrix representation. Hence, we may perform intermediate calculations using one particular representation, returning at the end to the form independent of a representation. To this end, let us consider a representation when the multiplier matrix is diagonal. Because of Eq. (11) with the initial condition $`M_{ij}(0)=\delta _{ij}`$, the multiplier matrix is diagonal if an only if the Jacobian matrix is diagonal too. In this case, the evolution equation (11) yields

$$M_{ii}(t)=\mathrm{exp}\left\{_0^tJ_{ii}(y(t^{}),t^{})𝑑t^{}\right\}.$$

From Eq. (14) we have

$$L_{ii}(t)=_0^t\mathrm{Re}J_{ii}(y(t^{}),t^{})𝑑t^{}.$$

Introducing the notation

$$K(t)\underset{i}{}\mathrm{Re}J_{ii}(y,t),$$
(16)
we get

$$\mathrm{Tr}\widehat{L}(t)=_0^tK(t^{})𝑑t^{}.$$

Without the loss of generality, we may assume that the state (1) consists of real functions since any complex function can always be treated as a pair of real functions. Hence the velocity field in Eq. (2) can also be considered as real. Then the eigenvalues of the Jacobian matrix (12) are either real or, if complex, come in complex conjugate pairs. Therefore

$$\underset{i}{}\mathrm{Re}J_{ii}(y,t)=\underset{i}{}J_{ii}(y,t)=\mathrm{Tr}\widehat{J}(y,t).$$

Thus, Eq. (16) can be written as

$$K(t)=\mathrm{Tr}\widehat{J}(y,t).$$
(17)
And for the entropy variation (15), we find

$$\mathrm{\Delta }S(t)=_0^tK(t^{})𝑑t^{}.$$
(18)
Expression (17) in dynamical theory is called the contraction rate . All consideration given above can be straightforwardly generalized to the case when the state (1) consists of functions $`y_i(x,t)`$, where $`x`$ is the set of space variables. This results not more than in a slight complication of notations. Then everywhere the variable $`x`$ appears together with the index $`i`$ as an additional continuous index, and the sums over $`i`$ are to be complimented by the integrals over $`x`$. The multiplier matrix (9) becomes a matrix with respect to indices $`i,j`$ as well as with respect to $`x,x^{}`$,

$$M_{ij}(x,x^{},t)\frac{\delta y_i(x,t)}{\delta y_j(x^{},0)},$$

with the initial condition

$$M_{ij}(x,x^{},0)=\delta _{ij}\delta (xx^{}).$$

Similarly, the Jacobian matrix (12) becomes a matrix with the elements

$$J_{ij}(x,x^{},y,t)\frac{\delta v_i(x,y,t)}{\delta y_j(x^{},t)}.$$

Employing the matrix notation with respect to continuous variables , we may repeat the same steps as above in deriving the probability distribution (13). As is mentioned, the sum over $`x`$ is to be understood as the corresponding integral. And the product over a continuous variable can be defined  as

$$\underset{x}{}f(x)\mathrm{exp}\mathrm{ln}f(x)𝑑x.$$

As a result, the pattern distribution (13) reads

$$p(\beta ,t)=\frac{1}{Z(t)}\mathrm{exp}\left\{\underset{i}{}\mathrm{ln}|M_{ii}(x,x,\beta ,t)|dx\right\}.$$

The contraction rate (17) becomes

$$K(t)=\underset{i}{}J_{ii}(x,x,y,t)𝑑x,$$
(19)
where $`y=\{y_i(x,t)\}`$. In this way, the contraction rate has always the form of $`\mathrm{Tr}\widehat{J}(y,t)`$, as in Eq. (17), where the trace has to be defined according to the representation of the Jacobian matrix.

Restoring in the contraction rate the dependence on the parameter $`\beta `$ labelling different patterns and using the entropy variation (18), we finally obtain the probability distribution of patterns

$$p(\beta ,t)=\frac{1}{Z(t)}\mathrm{exp}\left\{_0^tK(\beta ,t^{})𝑑t^{}\right\},$$
(20)
in which the normalization factor is

$$Z(t)=\mathrm{exp}\left\{_0^tK(\beta ,t^{})𝑑t^{}\right\}𝑑\beta .$$

Thus, each solution labelled by $`\beta `$ is equipped with the probability weight (20). Consequently, that pattern is preferred over the other which has a higher probability weight. This is equivalent, because of the form (20), to saying that one pattern is preferable over others if its local contraction

$$\mathrm{\Lambda }(\beta ,t)\frac{1}{t}_0^tK(\beta ,t^{})𝑑t^{}$$

is minimal with respect to $`\beta `$. The latter provides the ordering principle for classifying solutions and related patterns. The local contraction plays for nonequilibrium dynamical systems the role analogous to that of the free energy for equilibrium statistical systems.

## 3 Turbulent Photon Filamentation

To illustrate the probabilistic approach to pattern selection, developed in the previous section, let us consider spatial structures appearing in resonant samples when increasing the Fresnel number. In active media interacting with electromagnetic field, there appears a variety of spatiotemporal structures. For example, electric field in laser cavities can exhibit spatial states, such as solitons and vortices, which bear a close analogy with similar structures in liquids . In general, there is a direct correspondence between the Maxwell-Bloch equations for slowly varying field amplitudes and hydrodynamic equations for compressible viscous liquid . The Fresnel number for optical systems plays the same role as the Reynolds number for fluids. In the same way as when increasing the Reynolds number, the fluid becomes turbulent, the field dynamics exhibits chaotic behaviour when increasing the Fresnel number. Similarly to hydrodynamics, optical systems display spatial multi-stability with several coexisting distinct stable states for the same values of parameters . Spatiotemporal chaos in optical systems is characterized by the same fast decay of correlation functions as it occurs for turbulent fluids. This is why one calls such chaotic optical phenomena optical turbulence .

In active optical systems, having the standard cylindrical shape, the transition from regular behaviour to spatiotemporal chaos occurs with the increasing Fresnel number

$$F\frac{\pi R^2}{\lambda L},$$

where $`R`$ is the cylinder radius, $`L`$ is the characteristic length, and $`\lambda `$ is the radiation wavelength. Physical processes accompanying the route from a regular regime to chaotic one are similar in different active optical media.

At small Fresnel numbers $`F1`$, there exists the sole transverse central mode practically uniformly filling the medium. When the Fresnel number is around $`F1`$, the cavity can house several transverse modes seen as a regular arrangement of bright spots in the transverse cross-section. These regular spatial structures emerge from an initially homogeneous state with a break of space-translational symmetry. They are regular in space forming ordered geometric arrays, such a polygons, and they are regular in time being either stationary or periodically oscillating. These transverse structures are imposed by the cavity geometry and correspond to the empty-cavity Gauss-Laguerre modes. Such regular structures are well understood theoretically, their description being based on field expansions over the modal Gauss-Laguerre functions prescribed by the cylindrical geometry, and the theory being in reasonable agreement with experiments for lasers, e.g. for CO2 and Na2 lasers \[11–17\], and for active nonlinear media, as the photorefractive Bi12SiO20 crystal pumped by a laser \[18–20\]. Similar structures also arise in many passive nonlinear media, such as a Kerr medium . For Fresnel numbers up to $`F5`$, the number of bright spots is proportional to $`F^2`$. In the longitudinal cross-section this corresponds to the existence of bright filaments whose number follows the $`F^2`$ law as $`F`$ increases.

Around $`F10`$ there occurs a principal change of properties, from the existence of regular structures to a turbulent-type state \[21–23\], with the intermittent behaviour in the region $`5<F<15`$. This happens in lasers \[21–23\] as well as in photorefractive crystals \[18–20\]. For Fresnel numbers $`F>15`$, the arising spatial structures are very different from those associated with the empty-cavity modes. The modal expansion is no longer relevant and the boundary conditions have no importance. The medium looks as consisting of a large number of parallel independently flashing filaments, whose number is proportional to $`F`$, contrary to the case of small Fresnel numbers with the number of filaments proportional to $`F^2`$. The filaments are chaotically distributed in space, are not correlated with each other, and are aperiodically flashing in time. Such a spatio-temporal chaotic behaviour is characteristic of hydrodynamic turbulence, this is why the similar phenomenon in optics is commonly called the optical turbulence . Contrary to the regime of small $`F`$, when the regular spatial structures are prescribed by the geometry and boundary conditions imposing their symmetry constraints, the turbulent optical filamentation is strictly self-organized, with its organization emerging from intrinsic properties of the medium . This kind of optical turbulence has been observed in both photorefractive crystals \[18–20\] and lasers \[21–23,25–30\]. Especially accurate and thorough experimental studies for CO2 and Dye lasers have been accomplished in Refs. \[25–30\]. The independence of turbulent optical filamentation of boundary conditions and, hence, its purely self-organized nature are confirmed by the observation of this effect in the resonatorless discharge-tube superluminescent samples, such as lasers on Ne, Tl, Pb, N2, and N$`{}_{}{}^{+}{}_{2}{}^{}`$ vapors \[31–35\]. Since the optical turbulence is characterized by the formation of bright filaments with a high density of photons, this phenomenon may be named the turbulent photon filamentation.

Let us summarize the main features characterizing the regular and turbulent regimes in active optical media. At low Fresnel numbers $`F<10`$, the regular regime occurs whose typical characteristics are: (i) Bright filaments are regularly arranged in space forming an ordered structure seen in the transverse cross-section as a polygon made of bright spots. (ii) The spatially ordered structure is either stationary or periodically oscillating in time. (iii) The number of bright filaments is proportional to $`F^2`$.

At high Fresnel number $`F>10`$, the turbulent regime develops whose typical features are: (i) The bright photon filaments are chaotically distributed in space. (ii) The filaments aperiodically flash in time. (iii) The number of these uncorrelated filaments is proportional to $`F`$.

While for the low-Fresnel-number regime, with regular optical structures, a good theoretical understanding was achieved \[11–20\], for the high-Fresnel-number regime of the turbulent photon filamentation, there has been no persuasive theory developed. This problem was addresses in Refs. \[36–39\], where the consideration was based on a simple model, only the stationary case was analysed, and the minimum-energy arguments were employed. This model consideration showed that photon filaments can really be formed in resonant media due to an effective interaction between atoms through photon exchange. The estimates for the filament radius turned out to be in good agreement with experiment. However, the empirical approach of Refs. \[36–39\] cannot be accepted as completely satisfactory. This is mainly because of the two principal points. One is that a stationary case was addressed, while turbulence is rather a nonstationary phenomenon and must be treated by time-dependent evolution equations. The second point is that the minimal-energy condition was invoked, while there is no such a principle for nonequilibrium systems.

Here we demonstrate that the problem of turbulent photon filamentation can be successfully treated by the probabilistic approach to pattern selection formulated in Section 2. Our consideration is based on realistic evolution equations for resonant media. Such equations are necessary for calculating the contraction rate $`K`$ defining the probability distribution of patterns (20).

## 4 Equations for Resonant Medium

The system of $`N`$ resonant atoms interacting with electromagnetic field is described by the Hamiltonian

$$\widehat{H}=\widehat{H}_a+\widehat{H}_f+\widehat{H}_{af},$$
(21)
consisting of the following terms: The Hamiltonian of two-level resonant atoms

$$\widehat{H}_a=\frac{1}{2}\underset{i=1}{\overset{N}{}}\omega _0(1+\sigma _i^z),$$
(22)
with the transition frequency $`\omega _0`$, where $`\sigma _i^z`$ is a Pauli operator. The radiation-field Hamiltonian

$$\widehat{H}_f=\frac{1}{8\pi }\left(\stackrel{}{E}^2+\stackrel{}{H}^2\right)𝑑\stackrel{}{r},$$
(23)
with electric field $`\stackrel{}{E}`$ and magnetic field $`\stackrel{}{H}=\stackrel{}{}\times \stackrel{}{A}`$, where the vector potential $`\stackrel{}{A}`$ is assumed to satisfy the Coulomb gauge calibration $`\stackrel{}{}\stackrel{}{A}=0`$. The atom-field interaction Hamiltonian

$$\widehat{H}_{af}=\underset{i=1}{\overset{N}{}}\left(\frac{1}{c}\stackrel{}{j}_j\stackrel{}{A}_i+\stackrel{}{d}_i\stackrel{}{E}_{0i}\right)$$
(24)
has the standard dipole form, in which $`\stackrel{}{A}_i\stackrel{}{A}(\stackrel{}{r}_i,t)`$; the transition-current and transition dipole operators are

$$\stackrel{}{j}_i=i\omega _0\left(\stackrel{}{d}\sigma _i^+\stackrel{}{d}^{}\sigma _i^{}\right),\stackrel{}{d}_i=\stackrel{}{d}\sigma _i^++\stackrel{}{d}^{}\sigma _i^{},$$

where $`\stackrel{}{d}`$ is the transition dipole and $`\sigma _i^\pm `$ are the rising or lowering operators, respectively; and $`\stackrel{}{E}_{0i}`$ is a seed field .

To consider nonuniform systems with arising spatial structures, it is convenient to invoke the space-time representation of evolution equations, whose general description is given in the book  and a detailed analysis in Refs. . The equations are written for the average quantities

$$u(\stackrel{}{r},t)<\sigma ^{}(\stackrel{}{r},t)>,s(\stackrel{}{r},t)<\sigma ^z(\stackrel{}{r},t)>,$$
(25)
using the standard semiclassical and Born approximations. For the compactness of presentation, let us introduce the notation

$$f(\stackrel{}{r},t)f_0(\stackrel{}{r},t)+f_{rad}(\stackrel{}{r},t)$$
(26)
for an effective field acting on an atom. This field consists of the term

$$f_0(\stackrel{}{r},t)i\stackrel{}{d}\stackrel{}{E}_0(\stackrel{}{r},t),$$
(27)
due to the cavity seed field, and of the term

$$f_{rad}(\stackrel{}{r},t)=\frac{3}{4}i\gamma \rho \left[\phi (\stackrel{}{r}\stackrel{}{r}^{})u(\stackrel{}{r}^{},t)\stackrel{}{e_d}^2\phi ^{}(\stackrel{}{r}\stackrel{}{r}^{})u^{}(\stackrel{}{r}^{},t)\right]𝑑\stackrel{}{r}^{}$$
(28)
corresponding to the interaction of atoms through the common radiation field, where $`\rho `$ is the density of atoms, $`\stackrel{}{d}d_0\stackrel{}{e}_d,d_0|\stackrel{}{d}|`$, and

$$\phi (\stackrel{}{r})\frac{\mathrm{exp}(ik_0|\stackrel{}{r}|)}{k_0|\stackrel{}{r}|},k_0\frac{\omega _0}{c},\gamma \frac{4}{3}k_0^3d_0^2.$$

The seed field

$$\stackrel{}{E}_0(\stackrel{}{r},t)=\stackrel{}{E}_1e^{i(kz\omega t)}+\stackrel{}{E}_1^{}e^{i(kz\omega t)}$$
(29)
selects a longitudinal mode propagating along the axis $`z`$ which is the axis of a cylindrical sample. The resulting equations are

$$\frac{du}{dt}=(i\omega _0+\gamma _2)u+sf,$$

$$\frac{ds}{dt}=2(u^{}f+f^{}u)\gamma _1(s\zeta ),$$
(30)
$$\frac{d|u|^2}{dt}=2\gamma _2|u|^2+s(u^{}f+f^{}u),$$

where $`\gamma _1`$ and $`\gamma _2`$ are the level and line widths, respectively, and $`\zeta `$ is the pumping parameter defined by the intensity of nonresonant lamp pumping.

The sample has cylindrical shape typical of lasers. The radiation wavelength $`\lambda `$, cylinder radius $`R`$, and its length $`L`$ are related by the inequalities

$$\frac{\lambda }{R}1,\frac{R}{L}1.$$
(31)
We consider the quasiresonance situation when

$$\frac{|\mathrm{\Delta }|}{\omega _0}1,\mathrm{\Delta }=\omega \omega _0.$$
(32)
As usual, the relaxation parameters are assumed to be small, so that

$$\frac{\gamma _1}{\omega _0}1,\frac{\gamma _2}{\omega _0}1.$$
(33)
The solutions to Eqs. (30) are, in general, nonuniform. This is because atoms interact with each other thorough the radiation field (28) containing the function $`\phi (\stackrel{}{r})`$ fastly oscillating in space and diminishing with increasing $`|\stackrel{}{r}|`$. Hence the atoms, situated far from each other, do not interact and radiate independently, although in the longitudinal direction they may be correlated due to the cylindrical geometry of the sample and the propagation of radiated field along the axis $`z`$. This suggests to look for a solution of Eqs. (30) in the form of a bunch of filaments aligned along the cylinder axis. Let us imagine that the radiation field inside the sample is stratified into $`N_f`$ photon filaments stretched along the axis $`z`$. These filaments can have different radii $`r_n`$, with $`n=1,2,\mathrm{},N_f`$. The radiation is mainly concentrated inside the bright filaments, fading away outside them, so that at the radial distance $`b_n`$ the corresponding solutions are an order of magnitude smaller than at the axis of a filament. The relation between the characteristic lengths $`b_n`$ and $`r_n`$ can be written down if the profile of solutions inside a filament is known. If the filament profile can be approximated by the normal law $`\mathrm{exp}(r^2/2r_n^2)`$, with the filament radius $`r_n`$ being the standard deviation, then

$$b_n=\sqrt{2\mathrm{ln}10}r_n.$$
(34)
The filaments do not interact with each other, because of which their locations in the radial cross-section are random. This means that the filament axes are located at random points $`\{x_n,y_n\}`$, where $`n=1,2,\mathrm{},N_f`$.

The above discussion explains why it is reasonable to try to find the solutions of Eqs. (30) in the form of expansions

$$u(\stackrel{}{r},t)=\underset{n=1}{\overset{N_f}{}}u_n(\stackrel{}{r},t)\mathrm{\Theta }_n(x,y)e^{ikz},s(\stackrel{}{r},t)=\underset{n=1}{\overset{N_f}{}}s_n(\stackrel{}{r},t)\mathrm{\Theta }_n(x,y)$$
(35)
over filaments, where

$$\mathrm{\Theta }_n(x,y)\mathrm{\Theta }\left(b_n^2(xx_n)^2+(yy_n)^2\right)$$

is a unit-step function. Substituting the presentation (35) in Eqs. (30), we obtain a system of equations for $`u_n,s_n`$, and $`|u_n|^2`$. These equations can be simplified if we are not interested in the detailed internal structure of each filament but rather wish to find out their characteristic sizes. Then we may employ the mean-field approximation for the averages

$$u(t)\frac{1}{V_n}_{V_n}u_n(\stackrel{}{r},t)𝑑\stackrel{}{r},s(t)\frac{1}{V_n}_{V_n}s_n(\stackrel{}{r},t)𝑑\stackrel{}{r},$$

$$|u(t)|^2\frac{1}{V_n}_{V_n}|u_n(\stackrel{}{r},t)|^2𝑑\stackrel{}{r},$$
(36)
where the integration is over the volume $`V_n\pi b_n^2L`$ and the index $`n`$ in the left-hand sides, for short, is dropped.

To present the equations for the functions (36) in a compact form, we introduce two effective coupling parameters:

$$g\frac{3\gamma \rho }{4\gamma _2V_n}_{V_n}\frac{\mathrm{sin}[k_0|\stackrel{}{r}\stackrel{}{r}^{}|k(zz^{})]}{k_0|\stackrel{}{r}\stackrel{}{r}^{}|}𝑑\stackrel{}{r}𝑑\stackrel{}{r}^{}$$
(37)
and

$$g^{}\frac{3\gamma \rho }{4\gamma _2V_n}_{V_n}\frac{\mathrm{cos}[k_0|\stackrel{}{r}\stackrel{}{r}^{}|k(zz^{})]}{k_0|\stackrel{}{r}\stackrel{}{r}^{}|}𝑑\stackrel{}{r}𝑑\stackrel{}{r}^{}.$$
(38)
These parameters enter the definitions of the collective frequency and collective width, respectively,

$$\mathrm{\Omega }\omega _0+g^{}\gamma _2s,\mathrm{\Gamma }\gamma _2(1gs).$$
(39)
Then from Eqs. (30), employing the introduced notation, for the functions (36) we obtain the equations

$$\frac{du}{dt}=(i\mathrm{\Omega }+\mathrm{\Gamma })uis\stackrel{}{d}\stackrel{}{E}_1e^{i\omega t},$$

$$\frac{ds}{dt}=4g\gamma _2|u|^2\gamma _1(s\zeta )4\mathrm{Im}\left(u^{}\stackrel{}{d}\stackrel{}{E}_1e^{i\omega t}\right),$$
(40)
$$\frac{d|u|^2}{dt}=2\mathrm{\Gamma }|u|^2+2s\mathrm{Im}\left(u^{}\stackrel{}{d}\stackrel{}{E}_1e^{i\omega t}\right),$$

describing the time evolution of a filament. Such a reduction of initial equations for a bunch of filaments to the set (40) of equations for each of the filaments has become possible due to the absence of interactions between filaments. Different filaments can have different radii whose distribution is yet to be defined.

It is worth emphasizing that the presentation of solutions as the expansions (35) over filaments is rather general and includes as well the case when there are no separate filaments at all. Really, if it turns out that the most probable filament radius is close to the radius of the sample, $`R`$, this would actually mean that there are no several filaments, but the whole volume of the sample is filled by one filament.

## 5 Characteristics of Arising Filaments

To analyze Eqs. (40), we employ the scale separation approach  which is a generalization of the averaging technique  to the case of nonequilibrium statistical systems. To this end, we notice first of all that the functional variables in Eqs. (40) can be classified onto fast and slow in time. Due to the existence of the small parameters (33), the variable $`u`$ is fast as compared to the slow variables $`s`$ and $`|u|^2`$. The slow variables play the role of quasi-invariants for the fast function $`u`$. With $`s`$ being a quasi-invariant, the solution for the fast variable $`u`$, described by the first of Eqs. (40), writes

$$u(t)=u_0e^{(i\mathrm{\Omega }+\mathrm{\Gamma })t}+\frac{s\stackrel{}{d}\stackrel{}{E}_1}{\omega \mathrm{\Omega }+i\mathrm{\Gamma }}\left[e^{i\omega t}e^{(i\mathrm{\Omega }+\mathrm{\Gamma })t}\right],$$
(41)
where $`u_0=u(0)`$. This solution is to be substituted in the second and third of Eqs. (40) whose right-hand sides are to be averaged over time explicitly entering the fastly oscillating functions. In this way, we meet the quantity

$$\alpha \frac{\mathrm{Im}}{s\mathrm{\Gamma }}\underset{\tau \mathrm{}}{lim}\frac{1}{\tau }_0^\tau u^{}(t)\stackrel{}{d}\stackrel{}{E}_1e^{i\omega t}𝑑t$$
(42)
describing the influence of the seed field on atoms. Taking into account Eq. (41), we have

$$\alpha =\frac{|\stackrel{}{d}\stackrel{}{E}_1|^2}{(\omega \mathrm{\Omega })^2+\mathrm{\Gamma }^2}.$$
(43)
Recall that the seed field is, by definition, a weak cavity field whose role is to fix the axis of the cylindrical sample. This field does not pump the energy into the system which is standardly done by nonresonant lamps and which is described by the pumping parameter $`\zeta `$ in the evolution equations (30) or (40). The seed field selects a longitudinal mode imposing no restrictions on transverse modes. The selection of the latter has to be done by the internal properties of the evolution equations. The weakness of the seed field is to be understood in the sense of the smallness of the parameter (43),

$$\alpha 1.$$
(44)
In analysing further the evolution equations (40), it is convenient to make the transformation of the functional variable $`|u|^2`$ to

$$w|u|^2\alpha s^2.$$
(45)
Finally, following the scale separation approach , from the last two of Eqs. (40) we obtain the equations

$$\frac{ds}{dt}=4g\gamma _2w\gamma _1(s\zeta ),\frac{dw}{dt}=2\gamma _2(1gs)w$$
(46)
for the slow variables. These equations can be written in the form (2),

$$\frac{ds}{dt}=v_1,\frac{dw}{dt}=v_2,$$

with the velocity field

$$v_1=4g\gamma _2w\gamma _1(s\zeta ),v_2=2\gamma _2(1gs)w.$$

For the Jacobian matrix (12), we now have

$$J_{11}=\frac{v_1}{s}=\gamma _1,J_{12}=\frac{v_1}{w}=4g\gamma _2,$$

$$J_{21}=\frac{v_2}{s}=2\gamma _2w,J_{22}=\frac{v_2}{w}=2\gamma _2(1gs).$$

From here, we find the contraction rate (17),

$$K=\mathrm{Tr}\widehat{J}=\gamma _12\gamma _2(1gs).$$
(47)
The value of Eq. (47) depends on the characteristics of filaments through the coupling parameter (37) containing in its definition the volume $`V_n=\pi b_n^2L`$, which is the volume of a cylinder enveloping a filament. The radius of the enveloping cylinder, $`b_n`$, is related to the filament radius, $`r_n`$, by relation (34). For what follows, it is convenient to introduce the dimensionless quantity

$$\beta \frac{kb_n^2}{2L}=\frac{\pi b_n^2}{\lambda L}.$$
(48)
Thus, the contraction rate (47) is a function $`K=K(\beta ,t)`$ of the variable (48) and also a function of time entering through $`s=s(t)`$. For different filaments, the variable $`\beta `$ can take different values from the interval

$$0<\beta F,$$
(49)
where $`F`$ is the Fresnel number. The distribution of the variable $`\beta `$ in the interval (49) is given by the probability distribution (20). The maximum of the latter defines the most probable $`\beta `$ and, respectively, the most probable filament radius.

Let us consider the very beginning of the process, when $`t+0`$, in order to understand what are the most probable filaments to be formed. At this initial stage, the probability distribution (20) can be written as

$$p(\beta ,t)\frac{1}{Z(t)}\mathrm{exp}\left\{K(\beta ,0)t\right\}.$$

Hence the maximum of this distribution corresponds to the minimal contraction rate. The latter, according to Eq. (47), depends on $`\beta `$ through $`g=g(\beta )`$. We consider the standard laser-type setup, when the atoms at the initial time are not inverted, i.e. $`s_0<0`$, and the medium is pumped with a long nonresonant pulse whose action is described by the pumping parameter $`\zeta >0`$ in the evolution equations (30), (40), and (46). In this case, the minimum of the contraction rate $`K(\beta ,0)`$ is equivalent to the maximum of the coupling $`g(\beta )`$, which is defined by the equations

$$\frac{dg}{d\beta }=0,\frac{d^2g}{d\beta ^2}<0.$$
(50)
To solve these equations, we need to analyse the integral (37) giving $`g=g(\beta )`$.

The coupling $`g(\beta )`$ defined by the integral (37) can be presented in the form

$$g(\beta )=\frac{3\pi \gamma \rho L}{4\gamma _2k^2}\left[\pi \beta _0^{2\beta }\mathrm{Si}(x)𝑑x\right].$$
(51)
The procedure of reducing Eq. (37) to the form (51) is described in the Appendix. The extrema of $`g(\beta )`$ are given by the equation

$$\mathrm{Si}(2\beta )=\pi .$$
(52)
From several solutions of Eq. (52) we have to choose the absolute maximum of $`g(\beta )`$. The latter occurs at $`\beta =0.96`$, which, according to Eq. (48), gives $`b_n=0.55\sqrt{\lambda L}`$. Using the relation (34), we find the most probable radius of a filament

$$r_f=0.26\sqrt{\lambda L}.$$
(53)
In the typical experiments \[21–23,25–35\] observing the turbulent photon filamentation, the excitation was achieved by means of the quasistationary nonresonant pumping. Such a pumping can be treated as a quasistationary process if its duration is much longer than the characteristic time $`2\pi /\omega `$ of fast oscillations, which is always the case. The quasistationary pumping can be characterized by an effective pumping parameter $`\zeta `$ in the evolution equations. To describe the finite duration time of the pumping procedure, one can consider $`\zeta =\zeta (t)`$ as a slow function of time, slow in the sense of the inequality

$$\frac{2\pi }{\omega }\left|\frac{d\zeta }{dt}\right|1.$$

In this way, for the standard laser setup, the filament radius (53) defines the most probable pattern for the arising bunch of photon filaments.

The most probable number of filaments can be evaluated from the normalization condition

$$\frac{1}{V}s(\stackrel{}{r},t)𝑑\stackrel{}{r}=\zeta (t),$$
(54)
where the integration runs over the whole volume of the sample, $`V=\pi R^2L`$. To approximately calculate the integral (54), let us consider the moment of time when the filaments have already been formed and the population difference inside each filament of the radius $`r_f`$ has reached a value close to the pumping parameter $`\zeta `$. Then we may write

$$s(\stackrel{}{r},t)𝑑\stackrel{}{r}N_fV_f\zeta +(VN_fV_f)s_{out},$$

where $`V_f\pi r_f^2L`$ and $`s_{out}`$ is a population difference outside the filaments. In this case, independently of the value $`s_{out}`$, condition (54) yields

$$N_f=\left(\frac{R}{r_f}\right)^2.$$
(55)
For the filaments of the most probable radius (53), equation (55) gives

$$N_f=4.71F.$$

The number of filaments is proportional to the Fresnel number, which is in complete agreement with all experiments on the turbulent photon filamentation. The coefficient $`4.71`$ is also in good agreement with experiments \[18–20\].

## 6 Dynamics of Filament Flashing

The time evolution of the filaments is described by Eqs. (46). At the initial stage of the process, when $`\gamma _1t1`$, we may omit the term containing $`\gamma _1`$. Then the resulting system of nonlinear equations can be solved exactly yielding the solutions

$$s=\frac{\gamma _0}{g\gamma _2}\mathrm{tanh}\left(\frac{tt_0}{\tau _0}\right)+\frac{1}{g},w=\frac{\gamma _0^2}{4g^2\gamma _2^2}\mathrm{sech}^2\left(\frac{tt_0}{\tau _0}\right),$$
(56)
in which $`\gamma _0`$ is the radiation width, $`\tau _0`$ is the radiation time, so that

$$\gamma _0^2=\mathrm{\Gamma }_0^2+4g^2\gamma _2^2\left(|u_0|^2\alpha _0s_0^2\right),$$
(57)
where $`u_0u(0),s_0s(0),\alpha _0\alpha (0)`$, and

$$\mathrm{\Gamma }_0^2\gamma _2(1gs_0),\gamma _0\tau _01;$$
(58)
the delay time $`t_0`$ being given by the expression

$$t_0=\frac{\tau _0}{2}\mathrm{ln}\left|\frac{\gamma _0\mathrm{\Gamma }_0}{\gamma _0+\mathrm{\Gamma }_0}\right|.$$
(59)
Invoking the transformation (45), we get

$$|u|^2=\frac{\gamma _0^2}{4g^2\gamma _2^2}\mathrm{sech}^2\left(\frac{tt_0}{\tau _0}\right)+\alpha s^2.$$
(60)
The found solutions depend on the coupling $`g`$ both directly as well as through the characteristic parameters $`\gamma _0,\tau _0`$, and $`t_0`$. This means that the properties of filaments with different radii are essentially different, since $`g=g(\beta )`$.

At the later stage, the term with $`\gamma _1`$ cannot be dropped. Then Eqs. (46) cannot be solved exactly. But we may look at the long-time behaviour, when the system approaches a stationary regime. The latter can be achieved if the duration of the pumping pulse is longer than $`\gamma _1^1`$ as well as $`\gamma _2^1`$, which is usually the case when pumping a resonant medium.

Before the stationary regime is reached, the oscillation of solutions is described by the eigenvalues of the Jacobian matrix associated with Eqs. (46). These eigenvalues are

$$\lambda ^\pm =\frac{1}{2}\left[\gamma _1+\gamma _2(1gs)\right]\left\{[\gamma _12\gamma _2(1gs)]^232\gamma _2^2g^2w\right\}^{1/2}.$$
(61)
Because of the dependence of $`s(t)`$ and $`w(t)`$ on time, Eq. (61) shows that $`\lambda ^\pm (t)`$ is also time dependent, which means that the oscillation of solutions cannot be defined by fixed frequencies but occurs in an aperiodic way.

Al long times, such that $`\gamma _1t1`$ and $`\gamma _2t1`$, the system approaches a stationary regime. There are two types of stationary solutions to Eqs. (46), one of them is given by the fixed point

$$s_1^{}=\zeta ,w_1^{}=0;$$
(62)
and another, by the fixed point

$$s_2^{}=\frac{1}{g},w_2^{}=\frac{\gamma _1(g\zeta 1)}{4\gamma _2g^2}.$$
(63)
The stability of the stationary solutions can be determined from the Lyapunov analysis. To this end, we have to evaluate the Jacobian eigenvalues (61) at the corresponding fixed points (62) and (63), which results in the characteristic exponents

$$\lambda _1^+=\gamma _1,\lambda _1^{}=2\gamma _2(1g\zeta )$$
(64)
and, respectively,

$$\lambda _2^\pm =\frac{\gamma _1}{2}\left\{1\pm \left[1+8\frac{\gamma _2}{\gamma _1}(1g\zeta )\right]^{1/2}\right\}.$$
(65)
The real parts of Eqs. (64) and (65) define the Lyapunov exponents whose signs characterize the stability of the related fixed points. As far as the values of the characteristic exponents (64) and (65) depend on the coupling $`g=g(\beta )`$ which, according to Eq. (51), essentially depends on the parameter $`\beta `$, that is, on the radius of a filament, the stability properties for the filaments of different radii can be drastically different. For the parameter $`\beta `$ in the interval (49), the coupling (51) varies in a wide diapason. Thus, for small $`\beta `$, we have

$$g(\beta )\frac{3\gamma \rho \lambda ^2L}{16\gamma _2}\beta (\beta 1),$$

hence $`g(\beta )0`$ as $`\beta 0`$. For large $`\beta `$, Eq. (51) gives

$$g(\beta )\frac{3\gamma L}{16\pi \gamma _2\lambda }\rho \lambda ^3(\beta 1).$$

In the case of the wavelength $`\lambda a`$, much larger than the mean interatomic distance, one has $`\rho \lambda ^31`$. Since usually $`\gamma \gamma _2`$ and $`L\lambda `$, the coupling may reach quite large values $`g(\beta )1`$. Such a wide variation of the coupling $`g=g(\beta )`$ results in rather different characteristic exponents (64) and (65) for different filaments and, as a consequence, in essentially different stability properties of the latter.

All filaments can be separated into three types. One group consists of those filaments for which

$$g(\beta )\zeta <1.$$
(66)
Then the fixed point (62) is a stable node, while that (63) is a saddle point. The filaments with the radii satisfying condition (66) are characterized by the stationary solutions (62), which shows that these filaments after the time $`T_1\gamma _1^1`$ practically stop radiating coherently.

Increasing $`g(\beta )`$, one reaches the equality $`g(\beta )\zeta =1`$, when both fixed points (62) and (63) merge together becoming neutral. In the interval

$$1<g(\beta )\zeta 1+\frac{\gamma _1}{8\gamma _2},$$
(67)
the fixed point (62) is a saddle point and that (63) is a stable node. Hence, the filaments with the radii satisfying condition (67) are described by the stationary solutions (63). Such filaments continue radiating coherently in the stationary regime, although the level of coherence, characterized by the value of $`w_2^{}`$, is not high.

The third group of filaments corresponds to large coupling, such that

$$g(\beta )\zeta >1+\frac{\gamma _1}{8\gamma _2}.$$
(68)
Then the fixed point(62) remains a saddle point while that (63) becomes a stable focus. The latter means that the approach of the related solutions to the stationary point (63) is not monotonic but through a series of pulses. The filaments with the radii satisfying condition (68) radiate by bright flashes interrupted by the intervals of darkness. This intermittent behaviour, for finite times, is not periodic, as is seen from the form of Eq. (61), but as $`t\mathrm{}`$, there appears an asymptotic period. This can be noticed from the characteristic exponent (65) which, for the case (68), can be written as

$$\lambda _2^\pm =\frac{\gamma _1}{2}i\omega _{\mathrm{}},$$
(69)
where the asymptotic frequency is

$$\omega _{\mathrm{}}=\frac{\gamma _1}{2}\left\{8\frac{\gamma _2}{\gamma _1}[g(\beta )\zeta 1]1\right\}^{1/2}.$$
(70)
The value of the latter, because of the dependence on $`g(\beta )`$, is essentially different for different filaments.

## 7 Conclusion

An approach is developed for treating the problem of pattern selection. The approach is based on defining the probability distribution of patterns. This gives the ordering principle for the multiplicity of admissible solutions to nonlinear differential equations and, respectively, to the corresponding space structures. The maximum of the probability distribution naturally defines the most probable pattern.

This probabilistic approach is applied to the problem of turbulent photon filamentation in laser media. The most probable filament radius is found, and the number of filaments is evaluated. The found values are in very good agreement with experiment. Thus, the filament radii observed in experiments with dye and CO2 lasers \[25–30\] are $`r_f0.01`$ cm for dye lasers and $`r_f0.08`$ cm for CO2 lasers, which is in perfect agreement with formula (53). The filament radii measured in experiments with Ne, Tl, Pb, N2, and N$`{}_{}{}^{+}{}_{2}{}^{}`$ vapor lasers \[31–35\] are $`r_f0.01`$ cm, which again coincides with the result of Eq. (53) for the corresponding wavelengths and laser lengths. The number of filaments, given by Eq. (55) is $`N_f10^210^3`$ for different experiments \[31–35\]. The found dependence of the filament radius $`r_fF^{1/2}`$ and filament number $`N_fF`$ on the Fresnel number is the same as observed experimentally.

The turbulent photon filamentation provides a good example of how a rather complicated phenomenon can be successfully described by the probabilistic approach to pattern selection. The formulation of this approach in Section 2 is general, which permits one to apply it for different phenomena where the problem of pattern selection arises.

Acknowledgement

This work is supported by the Bogolubov-Infeld Grant of the State Agency for Atomic Energy, Poland and by a Grant of the São Paulo State Research Foundation, Brazil.

Appendix. Effective Coupling

The effective coupling $`g`$ is given by the integral (37) over the enveloping cylinder of radius $`b_n`$. Assuming that $`\lambda b_n`$, the integral (37) can be reduced to the form

$$g=\frac{3\pi \gamma \rho }{2\gamma _2}_0^{b_n}r𝑑r_{L/2}^{L/2}\frac{\mathrm{sin}(k_0\sqrt{r^2+z^2}kz)}{k_0\sqrt{r^2+z^2}}𝑑z,$$

where $`r`$ is the radial variable. Because of the quasiresonance condition (32), we have $`k_0k`$. Then, with the change of the variable $`xk(\sqrt{r^2+z^2}z)`$, we get

$$g=\frac{3\pi \gamma \rho }{2\gamma _2k}_0^{b_n}r𝑑r_{kr^2/L}^{kL}\frac{\mathrm{sin}x}{x}𝑑x.$$

Here, the integration limit $`kL`$, due to the inequality $`\lambda L`$, can be replaced by $`\mathrm{}`$. As a result, we have

$$g=\frac{3\pi \gamma \rho }{2\gamma _2k}_0^{b_n}\left[\frac{\pi }{2}\mathrm{Si}\left(\frac{kr^2}{L}\right)\right]r𝑑r,$$

where

$$\mathrm{Si}(x)_0^x\frac{\mathrm{sin}t}{t}𝑑t=\frac{\pi }{2}+_{\mathrm{}}^x\frac{\mathrm{sin}t}{t}𝑑t$$

is the integral sine. Introducing notation (48), we come to $`g(\beta )`$ given in Eq. (51). In the same way, for the coupling (38) we find

$$g^{}(\beta )=\frac{3\pi \gamma \rho L}{4\gamma _2k^2}_0^{2\beta }\mathrm{Ci}(x)𝑑x,$$

where

$$\mathrm{Ci}(x)_0^x\frac{\mathrm{cos}t}{t}𝑑t$$

is the integral cosine. The found expressions for the effective couplings can also be transformed by means of the integrals

$$\mathrm{Si}(x)𝑑x=x\mathrm{Si}(x)+\mathrm{cos}x,\mathrm{Ci}(x)𝑑x=x\mathrm{Ci}(x)\mathrm{sin}x.$$

This yields

$$g(\beta )=\frac{3\pi \gamma \rho L}{4\gamma _2k^2}\left[\pi \beta 2\beta \mathrm{Si}(2\beta )+1\mathrm{cos}(2\beta )\right],$$

$$g^{}(\beta )=\frac{3\pi \gamma \rho L}{4\gamma _2k^2}\left[\mathrm{sin}(2\beta )2\beta \mathrm{Ci}(2\beta )\right].$$


format: text

text:
# Light curves of jetted gamma-ray burst afterglows in circumstellar clouds

## 1 Introduction

The origin of gamma-ray bursts (GRBs) has puzzled us for more than three decades (see, e.g., Wijers, Rees & Mészáros 1997; Wijers 1998). In the study of GRBs, two fundamental questions pertinent directly to the central engine need to be solved urgently and unambiguously. One question is about the existence of jets in GRBs. Multi-wavelength observations of afterglows have provided several indications of jetted GRBs: theoretically, it is first shown analytically that the sideways expansion (Rhoads 1999; Sari, Piran & Halpern 1999) and edge effects of jets (Mészáros & Rees 1999) can lead to a marked break in afterglow light curves, and subsequent calculations indicate that one break appears but extends over at least one decade in time (Moderski, Sikora & Bulik 2000; Kumar & Panaitescu 2000; Huang et al. 2000a, b, c; Wei & Lu 2000). These effects seem to account for the observed light curves of the optical afterglows from GRB 990123 (Kulkarni et al. 1999; Castro-Tirado et al. 1999; Fruchter et al. 1999), GRB 990510 (Harrison et al. 1999; Stanek et al. 1999), GRB 990705 (Masetti et al. 2000a), GRB 991208 (Sagar et al. 2000a; Galama et al. 2000a), GRB 991216 (Halpern et al. 2000a), and GRB 000301C (Rhoads & Fruchter 2001; Masetti et al. 2000b; Jensen et al. 2000; Berger et al. 2000; Sagar et al. 2000b). Further evidence for a jet-like geometry in GRBs is the observed radio flares from GRB 990123 (Kulkarni et al. 1999) and GRB 990510 (Harrison et al. 1999). In addition, the polarization observed in some afterglows also provides an important signature for jetted GRBs (Wijers et al. 1999).

Another question is about the association of GRBs with star-forming regions. There is considerable evidence linking the progenitors of GRBs with massive stars. For example, the sources of the GRBs with known redshifts lie within the optical radii and central regions of the host galaxies rather than far outside the disks of the galaxies (Bloom, Kulkarni & Djorgovski 2000), which seems to rule out mergers of neutron-star binaries as the GRB central engine. Further evidence has been provided by the fact that the brightness distribution of GRBs is in agreement with the models in which the GRB rate tracks the star formation rate over the past 15 billion years of cosmic history (Totani 1997; Wijers et al. 1998; Kommers et al. 2000). The most direct evidence for the relation between GRBs and a specific type of supernova (i.e., hypernovae/collapsars) is the discovery of SN 1998bw in the error box of GRB 980425 (Galama et al. 1998) and the detection of a supernova-like component in the afterglows from GRB 980326 (Bloom et al. 1999) and GRB 970228 (Reichart 1999; Galama et al. 2000b). Finally, the recent discovery of a transient absorption edge in the X-ray spectrum of GRB 990705 (Amati et al. 2000) and the observations of X-ray lines from GRB 991216 (Piro et al. 2000) and GRB 000214 (Antonelli et al. 2000) provide new evidence that GRBs are related to the core collapse of massive stars. Based on these observational facts, it is natural to suppose that the environments of GRBs should be pre-burst winds and/or circumstellar clouds. If GRBs occur in such pre-burst winds, their afterglows will fade down rapidly (Dai & Lu 1998; Chevalier & Li 1999, 2000); if GRBs are surrounded by the circumstellar clouds, their fireballs (or jets) must evolve to the non-relativistic regime within a few days after the bursts, leading to a steepening of the afterglow light curves (Dai & Lu 1999, 2000; Wang, Dai & Lu 2000).

The purpose of this paper is to combine these two fundamental points to discuss their implications instead of directly proving their existence. We assume that a GRB comes from a highly collimated jet expanding in a circumstellar cloud, and discuss the resultant afterglow emission when the jet is spreading laterally. In section 2, we analyze the motivation of our work: a strong UV flash from the reverse shock produced by the interaction of the jet with the cloud may clear the dust out to about 30 pc only along the initial path of the GRB. As a result, most of the optical radiation from the high-latitude part of the jet may be strongly extincted by the dust outside the initial solid angle of the jet, but only the radiation from the part within the initial solid angle can be observed. One expects that this could produce a rapidly fading afterglow. In section 3, we further analyze the spectrum and light curve. In section 4, we present our numerical procedure and results. We carry out a detailed modelling for the R-band afterglow of GRB 000926 in section 5, and give a brief discussion and summary in the final section.

## 2 Motivation

In the standard afterglow model (for recent reviews see Piran  and van Paradijs, Kouveliotou & Wijers ), a GRB relativistic shell with a Lorentz factor of $`\eta `$ is assumed to interact with the ambient medium (circumstellar cloud) via two shocks: a reverse shock and a forward shock. The forward shock runs forward into the cloud while the reverse shock sweeps up the shell material. The observed prompt optical flash of GRB 990123 has been argued to come from a reverse shock (Akerlof et al. 1999; Sari & Piran 1999; Mészáros & Rees 1999). We believe that reverse shock emission should be common for all GRBs. The shocked cloud and shell materials are in pressure balance and are separated by a contact discontinuity. Since an optical flash is produced when the reverse shock crosses the shell, the shocked materials expand with a Lorentz factor of $`\gamma `$, which is approximated by the Blandford-MeKee’s self-similar solution. As argued by Waxman & Draine (2000), the reverse shock (with a Lorentz factor of $`\gamma ^{\mathrm{rs}}`$) may be mildly relativistic, i.e., $`|\gamma ^{\mathrm{rs}}1|1`$, and thus $`\gamma \eta `$. We assume that the power-law index of the electron distribution, $`p`$, is similar in the forward and reverse shocks. If the fractions of the thermal energy carried by electrons, $`ϵ_e`$, and magnetic field, $`ϵ_B`$, are also similar in the two shocks, then the typical electron Lorentz factor in the reverse shock $`\gamma _m^{\mathrm{rs}}610\xi ϵ_e`$ and the magnetic field $`B(32\pi \eta ^2nm_pc^2)^{1/2}`$, where $`\xi =3(p2)/(p1)1`$ and $`n`$ is the baryon number density of the cloud.

We consider synchrotron radiation from the reverse shock. To calculate the luminosity of the optical flash, one needs to know the emission spectrum, which is determined by two break frequencies: the peak frequency $`\nu _m^{\mathrm{rs}}`$ and the cooling frequency $`\nu _c^{\mathrm{rs}}`$. We first derive the observed peak synchrotron frequency:

$$\nu _m^{\mathrm{rs}}5.8\times 10^{12}\xi ^2ϵ_{e,1}^2ϵ_{B,6}^{1/2}\eta _{300}^2n_3^{1/2}\left(\frac{1+z}{2}\right)^1\mathrm{Hz},$$
(1)
where $`ϵ_{e,1}=ϵ_e/0.1`$, $`ϵ_{B,6}=ϵ_B/10^6`$, $`\eta _{300}=\eta /300`$, and $`n_3=n/10^3\mathrm{cm}^3`$. The cooling frequency can be approximately expressed as

$$\nu _c^{\mathrm{rs}}10^{19}ϵ_{B,6}^{3/2}\epsilon _{54}^{1/2}n_3^1t_1^{1/2}\left(\frac{1+z}{2}\right)^{1/2}\mathrm{Hz},$$
(2)
where $`t_1`$ is the observer’s time ($`t`$) in units of 10 s (Wang, Dai & Lu 2001). In addition, using equation (4) of Waxman & Draine (2000), we find the observed peak synchrotron flux density for the reverse shock

$$F_{\nu _m}^{\mathrm{rs}}40ϵ_{B,6}^{1/2}\eta _{300}^1n_3^{1/4}\epsilon _{54}^{5/4}t_1^{3/4}\left(\frac{1+z}{2}\right)^{3/4}\left(\frac{\sqrt{1+z}1}{\sqrt{2}1}\right)^2\mathrm{Jy},$$
(3)
where $`\epsilon =\epsilon _{54}\times 10^{54}\mathrm{ergs}\mathrm{sr}^1`$ is the shell energy per unit solid angle, $`z`$ is the redshift of the source, and a flat universe with zero cosmological constant and the Hubble constant $`H_0=65\mathrm{km}\mathrm{s}^1\mathrm{Mpc}`$ is assumed. Therefore, we get the spectrum

$$F_\nu ^{\mathrm{rs}}=\{\begin{array}{cc}(\nu /\nu _m^{\mathrm{rs}})^{(p1)/2}F_{\nu _m}^{\mathrm{rs}}\hfill & \mathrm{if}\nu _m^{\mathrm{rs}}<\nu <\nu _c^{\mathrm{rs}},\hfill \\ (\nu _c^{\mathrm{rs}}/\nu _m^{\mathrm{rs}})^{(p1)/2}(\nu /\nu _c^{\mathrm{rs}})^{p/2}F_{\nu _m}^{\mathrm{rs}}\hfill & \mathrm{if}\nu >\nu _c^{\mathrm{rs}}.\hfill \end{array}$$
(4)
From this spectrum, we can easily calculate the local-frame prompt luminosity ($`L_{17.5}`$) in the $`17.5`$ eV range (UV band). For example, if $`p=2.5`$, we find

$$L_{17.5}1\times 10^{50}ϵ_{e,1}^{3/2}ϵ_{B,6}^{7/8}\eta _{300}^{1/2}n_3^{5/8}\epsilon _{54}^{5/4}t_1^{3/4}\left(\frac{1+z}{2}\right)^{3/4}\mathrm{ergs}\mathrm{s}^1.$$
(5)
Please note that Waxman & Draine (2000) calculated the UV luminosity at $`p=2`$. Here we have extended their discussion to any given $`p`$.

We turn to discuss dust destruction in a circumstellar cloud with a typical radius of $`R_c30`$pc. Since the column density of the cloud is $`N_H=nR_c10^{23}n_3(R_c/30\mathrm{pc})\mathrm{cm}^2`$, the inferred visual extinction is $`50n_3(R_c/30\mathrm{pc})`$ if the cloud is assumed to have the same dust-to-gas ratio as the Galactic clouds, implying that an afterglow from the GRB would be completely extincted without dust destruction. Fortunately, as argued by Waxman & Draine (2000), a strong UV flash from the GRB can destroy its ambient dust by thermal sublimation. Using their equation (17), we find that the dust grains are completely sublimed out to a destruction radius $`R_d30(L_{17.5}/10^{50}\mathrm{ergs}\mathrm{s}^1)^{1/2}`$ pc, where the grain radius $`a=0.1\mu `$m is assumed. Thus, a strong UV flash with luminosity of $`10^{50}\mathrm{ergs}\mathrm{s}^1`$, if isotropic, can clear all the dust grains in a typical giant molecular cloud. Furthermore, Fruchter, Krolik & Rhoads (2001) show that if a GRB emits X-rays in a way similar to those observed by BeppoSAX, dust grains along the line of sight at a distance as large as $`100`$ pc in the host galaxy of the burst can be destroyed by these prompt X-rays due to grain heating and charging. These results were recently confirmed by Galama & Wijers (2001), who analyzed a complete sample of GRB afterglows, and found that the GRB environments have both high column densities ($`10^{23}\mathrm{cm}^2`$) of gas and low optical extinctions. However, since the GRB is believed to come from a jet with the initial half opening angle of $`\theta _00.1`$, the UV flash and prompt X-ray emission will clear the dust only along the initial path of the burst. Therefore, the optical radiation from the $`\theta \theta _0`$ part has low extinction but the optical radiation from the high-latitude ($`\theta >\theta _0`$) part could be absorbed by the dust outside the initial solid angle of the jet, implying that only an afterglow emission from the low-latitude of the jet would be observed if the line of sight is along the jet axis. This could lead to a rapidly fading optical afterglow. Motivated by this argument, we will discuss the afterglow emission from a spreading jet expanding in a circumstellar cloud in the next section.

## 3 Analytical Model: Spectrum and Light Curve

Rhoads (1999) has considered the evolution of a relativistic jet (with Lorentz factor of $`\gamma `$) that is spreading laterally at the local speed of sound $`c_s=c/\sqrt{3}`$ (but $`c_s=c`$ in Sari et al. ), so the half opening angle $`\theta _j\theta _0+\gamma ^1/\sqrt{3}`$. In this case, the dynamical transition takes place at $`\gamma \theta _0^1/\sqrt{3}`$. In the initial stage of the evolution, since $`\gamma >\theta _0^1/\sqrt{3}`$, the jet is spherical-like and its Lorentz factor decays as $`\gamma t^{3/8}`$. The resulting spectrum and light curve are well known (Sari, Piran & Narayan 1998):

$$F_\nu ^{\mathrm{fs}}\{\begin{array}{cc}\nu ^{(p1)/2}t^{3(p1)/4}\hfill & \mathrm{if}\nu _m^{\mathrm{fs}}<\nu <\nu _c^{\mathrm{fs}},\hfill \\ \nu ^{p/2}t^{(3p2)/4}\hfill & \mathrm{if}\nu >\nu _c^{\mathrm{fs}}.\hfill \end{array}$$
(6)
We next stress to discuss the afterglow emission for $`\gamma \theta _0^1/\sqrt{3}`$. In this spreading stage, the sideways expansion leads to an exponential decay of $`\gamma `$ with radius (Rhoads 1999; Sari et al. 1999). As a result, the jet’s radius $`Rt^0`$ and $`\gamma t^{1/2}`$. The typical synchrotron frequency decays as $`\nu _m^{\mathrm{fs}}t^2`$ and the cooling frequency $`\nu _c^{\mathrm{fs}}`$ is a constant. From the discussion in section 2, the total number of electrons radiating toward to the observer, those located in a cone of the initial opening angle, can be estimated as $`N_e=\pi \theta _0^2R^3n/3`$. The total specific luminosity emitted by these electrons, $`(1+z)\sigma _Tm_ec^2N_eB^{}\gamma /(3e)`$, is distributed over an area of $`\pi \gamma ^2D_L^2`$ at the luminosity distance $`D_L`$ from the source (where $`B^{}`$ is the magnetic field strength of the shocked cloud). Thus, the observed peak flux density is given by $`F_{\nu _m}^{\mathrm{fs}}=(1+z)\sigma _Tm_ec^2R^3n\theta _0^2B^{}\gamma ^3/(9eD_L^2)R^3\gamma ^4t^2`$. Therefore, we obtain the afterglow’s spectrum and light curve:

$$F_\nu ^{\mathrm{fs}}=\{\begin{array}{ccc}(\nu /\nu _m^{\mathrm{fs}})^{1/3}F_{\nu _m}^{\mathrm{fs}}\nu ^{1/3}t^{4/3}\hfill & \mathrm{if}\nu <\nu _m^{\mathrm{fs}},\hfill & \\ (\nu /\nu _m^{\mathrm{fs}})^{(p1)/2}F_{\nu _m}^{\mathrm{fs}}\nu ^{(p1)/2}t^{(p+1)}\hfill & \mathrm{if}\nu _m^{\mathrm{fs}}<\nu <\nu _c^{\mathrm{fs}},\hfill & \\ (\nu _c^{\mathrm{fs}}/\nu _m^{\mathrm{fs}})^{(p1)/2}(\nu /\nu _c^{\mathrm{fs}})^{p/2}F_{\nu _m}^{\mathrm{fs}}\nu ^{p/2}t^{(p+1)}\hfill & \mathrm{if}\nu >\nu _c^{\mathrm{fs}}.\hfill & \end{array}$$
(7)
It is easy to see that the temporal decay index of the high-frequency afterglow changes from $`\alpha _1=3(p1)/4`$ or $`(3p2)/4`$ to $`\alpha _2=p+1`$ because of the effects of sideways expansion and dust extinction. In conclusion, a marked break should appear in the light curve of the afterglow from a spreading jet expanding in the circumstellar giant cloud.

## 4 Numerical Results

To prove the analytical result in section 3, we perform detailed numerical calculations for the evolution of the afterglow emission. We use the model proposed by Huang, Dai & Lu (1999) and the calculational code developed by Huang et al. (2000c). This model has several advantages: (i) It is applicable to both radiative and adiabatic jets, and proper for both ultra-relativistic and non-relativistic stages. The model even allows the radiative efficiency $`ϵ`$ to change with time, so that it can trace the evolution of a partially radiative jet (Dai, Huang & Lu 1999). However, note that some authors (e.g., Mészáros, Rees & Wijers 1998) have argued that the jet should be adiabatic most likely, since it is unlikely that the radiative efficiency could reach 1; (ii) The model considers the lateral expansion of the jet. The evolution of the lateral speed (taken as the speed of sound) is given by a reasonable expression (see equation ). (iii) The model also takes into account many other effects, e.g., the cooling of electrons and the equal arrival time surfaces. Let $`M_{\mathrm{ej}}`$ be the initial ejecta mass and $`m`$ the swept-up cloud mass. The dynamical evolution of the jet is described by (Huang et al. 2000c)

$$\frac{dR}{dt}=\beta c\gamma (\gamma +\sqrt{\gamma ^21}),$$
(8)
$$\frac{d\theta _j}{dt}=\frac{c_s(\gamma +\sqrt{\gamma ^21})}{R},$$
(9)
$$\frac{dm}{dR}=2\pi R^2(1\mathrm{cos}\theta _j)nm_\mathrm{p},$$
(10)
$$\frac{d\gamma }{dm}=\frac{\gamma ^21}{M_{\mathrm{ej}}+ϵm+2(1ϵ)\gamma m},$$
(11)
$$c_s^2=\frac{\widehat{\gamma }(\widehat{\gamma }1)(\gamma 1)c^2}{1+\widehat{\gamma }(\gamma 1)},$$
(12)
where $`\beta =\sqrt{11/\gamma ^2}`$ and $`\widehat{\gamma }(4\gamma +1)/(3\gamma )`$ (Dai et al. 1999). If $`ϵ0`$ as in the following calculations, equation (11) turns out to express the conservation of energy: $`(\gamma 1)M_{\mathrm{ej}}+(\gamma ^21)m=\mathrm{constant}`$ (for a discussion see Huang et al.  and van Paradijs et al. ).

Figure 1 presents light curves of R-band afterglow emission ($`p=2.5`$) for two cases with dust extinction (solid line) and without dust extinction (dashed line). This figure clearly shows that the light curve with dust extinction is indeed more steepening than that without dust extinction at late times. Figure 2 exhibits $`\alpha d\mathrm{ln}F_R/d\mathrm{ln}t`$ as a function of time for these two cases. For the case without dust extinction, $`\alpha `$ increases from $`1.3`$ at initial one day to $`2.5`$ at later times. This change extends over one order of magnitude in time, which is consistent with the previous numerical results (Moderski et al. 2000; Kumar & Panaitescu 2000; Huang et al. 2000a, b, c). However, for the case with dust extinction, $`\alpha `$ increases from $`1.3`$ at initial one day to $`3.3`$ at later times, which is in approximate accord with the analytical result in section 3. Furthermore, the steepening is completed over a factor of $`3`$ in time, leading therefore to a sharper break in the light curve.

## 5 Afterglow of GRB 000926

GRB 000926 was detected on 2000 September 26.9927 UT by the Inter-Planetary-Network (IPN) group of spacecrafts Ulysses, Russian Gamma-Ray Burst Experimant (KONUS) and Near Earth Asteroid Rendezvous (NEAR) (Hurley et al. 2000). The burst lasted $`25`$ s and had a $`25100`$ keV fluence of $`2.2\times 10^5\mathrm{ergs}\mathrm{cm}^2`$. The redshift of the source was determined at $`z=2.0369\pm 0.0007`$ (Castro et al. 2000), yielding a luminosity distance of $`D_L=16.9`$ Gpc. The optical light curve of the afterglow from GRB 000926 fell off as $`t^{1.1\pm 0.2}`$ for the first $`t_{\mathrm{br}}=2.0\pm 0.4`$ days (break time) and subsequently steepened to $`t^{3.2\pm 0.4}`$ (Rol, Vreeswijk & Tanvir 2000). Furthermore, the resultant sharp break extended over a factor of $`3`$ in time, inferred from the fitted function of Rol et al. (2000). In addition, from the optical to X-ray observations, Sagar et al. (2000c) have found that the spectral index $`\overline{\beta }=0.82\pm 0.02`$ at $`t=2.26`$. Recently, Price et al. (2001) also fitted a late-time light curve for this afterglow, which slightly flattens as compared with Rol et al.’s (2000).

In our spreading jet model (with dust extinction), the spectral index $`\overline{\beta }=0.82\pm 0.02`$ requires $`p=2.6\pm 0.08`$, implying $`\alpha _1=1.2\pm 0.06`$ (in the slow cooling regime) prior to the break time $`t_{\mathrm{br}}`$ and subsequently $`\alpha _2=3.6\pm 0.08`$ (from equation ). These values are in excellent agreement with the observations. Figure 3 shows the observed R-band data for GRB 000926 and the theoretically calculated light curve based on the model described in section 4. It can be seen from this figure that our model provides a good fit to the optical afterglow data. More importantly, our model accounts successfully for the observed sharpness of the break that appears in the R-band light curve. We note that the previous spreading jet model (without dust extinction) seems inconsistent with the observed data. This is both because in the previous jet model the temporal decay index increases from $`1.3`$ to $`2.6`$ if $`p=2.6`$ and because the sharpness of the observed break cannot be produced in that model, as shown in the present paper and in other studies (e.g., Moderski et al. 2000; Kumar & Panaitescu 2000; Huang et al. 2000a, b, c).

## 6 Discussion and Conclusions

We have discussed the afterglow emission when a highly collimated jet is spreading laterally and expanding radially in a circumstellar cloud, and have analytically found that the flux decay of the observational optical afterglow changes from $`t^{3(p1)/4}`$ (or $`t^{(3p2)/4}`$) in the relativistic spherical-like phase to $`t^{(p+1)}`$ in the relativistic sideways-expansion phase because of strong dust extinction outside the initial solid angle of the jet. This provides a natural explanation for very rapidly fading afterglows. We have also performed numerical calculations based on the generic afterglow model proposed by Huang et al. (1999) and developed by Huang et al. (2000c). Our numerical results show that one break in the light curve of the optical afterglow extends over a factor of $`3`$ in time in our present model rather than one decade in time in the previous jet model. These results may provide a signature for GRBs in dense clouds. In addition, we have given a good fit to the observed R-band data of the afterglow from GRB 000926.

It is interesting to discuss two effects of dust on afterglows. First, thermal reradiation and scattering outside the destruction radius in the circumstellar giant cloud are expected to cause delayed IR emission (Waxman & Draine 2000; Esin & Blandford 2000), which has been proposed to account for the late-time afterglow feature of GRB 980326 and GRB 970228 instead of the supernova-like component explanation. At $`z=2`$ (e.g., GRB 000926), the time delay and the flux for the emission at $`2.2\mu `$m can be estimated by $`t_{\mathrm{IR}}1.5\times 10^7(\theta _0/0.1)^2(L_{17.5}/10^{50}\mathrm{ergs}\mathrm{s}^1)^{1/2}`$ s, and $`F_{2.2\mu \mathrm{m}}0.2(L_{17.5}/10^{50}\mathrm{ergs}\mathrm{s}^1)^{1/2}\mu `$Jy, respectively. This IR flux is expected to be detected in about six months after the burst. Second, since the visual extinction is $`50`$, the X-rays could be both scattered by the dust grains (e.g., Mészáros & Gruzinov ), whose scattering optical depth at the X-ray energy $`ϵ_\gamma `$ is $`\tau 15(ϵ_\gamma /1\mathrm{keV})^2`$, and efficiently heat the grains. Also, the soft X-rays could be absorbed in the giant dense cloud because of a high optical depth for photo-ionization. These effects may lead to the rapidly fading X-ray afterglow of GRB 000926.

## Acknowledgments

We wish to thank the referee, Dr. Ralph Wijers, for valuable comments, and X. Y. Wang for useful discussions. This work was supported partially by the National Natural Science Foundation of China (grants 19825109, 10003001 and 19773007), by the National 973 Project (NKBRSF G19990754), and by the Special Funds for Major State Basic Research Projects.


text:
# Raman spectrum and lattice parameters of MgB2 as a function of pressure

## Abstract

We report Raman spectra and synchrotron x-ray diffraction measurements of lattice parameters of polycrystalline MgB2 under hydrostatic pressure conditions up to 15 GPa. An anomalously broadened Raman band at 620 cm-1 is observed that exhibits a large linear pressure shift of its frequency. The large mode damping and Grüneisen parameter indicate a highly anharmonic nature of the mode, broadly consistent with theoretical predictions for the E2g in-plane boron stretching mode. The results obtained may provide additional constraints on the electron-phonon coupling in the system.

The newly discovered  high-temperature superconductor MgB2 has attracted considerable interest from theoretical and experimental points of view. Theory indicates that MgB2 can be treated as phonon mediated superconductor with very strong coupling . Calculations show that the strongest coupling is realized for the near-zone center in-plane optical phonon (E2g symmetry) related to vibrations of the B atoms . According to recent calculations, this phonon is very anharmonic because of its strong coupling to the partially occupied planar B $`\sigma `$ bands near the Fermi surface . The frequency of this phonon ranges from 460 to 660 cm-1 according to different computation techniques . The phonon density of states for MgB2 has been determined by neutron inelastic scattering , but the E2g mode could not be detected separately. Raman experiment  indicated a presence of a broad mode at 72 meV (580 cm-1) in agreement with calculations for the E2g mode. Transport, magnetic susceptibility, and specific heat measurements show a large isotope effect consistent with phonon mediated superconductivity .

Pressure is an important variable, that can be used to tune physical properties and compare the results with theoretical predictions. Pressure effects on superconductivity studied to 1.84 GPa  and 0.5 GPa  show a decrease of $`T_c`$ with the rate of 1.6 K/GPa and 1.11 K/GPa, respectively (see also Ref. ). Compressibility data have been obtained by neutron diffraction (to 0.62 GPa)  and synchrotron x-ray diffraction (to 6.15 GPa  and 8 GPa ). Based on theoretical calculations of the electronic density of states at the Fermi level, which show a very moderate decrease with pressure, the dominant contribution to the decrease of $`T_c`$ under pressure has been proposed to be due to an increase in phonon frequency .

In this Letter we present Raman measurements of the phonon mode including the first measured under pressure. We find that the E2g band is unusually broaden and shows a large positive pressure shift of the frequency. We also present x-ray diffraction data, which allow us to determine lattice parameters on the sample from the same batch in purely hydrostatic conditions to 12 GPa. As a result, we determined the E2g mode Grüneisen parameter, which is much larger than in case of ”normal” materials. We ascribe this to large anharmonic effects predicted by theory. The increase in phonon frequency measured in this work explains the reported TC drop with pressure.

Samples of Mg10B2 were similar to those used in Refs. . They are essentially in a powdered form consisting of aggregates of 30-50 $`\mu `$m linear dimensions, which is ideal for high-pressure experiments. Our experiments have been done with various types of diamond anvil cells. In Raman experiments a long piston-cylinder cell was used and Ne served as a pressure transmitting medium . Synthetic ultrapure diamonds were used as anvils to reduce background fluorescence. Raman scattering was excited in a 145 geometry (see Ref. ) to reduce further background from diamond Raman and that originating from spuriously reflected elastic light. The spectra were recorded with a single-stage spectrograph equipped with a CCD detector and holographic notch filters (150-5000 cm-1), although occasional measurements were also done with a conventional triple spectrometer to cover the lower frequency range. X-ray diffraction was measured with a two-side open diamond cell in an energy-dispersive configuration at X17C beamline of the National Synchrotron Light Source with 2$`\theta `$=10 . In the x-ray experiment we used helium as a pressure transmitting medium, which is purely hydrostatic to 12 GPa. Pressure was determined by the standard ruby fluorescence technique. All measurements were performed at room temperature.

Figure 1 presents the Raman spectra at different pressures. The broad band observed is a Raman excitation as shown by changing the excitation wavelength and by anti-Stokes measurements. It has also been checked that the signal originates from MgB2 because identical spectra were recorded by separate micro Raman measurements from individual micron-size grains (shown in Fig. 1 as the 0 GPa spectrum). Also, the Raman spectra contain a wide unstructured background component (presumably of electronic origin as in the cuprate HTSC materials ), which increases intensity at lower frequencies. Pressure leads to an increase in the frequency of the broad band without any appreciable change of its shape. The spectra can be fitted reasonably well with a combination of a linear background and a Gaussian peak. The frequency determined by this procedure is plotted as a function of pressure in Fig. 2. The pressure dependence is linear within the accuracy of the experiment. No essential pressure dependence of the mode damping was found (inset in Fig. 2).

Factor-group analysis predicts for MgB2 (space group $`P6/mmm`$, Z=1) B1g +E2g\+ A2u+E2u zone center optical modes, of which only E2g is Raman active. Thus, it is natural to assign the band observed at 620 cm-1 at ambient conditions to the E2g mode (see also Ref. ). The experimental frequency agrees well with theoretical calculations . The anomalously large linewidth (FWHM=300 cm-1) can be ascribed to large electron-phonon coupling , which will be described below.

The experimental pressure dependencies of lattice parameters determined by x-ray diffraction are shown in Fig. 3. Our data are in good agreement with Refs. , while the results of Ref.  show systematically larger lattice parameters and yet comparable compressibility. We calculated the bulk modulus K0 assuming ”normal” behavior and K$`{}_{}{}^{}{}_{0}{}^{}`$=4, which is typical for covalent and metallic bonding (our data do not allow us to fit data with two parameters K0 and K$`{}_{}{}^{}{}_{0}{}^{}`$). The result is 155(10) GPa in good agreement with Refs.. Similar calculations for in-plane and out of plane compressibilities give $`\beta _a`$=0.0016(2) GPa-1 and $`\beta _c`$=0.0030(2) GPa-1.

Thus, the mode Gr̈uneisen parameter $`\gamma `$=K0dln$`\nu `$/dP determined from our data equals 2.9$`\pm `$0.3. In the case of anisotropic crystals it would be more appropriate to scale the frequency shift of in-plane mode with the variation of interatomic bond distance or lattice parameter $`a`$ . The corresponding component of the Grüneisen parameter ($`\gamma `$=dln$`\nu `$/3dln$`a`$) is 3.9$`\pm `$0.4. These values are substantially larger than those expected for the phonon in a compound with covalent bonding , which should be dominant for this mode, where typically, $`\gamma `$ does not exceed 2. For example, for graphite $`\gamma `$=1.06  and for iron (with metallic bonding partially present in our case) $`\gamma `$=1.7 . Larger $`\gamma `$’s are normally related to increased anharmonicity of the particular normal vibration . It can also be a consequence of a soft mode behavior when the system is approaching (or departing from) a structural instability (e.g. Ref. ).

The proposed assignment of the Raman peak observed to the first-order phonon scattering is not the only possible. Alternatively, if the first order scattering is inherently weak, the observed Raman peak can in principle be second order (e.g., due to overtones and combinations of the zone boundary acoustical phonons). However, this interpretation does not seem plausible because no higher frequency peak corresponding to combinations of acoustic and optical or two optical modes is observed. Also, the observed excitation may not necessarily be of phonon nature, but in principle could be a magnetic excitation  (e.g. two-magnon peak, which is strongly dependent on interatomic distances with $`\gamma `$=3.5 ). However, we believe that the data available strongly suggest the first-order phonon interpretation because of the agreement with the calculated frequency  and the linewidth .

Theoretical calculations  suggest a scenario with the E2g phonon strongly coupled to electronic excitations. Our data show a very broad, strongly pressure dependent excitation, which is consistent with this idea. Conventional anharmonicity (not coupled to electronic degrees of freedom) is expected to exhibit some variation with pressure (see e.g. Ref. ), which is not the case. Within this picture, our data favor the coupling of the E2g phonon to the electronic subsystem.

Finally, we address the observed strong pressure dependence of $`T_c`$ . Assuming a pressure independent density of electronic states N(0), the averaged electron-ion matrix element $`I`$ and the Coulomb pseudopotential $`\mu ^{}`$, one can get d$`T_c`$/dP=-1.5 K/GPa with physically reasonable values of $`\mu ^{}`$=0.005-0.1 and electron-phonon coupling constant $`\lambda `$=0.65-1 . Thus, the pressure dependence of $`T_c`$ can be easily explained by an increase in phonon frequency as proposed in Refs. .

In conclusion, we observed a strongly broadened Raman band of MgB2 that shows anomalously large pressure dependence of its frequency. This band and its pressure dependence can be interpreted as the E2g zone center phonon, which is strongly anharmonic because of coupling to electronic excitations.

We acknowledge financial support of CIW, NSLS, NSF and Keck Foundation. Work at Ames laboratory was supported by the Director for Energy Research, Office of Basic Energy Sciences, U.S. Department of Energy.


text:
# Untitled Document

Quantum Cable as transport spectroscopy of 1D DOS of cylindrical quantum wires

Y. Xiang, Z. Y. Zeng, and L. D. Zhang

Institute of Solid State Physics, Chinese Academy of Sciences, P.O. Box 1129, Hefei,

230031, People’s Republic of China

Abstract

We considered the proposed Quantum Cable as a kind of transport spectroscopy of one-dimensional (1D) density of states (DOS) of cylindrical quantum wires. By simultaneously detecting the direct current through the cylindrical quantum wire and the leaked tunneling current into the neighboring wire at desired temperatures, one can obtain detailed information about 1D DOS and subband structure of cylindrical quantum wires.

PACS numbers: 73.23.-b, 73.23.Ad, 71.20. -b, 73.50.-h

* Correspondence to whom should be addressed.

Tel: 0551-5591444 (office), 5591453 (home)

Fax: 0551-5591434

E-mail:zyzeng@mail.issp.ac.cn

Nanometer structure and materials are of great interest as potential candidates for future electronic devices of greatly reduced size1. The fundamental physics underlying these devices is quantum coherence and quantum confinement2. Rational design and fabrication will require a better understanding of how the properties of these materials depends on, for example, their dimensionality and size. With the rapid development of nanofabrication and synthesis technology, a variety of low-dimensional nanostructures with unusual characteristics have been realized, such as quantum well, quantum wire, quantum dot as well as their composite structures3.

As a kind of quantum wire structure, nanotubes4 have received extensive attention both theoretically and experimentally. Most recently, the group of Iijima5 has successfully synthesized a new kind of coupled nanotube structure - coaxial nanocable, in which two conducting cylindrical layers are coupled through a very thin insulating layer. Motivated by this new nanostructure, we proposed Quantum Cable6 similar to coaxial nanocable structure, and have studied its energy subband spectrum and ballistic transport properties. Quantum Cable comprises two cylindrical quantum wires and a tunable potential barrier, which is sketched in Fig. 1 (a). It is achievable either by confining electrons in the nanocable-like potential wells, or by curling multilayer heterostructure into a solid cylinder or hollow cylinder with nanometer width. Since Quantum Cable is a coupled-quantum-wire structure, it is very akin to the usual coupled dual-quantum-well waveguides7 (DQWs) in two-dimensional electron gas (2DEG). On the other hand, Quantum Cable is different from the DQWs due to its cylindrical symmetry6. For this reason Quantum Cable is a unique nanostructure and expected to be used as some kind of quantum-effect device. In this letter, we considered Quantum Cable as a transport spectroscopy for detecting quasi-one-dimensional density of states (1D DOS) (bound states) in cylindrical quantum wires. The system will be studied is two GaAs cylinder wires surrounded by $`Ga_{0.7}Al_{0.3}As`$ for which the parameters are $`U_B=100`$ mev, $`m_1^{}=5.73\times 10^{32}kg`$ and $`m_2^{}=1.4m_1^{}`$. The widths of two cylindrical wires are both $`50`$ nm and the barrier width is set to be $`5`$ nm.

The bias configuration \[see Fig. 1 (b)\] is as follows: a small constant bias $`V_D`$ ( 100 $`\mu `$v) drives electrons’ direct flowing through one of the cylindrical quantum wires into the drain, while another variable bais voltage $`V`$ is applied to the side-wall of the same quantum wire to adjust the electrons’ Fermi potential. In the forthcoming calculations, we directly relate the applied side-wall voltage to the Fermi potential of electrons in quantum cylinder for simplicity. In fact, this needs self-consistent calculation. Note that such simplication will not lead to any confusion. The potential barrier is set to be good isolation for qauntum wires, yet allow for electron’s tunneling between two wires. Then tunneling current through the barrier can be measured simultaneously in another wire. From Landauer-Büttiker formula8, the direct current $`I_D`$ along the cylindrical wire from source to drain can be derived as

$$I_D=\frac{2e^2}{h}\underset{mn}{}_{E_{mn}}^{\mathrm{}}𝑑E[\frac{f(EE_FeV)}{E}]V_D,$$
(1)
and the tunneling current $`I_T`$ from one cylinder wire to another is given by

$$I_T=e\underset{mn}{}_{\mathrm{}}^+\mathrm{}𝑑Ev_{}(EE_{mn})g_{mn}^{1D}(EE_{mn})T_{mn}(EU_B)[f(EE_FeV)f(EE_F)],$$
(2)
where we have accounted for the contribution to the current from every occupied subband $`(mn)`$. $`E_F`$ is electrons’ Fermi energy as $`V=0`$, and will be set to be zero for convinience in our actual numeration, $`E_{mn}`$ is the energy to be determined of the $`(mn)th`$ subband and $`U_B`$ is the height of the tunneling barrier. $`f(E)=[1+e^{E/(k_bT)}]^1`$ is the Fermi distribution function with $`k_b`$ the Boltzmann constant. $`v_{}(EE_{mn})=(2E_{mn})^{1/2}/(2m_1^{})`$ is the velocity of tunneling electrons at subband (mn). $`g_{mn}^{1D}(EE_{mn})=\theta (EE_{mn})(EE_{mn})^{1/2}/(\pi \mathrm{})`$ is the 1D DOS. The transmission coefficient is approximated to be that through a 1D square barrier, i.e., $`T(EU_B)=[1+sh^2(\sqrt{2m_2^{}(U_BE)}R_B/\mathrm{})/(E/U_BE^2/U_B^2)]^1`$ $`(E<U_B)`$ . The subband energy $`E_{mn}`$ is followed from

$$J_n(\sqrt{2m_1^{}E}R_1/\mathrm{})Y_n(\sqrt{2m_1^{}E}R_2/\mathrm{})Y_n(\sqrt{2m_1^{}E}R_1/\mathrm{})J_n(\sqrt{2m_1^{}E}R_2/\mathrm{})=0,$$
(3)
where $`J_n(x),Y_n(x)`$ are the first-order,second-order Bessel functions, $`R_1`$, $`R_2`$ are the inner and outer radius of the hollow cylindrical quantum wire considered, respectively.

Once the subband energy is obtained, the 1D DOS of the cylindrical wire can be calculated. In Fig. 2 we present the calculated 1D DOS for a hollow cylindrical quantum wire of width $`50`$ nm , in the approximation of infinite potential well. It is worth noting that subband spectrums are the same in the cases of both infinite and finite confining potential well model, except that subband energy in the case of infinite potential well is slightly higher than that in the finite potential well case for not very small wire width9. From Fig. 2, one can find that, over the consecutive energy range, subband structure and 1D DOS of a hollow quantum cylinder are very similar to that of a single 2D quantum waveguide (SQW). Whereas irregular subband arrangement appears in the intersection between two energy regions with SQW-like subband structure. This feature arises from the mixing of subbands of lowest-order azimuthal quantum number $`m`$ and higher-order radial quantum number $`n`$ with that of large $`m`$ and small $`n`$. For solid cylinder, its energy subbands are organized somewhat irregularly6.

In Fig. 3, we give the direct current $`I_D`$ through a quantum cylinder and tunneling current $`I_T`$ leaked into another cylinder at absolute zero temperature, as the side-wall voltage bias is adjusted. It can be shown that the direct and tunneling currents decrease as $`V`$ is made smaller, until at $`V<2.2`$ mv, all current-carrying channels (subbands) are turned off and no direct and tunneling current is observed. If one examines carefully the detailed structure of direct and tunneling currents, one can see the perfect quantized steps in the direct current $`I_D`$, and the 1D DOS-similar peak oscillation in the tunneling current $`I_T`$. Moreover, the start of each step in $`I_D`$ lines perfectly up with the corresponding peak in $`I_T`$, as indicated by the arrows in Fig. 3. By counting the number of steps in $`I_D`$ and peaks in $`I_T`$, we find they are the same. In addition, two steps of conductance magnitude $`2e^2/h`$ can be discerned from the other steps of conductance unit $`4e^2/h`$. The former steps come from the contribution of the conducting subbands belonging to the azimuthal quantum number $`m=0`$, while the later result from the conducting subbands $`m0`$ which are doubly degenerated. This reflects the cylindrical symmetry of cylindrical quantum wires. As the side-wall voltage increases, electron’s Fermi potential sweep the transverse energy subband one by one and thus conducting channels, each of which carrying the same ammount of current, opens up progressively. Then we can observe quantum steps in $`I_DV`$ characteristic. As to the peak structure in the tunneling current $`I_T`$, it is the direct consequence of tunneling current being propotional to the 1D DOS for a given subband, as can be found in Eq. (2).

At nonzero temperature, it is expected that the plateaus in the direct current will be rounded and the peak in tunneling current will be broadened, until these steps and peaks are smeared out by thermal effect. Fig. 4 shows the currents $`I_D`$ and $`I_T`$ vs side-wall bias at temperatures $`T=0.2,1,3.5`$ K. It is evident that, with the increase of the temperature, the plateaus in $`I_D`$ are more rounded and the peaks in $`I_T`$ more broadened. As the thermal energy is comparative to the energy interval between two adjoining subbands, the original two distinct steps at low temperature will be blurred, and two narrow peaks will be mingled into one broadened peak. At around $`T=3.5`$ K, the step-fasion structure in $`I_D`$ as well as the peak oscillations of $`I_T`$ is washed out. This is consistent with the predicted smearing-out temperature by Bagwell and Orlando10. If the width of cylindrical quantum wires is decreased, intersubband energy spacing will increase, which implies a higher operatation temperature, though the quantum steps in direct current still will be smeared out at about $`T=3.5`$ K becuase of fixed height of current step. As a result, it is more convinient and more accuate to determine the 1D DOS and subband structure of cylindrical quantum wire from the leaked tunneling current than from the direct current into the drain. However, one can not gain detailed insight into the 1D DOS and subband structure of cylindrical quantum wire simply from the data of tunneling current. Since at some temperature, two narrow peaks in tunneling current corresponding to the neighboring subbands of small energy interval will be smeared into one broadened peak, at this moment one can notice this fact according to the slope of the direct current at corresponding side-wall bias, and make good resolution of such neighboring subbands . Therefore it is necessry to measure simultaneous the direct current in the drain and leaked current into another cylindrical quantum wire.

Based on the aforementioned observations, one can determine the energy subband structure and 1D DOS of cylindrical quantum wire from the simultaneous measurement of the direct current through the wire and the leaked current into the neighboring cylinder wire at accessible temperatures. To do so, one should first form Quantum Cable structure from a single cylindrical quantum wire using either fabrication or synthesis approach. Then apply bias voltages with the configuration mentioned above to the obtained Qauntum Cable. By simultaneous measuring the direct current and tunneling current at desired temperature, one can capture detailed information about the subband structure and 1D DOS of a cylindrical quantum wire. Experimentally, such transport spectroscopy is feasible, since two-terminal11 and multi-terminal12 measurement of probing electriccal transport in individual single-wall carbon nanotubes have been achieved.

This work is supported by a key project for fundamental research in the National Climbing Program of China, one of us (ZYZENG) would like to thank Prof. G. H. Li, Prof. G. W. Meng, Dr. Z. Cui for helpful discussions.

References

1. T. W. Ebbesen, Annu. Rev. Mater. Sci. 24, 235 (1994); C. -G. Wu and T. Bein, Science 266, 1013 (1994).
2. B. L. Altshuler, P. A. Lee, and R. A. Webb, Eds., Mesoscopic Phenomenon in Solids (North-Holland, Amsterdam, 1991).
3. M. A. Reed, ed., Nanostructured Systems (Academic Press, New York, 1992).
4. S. Iijima, Nature 354, 54 (1991); M. S. Dresseelhaus, G. Dresseelhaus, and P. C. Eklund, Science of Fullerenes and Carbon Nanotube (Academic, San Diego, 1996) and references cited.
5. X. Suenaga, C. Colliex, N. Demoncy, A. Loiseau, H. Pascard, and F. Willaime, Science 278, 653 (1997); Y. Zhang, K. Suenaga, C. Colliex and S. Iijima, Science 281, 973 (1998).
6. Z. Y. Zeng, Y. Xiang, and L. D. Zhang (unpublished).
7. J. A. del Alamo and C. C. Eugster, Appl. Phys. Lett. 56, 78 (1990); Phys. Rev. Lett. 67, 3586 (1991).
8. R. Landauer, Philos. Mag. 21, 863 (1970); M. Büttiker, Phys. Rev. Lett. 57, 1761 (1986).
9. N. C. Constantinou and B. K. Ridley, J. Phys.:Condens. Matter 1,2283 (1989).
10. P. F. Bagwell and T. P. Orlando, Phys. REv. B 40, 1456 (1989).
11. H. Dai, E. W. Wong, and C. M. Lieber, Science 272, 523 (1996); S. J. Tans, M. H. Devoret, H. Dai, A. Thess, R. E. Smalley, L. J. Geerliges, and C. Dekker, Nature 386, 474 (1997).
12. A. Bezryadin, A. R. M. Verschueren, S. J. Tans, and C. Dekker, Phys. Rev. Lett. 80, 4036 (1998).

Figure Captions

Fig. 1 (a) Schematic view of Quantum Cable structure, (b) bias configuration to the probing of direct and tunneling currents in Quantum Cable.

Fig. 2 1D DOS of hollow quantum cylinders with width $`50`$ nm in the approximation of infinite potential well.

Fig. 3 Variation of the direct current $`I_D`$ through the quantum cylinder and tunneling current $`I_T`$ as the side-wall voltage bias is altered, at absolute zero temperature.

Fig. 4 Variation of the direct current $`I_D`$ through the quantum cylinder and tunneling current $`I_T`$ as the side-wall voltage bias is altered, at nonzero temperatures.