text:
# On the meteor height from forward scatter radio observations

## 1 Introduction

A meteoroid enters the Earth’s atmosphere at hypersonic speed and it collides with air molecules. The high kinetic energy involved in the process determine the transformation of a solid body into a plasma, which can scatter radio waves and can emit light (meteor).

During sixties and seventies several works investigated the formation and evolution of the meteor, with particular attention to diffusion, in order to study mesospheric winds. A complete review of standard meteor science can be found in Ceplecha et al. (CEPLECHA1 (1998)). However, there are still some aspects not well understood about the physical properties of a meteor, specifically whether it is an ionized gas or a plasma. During past years, these two appellatives were often used as synonymous in meteor physics, even though they indicate two different states of the matter. In some studies, such as those about diffusion, specific plasma properties are taken into account (e.g. ambipolar diffusion); however in other studies, such as about radiowave scattering, the meteor is simply considered a long narrow column of ionized gas.

This can appear as a futile debate, but it hides important concepts. Specifically, a plasma has collective properties (e.g. Langmuir frequency) that an ionized gas has not.

A first attempt to study the meteor as a plasma was carried out by Herlofson (HERLOFSON (1951)). He investigated the proper oscillations in the meteor and their interaction with radio waves. But, at our knowledge, none continued his studies. Only in 1999 the question of collective oscillations in meteoric plasma was reprised (Foschini FOSCHINI1 (1999)). Perhaps, this gap may be explained by taking into account that, according to purposes of meteor astronomy, it was sufficient to use the approximation of the long narrow cylinder.

However, the meteoric plasma is something more complex than a reflecting rod and it is necessary to study it. There are several types of oscillations and instabilities, which can interact with radio waves. The scattering is not the only process: for example, fluctuations from equilibrium may lead to transformation of waves (longitudinal to transverse and vice versa). The question is: are such processes present in a meteoric plasma?

We think that the study of plasma collective oscillations may give new useful tools to understand the physics of meteors. Some basic concepts about meteoric plasma were settled in a previous paper (Foschini 1999), thereafter called Paper I. According to the theory exposed there, radio echoes can be divided into two classes and two subclasses. Then we have underdense and overdense echoes, according to whether the Langmuir frequency is higher or lower than the radio wave frequency. Overdense echoes totally reflect electromagnetic waves, but the presence of binary collisions among ions and electrons weaken the collective oscillations of the plasma, allowing the propagation of the waves, even though with strong attenuation. Therefore, we can divide the overdense echoes into two subclasses: type I, when there is total reflection; type II, when binary collisions allow the propagation. The division between overdense type I and II depends on the electron–ion collision frequency, which in turn depends on electron density and ion cross section. In Paper I, for the sake of simplicity, we considered potassium ion, that is the chemical element with lower ionization energy. In addition, recent studies show that potassium seems to be much more important in the evolution of meteor than previously thought (von Zahn et al. MURAD (1999)). With this assuption, the division between overdense echoes occurs at about $`10^{17}`$ m<sup>-3</sup>. It is worth noting that this border can be moved by considering other elements. But the calculation of particle distribution and evolution in a meteoric plasma will be object of other papers.

The overdense type I echoes derives from total reflection of radio waves (see Fig. 1 for an example). This allows to calculate the height of the meteor in an easy way, as shown in the Paper I. Here we want to present a statistical sample of several meteor showers, for which we have calculated the height. Data will be discussed and compared with available data in literature.

## 2 A simple formula for meteor height

We shortly recall how to calculate the meteor height, as described in Paper I. First, we have to take into account that the plasma has not a definite boundary and then, the incident electromagnetic wave penetrates a little into it before reaching the density necessary to allow total reflection. We can consider this something similar to skin effect in metals.

We can consider a simple geometry, as shown in Fig. 6 of Paper I, and then use the definition of the attenuation $`a`$ in decibel units:

$$a=10\mathrm{log}\frac{|E_i|^2}{|E_r|^2}[\mathrm{dB}]$$
(1)
where subscripts $`i`$ and $`r`$ stand for incident and reflected wave.

¿From the solution of Maxwell’s equations we obtain that, for overdense meteors type I, the electric field is:

$$E=E_0e^{i(kr\omega t)}$$
(2)
where the wave vector in Eq. (2) has the form:

$$k=\beta +i\alpha $$
(3)
We refer to Paper I for explanation of symbols, even though they are commonly used in literature about electromagnetic fields.

We substitute Eqs. (2) and (3) in Eq. (1) and, taking into account that the amplitude of a totally reflected wave is equal to the amplitude of the incident wave, we can obtain an attenuation value of about $`a=20\alpha l\mathrm{log}e`$, where $`l`$ is the path of the wave into the plasma:

$$l=\frac{2\delta }{\mathrm{cos}\varphi }$$
(4)
where $`\varphi `$ is the incidence angle and $`\delta =1/\alpha `$ is the penetration depth. Then, Eq. (1) becomes:

$$a=\frac{40\mathrm{log}e}{\mathrm{cos}\varphi }\frac{17.36}{\mathrm{cos}\varphi }=17.36\mathrm{sec}\varphi [\mathrm{dB}]$$
(5)
In the case of overdense type I (total reflection) the attenuation is simply a function of the angle of incidence.

The Eq. (5) refers to an idealized case. When we deal with real meteors and radio waves, we have to take into account of several factors, i.e. antenna gains, losses in radio receiver and trasmitter, atmospheric absorption, and the distance of reflecting point from trasmitter and receiver. Strictly speaking, Eq. (5) can be considered as the “meteor cross section” in the radar equation.

We can consider common factors in radar theory, as described in Kingsley & Quegan (KINGSLEY (1992)). By means of commonly used values for forward scatter radar, we obtain that the attenuation recorded with our receiver is:

$$a=20\mathrm{log}V_R2[\mathrm{dB}]$$
(6)
where $`V_R`$ is the received signal amplitude \[V\]. ¿From the amplitude of the reflected wave, we can calculate the incidence angle with the Eq. (5).

Therefore, we can calculate the meteor height by considering the geometry of Fig. 2. We can see that:

$$q_1=\frac{L/2}{\mathrm{tan}\varphi }q_2$$
(7)
Taking into account that the Earth’s mean radius $`R_{}`$ is much larger than $`L`$, we can calculate $`q_2`$:

$$q_2=R_{}q_3=R_{}\sqrt{R_{}^2(L/2)^2}\frac{L^2}{8R_{}}$$
(8)
By substituting Eq. (8) in Eq. (7), we obtain the height of the reflection for an overdense meteor type I:

$$q_1\frac{L}{2}\left(\frac{1}{\mathrm{tan}\varphi }\frac{L}{4R_{}}\right)$$
(9)
where $`L=2R_{}\mathrm{sin}\alpha `$. The angle $`\alpha `$ is the half angular distance between receiver and transmitter.

The distance $`L/2`$ plays an important role in the derivation of the above formula. Indeed, from a geometric point of view, the specular reflection in a forward scatter system occurs when the meteor trail lies along a tangent to an ellipsoidal surface, with the transmitter and the receiver stations in foci (Forsyth & Vogan FORSCAT (1955)). This condition is fulfilled by different values of distances of the reflecting point from the source and from the receiver. If we do not know the path source–meteor–receiver, this introduces an uncertainty of about 40 km in the height of the reflecting point (for our system).

One way to overcome this problem is to set up a third station, but if it is not possible, as in our case, we can reduce uncertainties by making heuristic considerations. Indeed, as explained by Forsyth & Vogan (FORSCAT (1955)), the forward scatter system is most sensitive to meteor trails which are nearly horizontal and directed along the transmission line. A meteor perpendicular to the source–receiver line gives an echo that is about five times lower than the case of parallel direction (Forsyth & Vogan FORSCAT (1955)). Therefore, the choice of the reflecting point located closely to the middle of the transmission path appears to be reasonable and, as we shall see, is justified by facts.

## 3 Observations

The general principle of meteor observation by forward scattering of radio waves is the following: a VHF radio receiver (30–100 MHz) is located at a large distance (about 600–2000 km) from a transmitter at the same frequency. Direct radio communication is not possible, owing to the Earth’s curvature, but the meteor allows the commuication over the horizon, by reflecting the transmitted signal.

In this study, observations were carried out mainly by M. de Meyere. His radio receiver is located in Deurle, Belgium (longitude $`3^{}37^{}`$ E, latitude $`51^{}00^{}`$ N), while the trasmitter is located in Sofia (Bulgary). It transmits radio signals at 66.50 MHz all over the day with 10 kW power. In this case, values of $`L`$ and $`\alpha `$ are, respectively, $`1751`$ km and $`7.9^{}`$.

The receiving station consists of a crossed Yagi antenna (4 elements), that is linked to a computer with a digital acquisition interface (150 samples per second, 8 bit resolution). Data are recorded and stored into a file.

Observations of several meteor shower were carried out during several years (see Table 1), with a total number of meteors recorded equal to 13401. The overdense I meteors are 131 (about 1% of the total) and the height distribution for analysed showers are shown in Figs. 36 (each bin is 1 km wide).

Values in Table 1 refer to all echoes recorded during shower days. In order to evaluate also the sporadic background, we have analysed some days in February, without any shower. We have found a mean value for sporadic overdense type I meteors of about 0.2 meteor per hour, so that the contribution of background can be considered negligible.

Measured heights are in the range between 70 and 110 km, in good agreement with typical meteor heights (60–110 km), even though the large part of meteors are in the range 95–110 km. The peaks of distributions are not centered, but are located toward right. However, this seems to be an effect dues to a low number of data. Indeed, the best fit for the observations, calculated with the $`\chi ^2`$ test, is a gaussian distribution (therefore $`\sigma `$ is calculated with standard methods for this type of distribution). In Table 2 values of mean speed (Allen ALLEN (1973)) and height of overdense I meteor are shown.

## 4 Analysis and discussion

Data obtained here show that the mean height is independent from the entry speed of the meteoroids (see Table 2). On the other hand, it is known that the height depends on entry speed of meteoroids: for example, Greenhow & Lovell (GREENHOW (1960)) wrote that the highest speed sporadic meteors, moving at 60–70 km/s, ionize at a mean height of 100 km, whereas those with minimum speed (11.2 km/s) reach a mean height of 85 km.

The theory of radio meteor height was elaborated by Kaiser (1954a , b) and recently Belkovich et al. (BELKOVICH (1999)) proposed some changes, in order to take into account the fragmentation. Kaiser found that the width of the height distribution depends on the atmospheric scale height and the mass distribution of incoming meteoroids. The mean height depends strongly on meteor speed, through two coefficients named $`k_1`$ and $`k_2`$, and depends also on the probability of ionization. Kaiser’s theory refers to the point of maximum ionization, but it is known that in experimental radio observations the height of reflecting point does not necessarily lie in the point of maximum ionization (Greenhow & Lovell GREENHOW (1960)).

If we observe a meteor shower with a given mean speed and mass distribution, the height distribution is related to the length of ionization curve. The point of maximum ionization corresponds to the most probable height. It is worth noting that Kaiser’s theory refers to underdense meteors. Only McKinley (MCKINLEY (1961)) referred to overdense meteors and found no clear dependence on speed. Moreover, he found a two peaks distribution: the main peak is located at about 95 km, and the second one at 106 km.

We can try to explain the differences between our results and data available in literature. The first reason is that we try to analyse overdense meteors, while the large part of published data refer to underdense meteors. It is very interesting to note that for bright meteors, which are surely overdense, the radio height of maximum echo duration is well above the mean height of maximum light, obtained from photographic data (Millman & McKinley MILLMAN (1963)). On the other hand, the situation is reversed for faint meteors.

The reason for this difference is that overdense type I meteors reflect totally the incoming electromagnetic wave. Total reflection is allowed only when plasma frequency is higher than the radio frequency and the collision frequency in the plasma is negligible (see Paper I). These conditions are achieved indipendently from speed of incoming of meteoroid, but it depends on the mass and chemical composition of the body. Once an overdense type I meteor is created, the signal amplitude of the reflected wave remains constant until the collision frequency in the plasma or recombination and attachment processes subtract energy to the plasma frequency. We can say that collective properties of the plasma, which generate the long plateau of overdense type I meteor, “hide” in some way some properties of the incoming meteoroid.

Concerning the two peaks found by McKinley (MCKINLEY (1961)), we note that our distributions show only the secondary peak. This can be explained by taking into account that while McKinley made no distinction between overdense meteors, we have considered only overdense type I meteors.

## 5 Concluding remarks

We have carried out an analysis of several overdense radio echoes, recorded during last years by a radio observer located in Belgium. We have analysed a particular class of overdense meteors (type I) and measured height distributions are in good agreement with previous results obtained by McKinley (MCKINLEY (1961)), even though only for the secondary peak. We suppose that the first peak in McKinley’s work should be due to overdense type II meteors, while the secondary peak, recorded also by our system, appeared to be due to overdense type I meteors.

We think that collective properties of the meteoric plasma (Langmuir oscillations) hide some characteristics of the originary cosmic body, specifically there is no clear dependence on speed. Further study, mainly theoretical and able to take into account collective properties of plasma, are required to assess the particle dynamics in the meteor.

It should be noted that our studies were carried out with an amateur forward scatter system and we have no full control on it. Moreover, heuristic considerations were introduced in order to minimize uncertainties, but results showed that they were justified by facts. The agreement with previous works, with other techniques, supports our conclusions. The future availability of a full forward scatter system would be of great help in more detailed studies.

###### Acknowledgements.

LF wishes to thank P. Farinella for useful discussions. Authors are also grateful to the referee, P. Pecina, for useful comments. This research has made use of *NASA’s Astrophysics Data System Abstract Service*.


format: text

text:
# 1 INTRODUCTION

## 1 INTRODUCTION

There is a great current interest in dynamics described by the complex Ginzburg-Landau equation (CGLE) , which is a generic nonlinear model with various physical applications, including binary-fluid thermal convection , semiconductor lasers , nonlinear optical fibers (see Ref.  and references therein), etc. In most applications, the boundary conditions (BC) supplementing the CGLE are an inherent type of the underlying physical model. In particular, in the case of the thermal convection and laser cavities, the simplest and most fundamental type of BC are the periodic conditions, as they make it possible to study the intrinsic dynamics of CGLE avoiding obscuring effects of wave reflection from edges.

It should be stressed that the dynamics of the CGLE with periodic BC is principally different from that of an infinite system, as various solutions (even chaotic) which are stable under periodic BC may be easily destabilized in a much longer system . Another important difference from the infinite limit is that the equation subject to periodic BC is equivalent to an infinite chain of evolution equations for amplitudes of the corresponding Fourier expansion. Using this approach, it has been proved that a finite-dimensional attractor exists in the system’s phase space  (which is true even in the case of two-dimensional CGLE ). Parallel to this, different technical versions of an approximation based on a finite-mode (Galerkin) truncation of the above-mentioned infinite chain of equations were developed, helping to understand the onset of dynamical chaos in the CGLE . Although it was argued that the Galerkin truncation does not always describe adequately the transition to chaos in CGLE , it has been demonstrated that, generally, the finite-mode truncation of the CGLE with periodic BC provides for quite accurate results .

For some particular cases, a chain of bifurcations leading to chaos in CGLE under periodic BC was found by means of direct numerical simulations : the chain features transitions limit-cycle $``$ two-torus $``$ three-torus $``$ chaos, $`n`$-torus being a quasi-periodic dynamical regime with $`n`$ incommensurate frequencies. Many other nonlinear modes are unstable at onset (i.e., when they appear), but they may be later restabilized through higher-order bifurcations .

In this paper, we study analytically and numerically the dynamical behavior and, in particular, stable solutions of the cubic CGLE, subjected to periodic BC, without diffusion, while diffraction (which formally corresponds to an imaginary part of the diffusion coefficient) is present. This situation is typical for laser cavities and other nonlinear optical problems in the spatial domain (see Ref.  and ), as light beams are subject to transverse diffraction but not diffusion.

A general form of the one-dimensional CGLE with cubic nonlinearity and without diffusion is

$$iu_t+\frac{1}{2}u_{xx}\left(1+iR\right)\left|u\right|^2u=iu.$$
(1)
Here, $`u=u(x,t)`$ is the complex field which depends on the time $`t`$ and the spatial coordinate $`x`$, $`R`$ is a real positive parameter that represents cubic losses and/or gain saturation. The coefficients of the spatial dispersion, the Kerr nonlinearity and the linear gain (the latter one accounts for the term on the right-hand side of Eq. (1)) are all normalized to 1. In nonlinear optics, where Eq. (1) finds its most important application, it describes the spatial evolution of the electromagnetic field envelope $`u(t,x)`$ in a planar waveguide or laser cavity with linear gain, cubic loss and Kerr nonlinearity. In this case, the evolutional variable $`t`$ is, in fact, not time but the propagation distance (which is frequently designated $`z`$), and the second term in Eq. (1) accounts for the transverse diffraction .

Equation (1) also describes the temporal or spatial evolution of the electric field in a (1+1)-dimensional semiconductor laser medium with gain saturation. The second term in Eq. (1) then accounts for the diffraction in space or dispersion in the temporal domain, and the cubic conservative term represents a normalized anti-guiding or linewidth-enhancement factor.

In either case, we supplement Eq. (1) with periodic BC, , where $`L`$ is the system’s length (in the transverse direction, in the case of the planar waveguide). In its literal form, Eq. (1) with periodic BC applies to a planar waveguide closed into a cylindrical surface. Such a configuration can be easily implemented in Vertical Cavity Surface-Emitting Lasers (VCSEL) . However, the most widespread application of Eq. (1) subject to periodic BC is the description of a ring cavity laser, in which some ingredients of the system (e.g., the saturable amplifier and nonlinear waveguide) may be separated. In this context, Eq. (1) has the meaning of an average master equation that govern the evolution of the electromagnetic field in the cavity. Thus, the model is controlled by two dimensionless parameters, $`R`$ and $`L`$. Despite a number of results obtained in the above-mentioned works for CGLE of a more general type with periodic BC, the case of zero diffusion is fundamentally important for applications to nonlinear optics and must be considered separately.

Equation (1) has an obvious continuous-wave (CW), i.e., stationary spatially uniform solution,

$$u_0=R^{1/2}\mathrm{exp}\left(iR^1t\right).$$
(2)
A straightforward stability analysis, taking into account the periodic BC, shows that the CW solution (2) is unstable to small perturbations if the length $`L`$ exceeds a critical value:

$$L_{\mathrm{cr}}=\pi \sqrt{R}.$$
(3)
A primary objective of this work is to study the model’s dynamics in the case $`L>L_{\mathrm{cr}}`$. To this end, we solved Eq. (1) numerically, using the uniform solution (2) perturbed by a small-amplitude harmonic with the fundamental wave-number ($`2\pi /L`$) as initial conditions. We have found that the long-time evolution of the solutions exhibits a wide variety of behaviors (stable, oscillating, chaotic), depending on values of the parameters $`R`$ and $`L`$.

The rest of the paper is organized as follows. In section 2, we display a phase diagram, which summarizes the dynamical behavior of the system in terms of the ($`R,L`$) parametric plane. In this section, we also find, in an asymptotic analytical form, the first spatially nonuniform stable stationary solution which sets in after the stability loss of (2), and compare it to numerical results. Although the procedure leading to the analytical solution is quite simple, the solution itself appears to be a new one. In section 3, we study higher-order stable stationary solutions which have a multi-humped structure. In section 4, we study various oscillatory dynamical regimes and section 5 concludes the paper.

All these dynamical regimes admit straightforward interpretation in terms of field patters in the laser cavity which is modeled by Eq. (1). In fact, as it was already mentioned above, we study in detail two vast classes of regular stable solutions: stationary and oscillatory (we also find chaotic dynamical regimes, but, from the viewpoint of applications, they are interesting mainly in the sense of choosing parameter values in order to avoid chaos). The physical relevance of both types of regular dynamical regimes is quite obvious. In particular, we will find different parametric regions that feature stable stationary states with one or several humps (local maxima). This finding has straightforward applications in term of the laser cavity, making it possible, for instance, to simultaneously generate several output light beams. Stable oscillatory regimes are relevant too, as they can help the laser to generate a temporally modulated (pulsed) beam with various types of the modulation, depending on the quasi-harmonic or strongly anharmonic character of the oscillations (we find stable solutions of both types).

## 2 A DIAGRAM OF DYNAMICAL REGIMES AND BASIC STATIONARY STATES

Figure 1 displays a diagram of different dynamical regimes in the parametric plane ($`L,R`$) generated by numerical simulations of Eq. (1). The figure shows several regions of regular (stationary or oscillatory) solutions, which are separated by chaotic layers. Regular solutions existing in different regions exhibit different spatial patterns. Each region is divided into sub-regions of stationary and oscillating solutions.

In region I, which is , the uniform solution (2) is stable. Direct simulations demonstrate that, indeed, there is no other attractor in this region. The first stable stationary non-uniform solution found beyond the boundary (3) is shown in Fig. 2. Close to the boundary, this solution can be approximated analytically as a combination of a uniform field and the first two spatial harmonics:

$$u(x,t)=R^{1/2}\mathrm{exp}\left(iR^1ti\omega _2t\right)\left[1+a_0+a_1\mathrm{cos}\left(kx\right)+a_2\mathrm{cos}\left(2kx\right)\right]$$
(4)
where $`k=2\pi /L`$, $`a_1`$, $`a_0`$ and $`a_2`$ are, respectively, small amplitudes of the fundamental, zeroth, and second harmonics, $`\omega _2`$ is a frequency shift, and it is expected that $`a_0`$, $`a_2`$ and $`\omega _2`$ are $`a_1^2`$.

To proceed with the analytical consideration of the solution (4), we eliminate $`a_0`$, $`a_2`$ and $`\omega _2`$, substituting Eq. (4) into Eq. (1) and equating terms at order $`a_1^2`$. We notice that the ansatz (4) allows one to assume that $`a_0`$ is purely real. Obviously, is real too, while $`a_1`$ and $`a_2`$ may be complex:

$$a_1a_{11}+ia_{12},a_2a_{21}+ia_{22}.$$
(5)
Substituting (5) into (4) and solving for both real and imaginary parts, we obtain the following expressions at order $`a_1^2`$:

$`\omega _2`$ $`=`$ $`\left(1+{\displaystyle \frac{1}{R^2}}\right)a_{11}a_{12},a_0={\displaystyle \frac{3}{4}}a_{11}^2{\displaystyle \frac{1}{4}}a_{12}^2{\displaystyle \frac{1}{2R}}a_{11}a_{12},`$
$`a_{21}`$ $`=`$ $`{\displaystyle \frac{1}{4}}a_{11}^2+{\displaystyle \frac{1}{12}}a_{12}^2{\displaystyle \frac{R}{6}}a_{11}a_{12},a_{22}={\displaystyle \frac{R}{4}}a_{11}^2+{\displaystyle \frac{R}{12}}a_{12}^2\left({\displaystyle \frac{R^2}{24}}{\displaystyle \frac{1}{8}}\right)a_{11}a_{12}.`$ (6)
Next, we equate coefficients in front of the fundamental harmonic (obtained after the substitution of (4) into (1)) to find the amplitude $`a_1`$ at order $`a_1^3`$. It is convenient to introduce a small positive parameter $`\delta `$ (overcriticality) which measures proximity of the system to the bifurcation point (3),

$$\delta 4R^1k^2$$
(7)
(recall $`k=2\pi /L`$). The amplitude a1 can then be found as a function of $`\delta `$,

$$a_{12}=Ra_{11},a_{11}^2=\frac{3R}{R^49R^2+30}\delta .$$
(8)
Figures 3 and 4 present comparison between the analytical solution (4) - (8) and a numerical one for different values of $`R`$ and $`L`$. Up to relatively large values of the overcriticality ($`\delta =0.30`$ in Fig. 3), fairly good agreement between the numerical and analytical solutions is evident. A very good agreement was also found between the analytically calculated frequency shift and the shift that was extracted from numerical simulations. With further increase of $`\delta `$, a deviation between the analytical and the numerical solutions becomes visible (at $`\delta 0.40`$ in Fig. 4).

The parameter range, in which this stationary solution is stable, is a part of region II in Fig. 1. Although the region is narrow, it has an internal structure, being divided into three subregions, hosting the above stationary solution and oscillatory ones. A detailed consideration of the subregions is given in section 4. Above the region II (at larger $`L`$ and $`R<2.6`$), the simulations reveal solutions which are chaotic in time.

## 3 HIGHER-ORDER STABLE STATIONARY SOLUTIONS

Another region partly filled with stable stationary solutions is found as $`L`$ is increased, namely, region III in Fig. 1. In contrary to the solutions considered above, the profile of the solution in the region III, a typical example of which is depicted in Fig. 5, has two maxima (with different values) and two minima, being similar to solutions discussed, in the framework of a similar model, in Ref. . In this region, the amplitude of the second harmonic is not small in comparison with that of the fundamental harmonic, giving rise to the double-humped profile of the solution. A minimum value of the nonlinear-loss coefficient for which region III exists is $`R=1.06`$, the corresponding value of the system’s length being $`L=5.9`$. Regions II and III merge at $`R>2.6`$ and form a single region (see Fig. 1). The dotted line inside the united II - III region in Fig. 1 indicates a border between the above-mentioned single-humped and double-humped stationary solutions. This border is a direct continuation of the upper border of region II at lower values of $`R`$, strongly suggesting that the transition between the single- and double-humped patterns is caused by destabilization of the single-humped solution. In fact, region III, as well as II, is also divided into several subregions exhibiting stationary and oscillating regular solutions.

Above region III, there is a well-defined stripe in which the behavior of the system is chaotic. For higher values of $`L`$, still another range of stable stationary solutions, region IV in Fig. 1, is revealed. Stationary solutions in this region also have two lobes, like in region III, but its second-harmonic component is larger, giving rise to a profile which is similar to the solution discussed in Ref.  (a cnoidal-wave-type, see Fig. 6). It should be noted, however, that region IV is rather narrow.

The next stability range for stationary solutions is found directly above region IV, at $`R>1.72`$ and $`L>11.8`$ (region V in Fig. 1). In this region, stable solutions have four maxima (two pairs of local maxima with of different heights), see Fig. 7. It seems plausible that this solution was obtained via a doubling bifurcation of the double-humped solution found above in region II, especially since the minimal system length for this region ($`L=11.8`$) is exactly twice the minimal length for region III ($`L=5.9`$).

This pattern of changing the structure of stationary solutions develops at larger values of $`L`$, with stability regions featuring stationary solutions that appear to be produced by tripling and quadrupling the basic stationary solution found in region III (see Fig. 5). As $`L`$ is increased, reaching region VI in Fig. 1, a new stable stationary solution of Eq. (1) appears, characterized by local maxima with three different amplitudes, see Fig. 8. The minimum value of $`R`$ for which this solution appears is $`R=1.55`$, at $`L=18`$. This solution evolves gradually from the four-lobed solution found in region V and described above, see Fig. 7. As $`L`$ is increased, the two small maxima move outwards; simultaneously, two additional humps evolve between the small maxima and the central high maxima (see Fig. 9). Eventually, these two humps evolve into two additional maxima.

## 4 OSCILLATORY SOLUTIONS

As it was mentioned before, each region containing stable regular solutions is actually divided into several subregions that exhibit a wide variety of dynamical behaviors. As an example, we present here a detailed description of the intrinsic structure of region II.

Figure 10 illustrates the intrinsic structure of region II for $`0.2<R<2.2`$. Although this region is relatively narrow, it is divided into three subregions. In the first subregion (IIa), the established solution is always the stationary one that was described in section 2. For larger values of $`L`$ (subregion IIb), the solution does not remain stationary; instead, it becomes an oscillatory function of $`t`$ (see Fig. 11).

Further increase of $`L`$ results in a change of the character of the oscillations in the subregion IIc: while the oscillations in the subregion IIb were quasi-harmonic, in the subregion IIc they acquire a more complex anharmonic structure. Clear characterization of the oscillations of different types is provided by examining the corresponding time dependence of the average field power, which is defined in a straightforward way,

$$P(t)=L^1_0^L\left|u(x,t)\right|^2𝑑x,$$
(9)
and also by the temporal Fourier transform (power spectrum) of the time-dependent power (9). In particular, in Fig. 12 it is easy to see distinction between the aforementioned quasi-harmonic and strongly anharmonic oscillations, in terms of the dependences $`P(t)`$.

Figure 13 shows the power spectrum of $`P`$ obtained at $`R=0.7`$, while varying $`L`$ between $`L=2.9`$ and beyond the upper border of region II ($`L=3.15`$). In this figure, the dc (zero-harmonic) component of the Fourier transform was removed in order to stress the dynamical behavior. As $`L`$ crosses the upper border of region IIa, two symmetric harmonics with opposite frequencies appear, indicating that the solution is oscillating quasi-harmonically. As $`L`$ is increased further, the frequency of the oscillations drops to a minimum value, after which several higher-order harmonics appear, and the dominant frequency increases with $`L`$. At the point with the lowest value of the oscillation frequency, a transition between two different dynamical regimes, viz., quasi-harmonic oscillations whose frequency decreases with $`L`$ and complex multi-frequency oscillations whose dominant frequency increases with $`L`$, takes place.

Further increase of $`L`$ eventually results in a transition to dynamical chaos. A similar behavior is observed when decreasing $`R`$ at constant $`L`$.

It should be emphasized that all the other stability regions marked in Fig. 1 are also divided into subregions which exhibit various oscillatory dynamical regimes, see examples shown in Fig. 14 for region V, and in Fig. 15 for region VI. However, a detailed analysis of these subregions is beyond the scope of the present paper.

## 5 CONCLUSIONS

We have studied the dynamical properties of the solutions of the complex cubic Ginzburg-Landau equation without diffusion, which is subject to periodic boundary conditions. This equation is a fundamental model of nonlinear-optical cavities characterized by saturable gain. Systematic investigation of its dynamics is necessary, as it is essentially different from both a more general Ginzburg-Landau equation with its entire coefficients complex (i.e., including also diffusion), and from its infinite-length counterpart.

We found that the model exhibits a complex structure, including regions of stable stationary, oscillatory, and chaotic solutions in its parameter plane. For a case of strong nonlinear loss, the system features a multi-layered structure in the parameter plane, with a different dynamical behavior in each layer. A boundary between adjacent layers implies the existence of a bifurcation destabilizing the solution in one of layers. This could be easily seen by examining the boundary between regions II and III in Fig. 1, i.e., the upper boundary of region II at $`R>2.6`$: it is a direct continuation of the upper boundary of region II which, at $`R<2.6`$, separates it from a chaotic zone. It should be stressed that we did not observe hysteresis (overlapping) between stable solutions of different types.

We have found an asymptotic analytical single-humped solution for the case when the system is close to the first bifurcation destabilizing the spatially uniform CW solution. Even when the system is not really close to the bifurcation point, the analytical solution is found to be in good agreement with numerical results. Numerically, we have also found a series of stable stationary solutions of a multi-humped type.

One of the authors (B.A.M.) appreciates a grant ”Window-on-Science”, provided by the European Office for Research and Development of the US Air Force.

Figure Captions

Fig. 1. The dynamical phase diagram of the system, showing regions of regular and chaotic behavior in the parameter plane ($`R,L`$), see detailed description of the regions in the text.

Fig. 2. A stable stationary nonuniform solution found for $`R=1.0`$ and $`L=3.4`$ (inside region II of Fig. 1), which has a single maximum and a single minimum. In this figure and in Figs. 3 - 9 below, shown is $`|u(x)|`$ in the interval $`0<x<L`$.

Fig. 3. Comparison between the analytical (circles) and the numerical (solid) stationary solutions for $`R=0.7`$ and $`L=2.7`$. At these values of $`R`$ and $`L`$, Eq. (7) gives the value of the overcriticality parameter $`\delta =0.30`$, and Eq. (3) yields the critical value beyond which the nontrivial stationary solution may exist, $`L_{\mathrm{cr}}=2.628`$.

Fig. 4. Comparison between the analytical (dashed) and numerical (solid) stationary solutions for $`R=1.4`$ and $`L=4.01`$. At this value of $`R`$, one has $`L_{\mathrm{cr}}=3.717`$ and $`\delta `$ = 0.40.

Fig. 5. A typical profile of the stable stationary solution found in region III, for $`R=2.0`$ and $`L=6.0`$. The solution exhibits two maxima and two minima.

Fig. 6. The profile of the stable stationary solution found solution in region IV for $`R=3.2`$ and $`L=11.8`$, featuring a cnoidal-wave-like solution with two maxima.

Fig. 7. The profile of the stable stationary solution found in region V for $`R=2.5`$ and $`L=12.0`$, which exhibits a double-period solution having local maxima and minima with different amplitudes.

Fig. 8. The profile of the stable stationary solution obtained in region VI for $`R=1.8`$ and $`L=18.8`$, with three maxima having different amplitudes.

Fig. 9. The evolution of the stable stationary solution observed when crossing from region V into VI, for $`R=2.0`$ and $`L=13.2`$ (solid), $`17.0`$ (dotted), $`19.0`$ (dashed), $`21.0`$ (dashed-dotted).

Fig. 10. The inner structure of region II which reveals three subregions that feature, respectively, stable stationary patterns and quasi-harmonic or strongly anharmonic oscillatory solutions, as it is described in the text in detail.

Fig. 11. The oscillatory solution as a function of $`t`$, found in region IIb (for $`R=0.3`$ and $`L=1.85`$).

Fig. 12. Comparison between the time dependence of the average field power $`P(t)`$ in regions IIb (for $`R=0.7`$ and $`L=2.97`$, solid line) and IIc (for $`R=0.7`$ and $`L=3.02`$, dashed line).

Fig. 13. The evolution of the power spectrum of $`P(t)`$ inside the interval $`2.9<L<3.15`$ at fixed $`R=0.7`$, indicating at the existence of two different types of oscillatory regimes.

Fig. 14. Examples of the temporal evolution of the field corresponding to different types of oscillatory dynamical regimes in region V: (a) $`R=4.0`$, $`L=17.5`$; (b) $`R=3.0`$, $`L=17.4`$.

Fig. 15. Field evolution in time in region VI, for $`R=1.7`$ and $`L=16.1`$.


text:
# An Upper Limit on the Mass of a Primordial Star due to the Formation of an HII Region: The Effect of Ionizing Radiation Force

## 1 Introduction

Stars are supposed to be born in the dense regions of gas clouds. The theory of the gravitational fragmentation of the gas clouds describes how these dense regions are formed inside the parental gas clouds. The typical mass scale of the dense region can be predicted by the Jeans mass of the cloud at the time of the fragmentation. However, a perturbation with a mass scale larger than this Jeans mass is always unstable. Therefore a dense region that is more massive than the typical Jeans mass can form depending on the initial fluctuation spectrum in the cloud. In this sense, the upper limit for the mass of stars cannot be obtained from the analysis of the gravitational fragmentation.

Once, the origin of the observed upper limit of stellar mass (around $`100M_{\mathrm{}}`$) was attributed to the instability of massive stars due to the $`\epsilon `$-mechanism (e.g., Schwarzshild & Härm 1959). Since later studies revealed that this mechanism leads to moderate mass-loss rather than disruption of the stars (e.g., Appenzeller 1970), formation process, instead of the stability of massive stars, has been believed to limit the upper bound for the stellar mass (e.g., Nakano, Hasegawa, & Norman 1995).

The star formation process is an accretion of ambient matter onto a protostar (stellar core) forming inside a protostellar cloud. In this scenario, the final mass of a star depends on how much mass the protostar can acquire. Although it is not yet clear what mechanism stops the accretion, increasingly strong feedback from a massive forming star is likely to halt the accretion.

Larson & Starrfield (1971) were first to propose that a firm upper limit on the stellar mass is provided by the formation of an HII region around a forming star, as well as by the radiation force acting on dust grains in the accreting envelope. The formation of an HII region prevents further accretion in the following way: when the HII region reaches outer layers of a protostellar cloud, the temperature and then the pressure support surge by some large factor, so that the further infall of material is immediately halted. Larson & Starrfield (1971) concluded that both of those mechanisms set an upper limit of about 50 M on the masses of Pop I stars forming out of the present interstellar medium. The radiation force onto dust grains was found by subsequent studies to set more stringent mass upper limit (Appenzeller & Tscharnuter 1974; Kahn 1974; Yorke & Krügel 1977; Wolfire & Cassinelli 1987). Among them, Wolfire & Cassinelli (1987) suggested that even stars as massive as about 10 M cannot form by accretion owing to this mechanism unless dust is significantly depleted.

However, for so-called Pop III stars forming out of metal-free gas, this mechanism of the radiation force does not work because of the absence of dust grains. The $`\epsilon `$-mechanism drives even milder mass-loss for those stars than Pop I massive stars (Baraffe, Heger, & Woosley 2001). According to recent simulations of the fragmentation of primordial clouds, the mass scale of fragments (protostellar clouds) is as large as $`10^3M_{\mathrm{}}`$ (Bromm, Coppi, & Larson 1999; Abel, Bryan, & Norman 2000; Tsuribe 2001). Therefore, the growth of the HII region plays a crucial role in limiting the maximum mass of a Pop III star, if the upper limit by this mechanism is less than the mass scale of fragments. This motivates us to study here the upper mass limit of Pop III stars due to the formation of an HII region.

In discussing the formation of an HII region, often the free-fall assumption has been imposed on the flow in the accreting envelope (e.g., Yorke 1986). The role of the momentum transfer to the gas due to ionizing radiation has been neglected. In the context of galaxy formation, however, its significance in dynamical and thermal evolution of the intergalactic medium has been pointed out by Haehnelt (1995). He showed, for the collapse of subgalactic clouds of $`10^{10}M_{\mathrm{}}`$, even a radiation-pressure-driven bounce is possible. Taking this effect into account, in this paper, we study the formation of an HII region in the accreting envelope of a forming massive star. Although this effect is negligible as long as the radius of the HII region is small, it will turn out that inclusion of the ionizing radiation force alters the later evolution of the HII region. In the course of the expansion of the HII region, the flow becomes slower than the free-fall rate owing to the ionizing radiation force. The growth of the HII region is strongly suppressed by this effect. Consequently, the mass upper limit imposed by the expansion of the HII region is increased by a large factor. In particular, the formation of an HII region does not set any firm mass upper limit for Pop III stars.

In Section 2, we briefly summarize the relevant aspects of the formation of an HII region in a free-falling envelope around a newborn star. The effect of ionization radiation force is included in Section 3. Finally, we provide a summary and discussion in Section 4.

## 2 Formation of an HII Region in a Free-Falling Envelope

In this section, we review the expansion law of an HII region around a newly formed star (e.g., Yorke 1986) and discuss the resultant upper mass limit of a forming star assuming that the accreting flow is free-falling.

Suppose an accreting flow onto a star of mass $`M_{}`$ is taking place at an accretion rate $`\dot{M}`$. The mass accretion rate $`\dot{M}`$ is determined by the temperature of the protostellar cloud and does not change very much except during the very late stage of the accretion phase where most of the mass of the envelope has accreted onto the central star. The relation between $`\dot{M}`$ and the sound speed $`c_\mathrm{s}`$ in the protostellar cloud is given by

$$\dot{M}=C_{\mathrm{acc}}\frac{c_\mathrm{s}^3}{G},$$
(1)
where $`C_{\mathrm{acc}}`$ is a non-dimensional number (see, e.g., Whitworth & Summers 1987). When the numerical value of $`\dot{M}`$ is needed, we use $`10^3M_{\mathrm{}}\mathrm{yr}^1`$ corresponding to $`C_{\mathrm{acc}}=1`$ and 400K which is a typical temperature of metal-free protostellar clouds (e.g., Yoneyama 1972; Carlberg 1981; Palla, Salpeter, & Stahler 1983).

Here and throughout this paper, we assume that the accreting flow is spherically symmetric and steady. Then, the continuity equation yields

$$\dot{M}=4\pi r^2\rho u,$$
(2)
$`\rho `$ and $`u`$ are the density and velocity in the accreting flow at radius $`r`$. In addition, we assume in this section that the accreting envelope is free-falling:

$$u=\sqrt{\frac{2GM_{}}{r}}.$$
(3)
Using equations (2) and (3), we obtain the density distribution

$$\rho =\frac{\dot{M}}{4\pi \sqrt{2GM_{}}}r^{3/2}.$$
(4)
In the HII region, we assume that the ionization and recombination are in equilibrium and that all the re-emitted ionization photons are consumed locally (a so-called on-the-spot approximation). Suppose $`Q`$ ionizing photons are emitted per unit time from a sphere of radius $`R_{\mathrm{in}}`$ surrounding the star. We regard the stellar surface at $`R_{\mathrm{in}}`$ as the base of the HII region for the moment (see discussion in §4). For the sake of simplicity, we assume that the gas is purelycomposed of hydrogen and is fully ionized inside the HII region (i.e., the number density of H+ $`n(\mathrm{H}^+)=n(e)=\rho /m_\mathrm{p}`$), hereafter. We obtain the radius $`R_{\mathrm{II}}`$ of the HII region from the ionization equilibrium:

$$Q=_{R_{\mathrm{in}}}^{R_{\mathrm{II}}}\alpha n(\mathrm{H}^+)n(e)𝑑V+\frac{\dot{M}}{m_\mathrm{p}},$$
(5)
where the recombination coefficient $`\alpha =\alpha _\mathrm{B}`$ from the on-the-spot approximation. The first term on the right hand side means the recombination rate in the HII region. The second term means the neutral flux entering into the HII region. Note that no HII region forms if the second term on the right-hand side is larger than $`Q`$. In this paper, we neglect the $`\dot{M}/m_\mathrm{p}`$ term, assuming $`Q\dot{M}/m_\mathrm{p}=3.8\times 10^{46}\dot{m}_3(\mathrm{s}^1)`$, where $`\dot{m}_3\dot{M}/10^3M_{\mathrm{}}\mathrm{yr}^1`$. Using the relation (4), we obtain

$$Q=\frac{\alpha \dot{M}^2}{8\pi Gm_{\mathrm{p}}^{}{}_{}{}^{2}M_{}}\mathrm{ln}(R_{\mathrm{II}}/R_{\mathrm{in}}).$$
(6)
Then,

$$R_{\mathrm{II}}=R_{\mathrm{in}}\mathrm{exp}(Q/Q_{\mathrm{crit}}^{(\mathrm{FF})}),$$
(7)
where

$`Q_{\mathrm{crit}}^{(\mathrm{FF})}`$ $`=`$ $`{\displaystyle \frac{\alpha \dot{M}^2}{8\pi Gm_{\mathrm{p}}^{}{}_{}{}^{2}M_{}}}`$ (8)
$`=`$ $`1\times 10^{51}(M_{}/100M_{\mathrm{}})^1\dot{m}_{3}^{}{}_{}{}^{2}(\mathrm{s}^1).`$ (9)
This shows that if $`Q`$ exceeds $`Q_{\mathrm{crit}}^{(\mathrm{FF})}`$, the HII region expands exponentially. Soon after, the HII region reaches an outer layer whose gravitational binding energy is lower than the thermal energy of ionized matter. This results in the immediate expulsion of the outer layer and the halting of the accretion.

When does the ionizing photon emissivity $`Q`$ of the star reach the critical value $`Q_{\mathrm{crit}}^{(\mathrm{FF})}`$? To be concrete, we consider here the case of Pop III stars. The ionizing photon emissivity

$$Q1.6\times 10^{50}(M_{}/100M_{\mathrm{}})(\mathrm{s}^1)$$
(10)
for very massive ($`M_{}300M_{\mathrm{}}`$) Pop III stars in the zero-age main sequence (ZAMS; Bromm, Kudritzki, & Loeb 2001). Comparing (10) with (9), we obtain the critical stellar mass above where $`Q`$ exceeds $`Q_{\mathrm{crit}}^{(\mathrm{FF})}`$:

$$M_{}3\times 10^2\dot{m}_3M_{\mathrm{}}.$$
(11)
The formation of a more massive star is increasingly difficult because of the growth of the HII region.

Recall here that the exponential expansion law of the HII region (7), which plays a key role in determining the upper limit of stellar mass, is a direct consequence of the density distribution (4) and then the free-fall assumption (3). However, the gas is subject to momentum transfer due to ionizing radiation and becomes slower than the free-fall rate within the HII region. This modifies the density law (4) and the expansion law of the HII region (7). Taking this effect into account, we study, in the next section, the growth of an HII region in the accreting envelope of a newborn massive star.

## 3 Effect of Radiation Force on Expansion of an HII Region

Hydrogen atoms in the HII region are ionized both by ionizing photons emitted directly by the star, and by re-emitted ionizing photons due to the recombination to the ground level. While the former is radially outward-directed, the re-emitted photons become almost isotropic and contribute little to the radiation force in the radial direction. For the calculation of radiation force, we must subtract the contribution from the ionization by locally emitted photons from the total ionization rate $`\beta n(\mathrm{H})`$. The rate of the former equals the recombination rate to the ground level $`\alpha _1n(\mathrm{H}^+)\mathrm{n}(\mathrm{e})`$ (per unit volume) according to the on-the-spot approximation, while the latter equals the total recombination rate $`\alpha _\mathrm{A}n(\mathrm{H}^+)n(e)`$ by ionization equilibrium. Hence, the radiation force by ionization radiation per unit mass is given by

$$f_{\mathrm{ion}}=\frac{h\nu _{\mathrm{ion}}}{c}\alpha n(\mathrm{H}^+)n(e)/\rho ,$$
(12)
where, $`\alpha =\alpha _\mathrm{B}`$ again. Since we are considering stellar UV sources of thermal spectra, the mean energy of ionizing photons is $`h\nu _{\mathrm{ion}}h\nu _0=13.6`$eV. For later convenience, we define the Haehnelt’s radius

$$R_\mathrm{H}=\frac{\alpha h\nu _{\mathrm{ion}}}{4\pi cGm_\mathrm{p}^2}26\mathrm{p}\mathrm{c},$$
(13)
first introduced by Haehnelt (1995), which is a fundamental length scale for the balance between ionizing radiation force and gravity.111 Comparing radiation force $`f_{\mathrm{ion}}=4\pi GR_\mathrm{H}\rho `$ from a point source with gravity $`f_{\mathrm{grav}}=G(4\pi \rho r^3/3)/r^2`$ exerted by homogeneously distributed matter with density $`\rho `$ around the source, Haehnelt (1995) found the radiation force dominates the gravity at $`r<3R_\mathrm{H}`$. Finally, the expression of the ionizing radiation force becomes

$$f_{\mathrm{ion}}=4\pi GR_\mathrm{H}\rho .$$
(14)
Note that we have neglected other photon extinction processes such as electron scattering, photodissociation of H-, etc. This assumption is valid as long as the HII region outside $`R_{\mathrm{in}}`$ is optically thin to those extinction processes.

We consider here an HII region with the steady radius $`R_{\mathrm{II}}`$ determined by the ionization equilibrium (5) as in the previous section. The I-front moves outward slowly as the stellar mass, and then the ionization photon emissivity from the star, increase. If the flow velocity across the ionization front (I-front) is more than twice the sound speed of ionized matter, the I-front becomes a (weak) R-type (Kahn 1954). Since we consider the steady radius of the HII region, the velocity across the I-front equals that of the flow at the edge of the HII region. As long as the HII region is not large, the flow velocity is supersonic and fast enough to keep the I-front an R-type. In fact, as we will show later, the halting of accretion occurs while the I-front is still R-type. Hence, we treat here only the HII region bounded by R-front and do not consider later evolution.

Including the ionization radiation force, the momentum equation for the steady flow yields

$$u\frac{du}{dr}=\frac{GM_{}}{r^2}+4\pi GR_\mathrm{H}\rho $$
(15)
in the HII region ($`R_{\mathrm{in}}rR_{\mathrm{II}}`$). The thermal pressure term can be neglected since we consider only supersonic flow.

Using the density-velocity relation in the steady accretion

$$\rho =\frac{\dot{M}}{4\pi r^2}(u)^1,$$
(16)
and introducing the kinetic energy density of the flow $`Ku^2/2`$, we obtain from equation (15)

$$\frac{dK}{dr}=\frac{GM_{}}{r^2}\left(1\sqrt{\frac{K_{\mathrm{crit}}}{K}}\right).$$
(17)
Here we define the critical kinetic energy $`K_{\mathrm{crit}}`$ and the critical flow velocity $`u_{\mathrm{crit}}`$:

$$K_{\mathrm{crit}}=\frac{1}{2}u_{\mathrm{crit}}^2=\frac{1}{2}(R_\mathrm{H}/t_{\mathrm{acc}})^2,$$
(18)
and the accretion timescale

$$t_{\mathrm{acc}}=M_{}/\dot{M}.$$
(19)
Given the value of $`K(=K_{\mathrm{II}})`$ at the outer boundary $`R_{\mathrm{II}}`$ of the HII region, we can solve the momentum equation (17) inward. Depending on whether the initial value $`K_{\mathrm{II}}`$ is larger than, equal to, or smaller than the critical value $`K_{\mathrm{crit}}`$, solutions of the equation (17) are classified into three types. In Figure 1, the solutions for $`K_{\mathrm{II}}/K_{\mathrm{crit}}=1.1,1.01,1,0.99`$, and 0.9 are shown. If $`K_{\mathrm{II}}`$ is larger than the critical value $`K_{\mathrm{crit}}`$, the right hand side of (17) is always negative. Then $`K`$ becomes larger and larger as the gas flows inward: The flow is accelerated. We call these types of solutions accelerating solutions. If $`K_{\mathrm{II}}=K_{\mathrm{crit}}`$, $`K`$ remains always at this value (the uniform velocity solution),

$$u(r)=u_{\mathrm{crit}}=\frac{R_\mathrm{H}}{t_{\mathrm{acc}}}$$
(20)
In other words, the ionizing radiation force exactly balances the gravity everywhere in the HII region. From equations (2), (18) and (19), the density distribution of the uniform velocity solution is given by

$$\rho =\frac{M_{}}{4\pi R_\mathrm{H}}r^2.$$
(21)
Similarly, if $`K_{\mathrm{II}}<K_{\mathrm{crit}}`$, the flow continues to decelerate (decelerating solutions). As can be seen in Figure 1, the deceleration is so rapid that the accretion is halted by the radiation force without reaching the star shortly after the critical flow is attained.

As the outer boundary condition, we assume that the flow is free-fall outside the ionizing region. We further assume that the flow velocity is continuous at the I-front, since velocity jump at the R-front is small. Then,

$$K_{\mathrm{II}}=\frac{GM_{}}{R_{\mathrm{II}}},$$
(22)
the kinetic energy $`K_{\mathrm{II}}`$ at the outer boundary of the HII region is related to the radius $`R_{\mathrm{II}}`$ of the HII region.

The flow equation (17) can be solved analytically. The solution satisfying the boundary condition (22) is

$$\frac{GM_{}}{r}=K+2\sqrt{K_{\mathrm{crit}}}(\sqrt{K}\sqrt{K_{\mathrm{II}}})+2K_{\mathrm{crit}}\mathrm{ln}(\frac{\sqrt{K}\sqrt{K_{\mathrm{crit}}}}{\sqrt{K_{\mathrm{II}}}\sqrt{K_{\mathrm{crit}}}})(\mathrm{if}K_{\mathrm{II}}K_{\mathrm{crit}}),$$
(23)
and

$$K=K_{\mathrm{crit}}(\mathrm{if}K_{\mathrm{II}}=K_{\mathrm{crit}}).$$
(24)
In the limit of small radiation force (i.e., $`K_{\mathrm{crit}}0`$), only the first term on the right hand side of equation (23) remains and the flow becomes free-fall. The second and third terms are corrections due to the radiation force.

From the above characteristics of the flow equation (17), the growth of the HII region can be described as follows: As long as the the forming star is not massive enough to emit a large quantity of ionization photons, the HII region around it is compact. Since the radius $`R_{\mathrm{II}}`$ is small, $`K_{\mathrm{II}}>K_{\mathrm{crit}}`$ in this stage. The flow inside the HII region is hardly affected by the ionizing radiation force. As the star grows in mass and begins to emit a large quantity of ionizing photons, the radius of the HII region becomes larger and larger. When the radius of the HII region reaches the critical value

$$R_{\mathrm{IIcrit}}=\frac{2GM_{}}{(R_\mathrm{H}/t_{\mathrm{acc}})^2}=0.59\times 10^3R_{\mathrm{}}(M_{}/100M_{\mathrm{}})^3\dot{m}_{3}^{}{}_{}{}^{1},$$
(25)
a uniform velocity flow is attained inside the HII region. Further growth of the HII region would result in non-steady accretion and the accretion would eventually halt. We cannot discuss the later evolution of the accreting flow within the framework of steady accretion. According to hydrodynamical calculations for the accretion being halted by the radiation force onto dust, the accretion rate begins to decrease with some oscillations when the luminosity reaches the critical value (Yorke & Krügel 1977). We speculate that a similar phenomenon would occur also in our case.

We justify here the assumption that the I-front remains in R-type until it reaches the critical radius $`R_{\mathrm{IIcrit}}`$. When the radius of the HII region reaches $`R_{\mathrm{IIcrit}}`$, the flow velocity in the HII region is $`u_{\mathrm{crit}}`$ (see eq. 20). The I-front is R-type if

$$|u_{\mathrm{crit}}|>2c_{\mathrm{sII}},$$
(26)
where $`c_{\mathrm{sII}}`$ is the sound speed of the ionized matter. Then, our assumption of R-front is justified as long as

$$t_{\mathrm{acc}}<\frac{R_\mathrm{H}}{2c_{\mathrm{sII}}}2\times 10^6\mathrm{yr},$$
(27)
where for $`c_{\mathrm{sII}}`$ we used the value at 104K. The right hand side of equation (27) is as large as the lifetime of a massive star. Recall that the accretion time $`t_{\mathrm{acc}}`$ means the time needed to form a star of mass $`M_{}`$ at an accretion rate $`\dot{M}`$. Since the halting must occur within the stellar lifetime, the condition (27) is always met. Therefore, in our case, we can conclude that the I-front remains R-type until the halting of the accretion. If the mass accretion rate decreases significantly, the formally defined accretion timescale $`t_{\mathrm{acc}}`$ (eq.) can be larger than the stellar lifetime, and hence, the condition (27) would be violated within the stellar lifetime. However this decrease of the mass accretion rate is expected only in the very late stage of the mass accretion phase where the mass of the envelope is significantly depleted. Therefore the mass of the star eventually becomes on the order of the mass of the parental cloud. That is, the halting mechanism due to the ionizing photons considered here is not important in this case.

Next, we estimate the ionizing photon emissivity needed for attainment of the uniform velocity solution. Substituting the density distribution at the critical flow (21) into the equation of ionization equilibrium (5), we obtain the critical ionization photon emissivity

$`Q_{\mathrm{crit}}`$ $`=`$ $`{\displaystyle \frac{\alpha \dot{M}^2}{8\pi m_\mathrm{p}^2K_{\mathrm{crit}}}}({\displaystyle \frac{1}{R_{\mathrm{in}}}}{\displaystyle \frac{1}{R_{\mathrm{IIcrit}}}})`$ (28)
$`=`$ $`({\displaystyle \frac{h\nu _{\mathrm{ion}}}{c}})^1{\displaystyle \frac{GM_{}^2}{R_\mathrm{H}}}({\displaystyle \frac{1}{R_{\mathrm{in}}}}{\displaystyle \frac{1}{R_{\mathrm{IIcrit}}}}),`$ (29)
above in which the flow is decelerated in the HII region. We have neglected the accretion term $`\dot{M}/m_\mathrm{p}`$ as in §2.

In equation (28), if the right hand side is dominated by the first term (i.e., $`R_{\mathrm{in}}R_{\mathrm{IIcrit}}`$),

$$Q_{\mathrm{crit}}=(\frac{h\nu _{\mathrm{ion}}}{c})^1\frac{GM_{}^2}{R_{\mathrm{in}}R_\mathrm{H}}=0.64\times 10^{53}(R_{\mathrm{in}}/10R_{\mathrm{}})^1(M_{}/100M_{})^2.$$
(30)
As a typical value of $`R_{\mathrm{in}}`$, we take 10$`R_{\mathrm{}}`$, which is approximately the radius of ZAMS Pop III stars of 500$`M_{\mathrm{}}`$ (El Eid, Fricke, & Ober 1983).

Although we describe both $`Q_{\mathrm{crit}}`$ (eq. ) and $`Q_{\mathrm{crit}}^{(\mathrm{FF})}`$ in §2 (eq. ) as the critical ionizing photon emissivities, mechanisms of halting the accretion are different. In §2, the accretion is halted by the gas pressure in the HII region, whereas here what halts the accretion is the radiation pressure of ionizing photons. The inclusion of ionization radiation force results in enhancement of $`Q_{\mathrm{crit}}`$ because of an increase in the density (and hence, the recombination rate) due to the reduced flow velocity.

So far, we have only considered ionizing radiation as a source of radiation force. However, for massive stars the luminosity approaches the Eddington limit, and then the radiation force due to electron scattering becomes important. We can include this effect in our scheme without difficulty. If we use $`G_{\mathrm{eff}}=G(1\mathrm{\Gamma })`$ and $`R_{\mathrm{H},\mathrm{eff}}=R_\mathrm{H}/(1\mathrm{\Gamma })`$ (where $`\mathrm{\Gamma }=L/L_{\mathrm{Edd}}`$) instead of $`G`$ and $`R_\mathrm{H}`$, the flow equation including the radiation force due to electron scattering retains the same form as equations (15) and (17). Note that $`G`$ in (22) is not replaced since the gravity outside the HII region is not altered. After the same procedure as that leading to (28), the critical emissivity of ionization photons becomes

$$Q_{\mathrm{crit}}=(\frac{h\nu _{\mathrm{ion}}}{c})^1\frac{GM_{}^2}{R_\mathrm{H}}[\frac{(1\mathrm{\Gamma })^2}{R_{\mathrm{in}}}\frac{1}{R_{\mathrm{IIcrit}}}].$$
(31)
The critical ionizing photon emissivities $`Q_{\mathrm{crit}}`$ given by equations (28) and (31), namely the critical values with and without the effect of the electron scattering, are both depicted in Figure 2 as functions of the central stellar mass. In this figure, we have used $`\dot{M}=10^3M_{\mathrm{}}\mathrm{yr}^1`$ and $`R_{\mathrm{in}}=10R_{\mathrm{}}`$ as before. Also shown are the ionizing photon emissivity $`Q`$ of ZAMS Pop III stars and the critical value $`Q_{\mathrm{crit}}^{(\mathrm{FF})}`$ obtained by assuming the free-fall (eq. ). While $`Q`$ exceeds $`Q_{\mathrm{crit}}^{(\mathrm{FF})}`$ at about 300$`M_{\mathrm{}}`$ as already mentioned in §2, $`Q_{\mathrm{crit}}`$ is about two orders of magnitude larger than the actual value $`Q`$ even at $`10^3M_{\mathrm{}}`$ for our assumed value of $`R_{\mathrm{in}}=10R_{\mathrm{}}`$, even though the electron scattering reduces $`Q_{\mathrm{crit}}`$ considerably. Recalling that $`Q_{\mathrm{crit}}R_{\mathrm{in}}^1`$, we can conclude that $`Q_{\mathrm{crit}}>Q`$ for $`M_{}<10^3M_{\mathrm{}}`$ unless $`R_{\mathrm{in}}`$ becomes as large as $`10^3R_{\mathrm{}}`$ and $`Q`$ remains the same value. We consider that this situation does not occur in reality as long as we regard $`R_{\mathrm{in}}`$ as the stellar surface. Therefore, the growth of a HII region does not prevent the formation of such massive stars as $`10^3M_{\mathrm{}}`$.

## 4 Summary and Discussion

To estimate the mass upper limit of metal-free stars imposed by the formation of an HII region, we have found solutions for spherically symmetric and steady accretion flow onto a star emitting ionizing radiation.

The behavior of the flow is determined by the velocity of the flow entering the HII region. If the flow velocity at the edge of the HII region goes beyond (below) a critical value, the flow is always accelerating (decelerating, respectively) in the HII region. When this value is equal to the critical value, the velocity remains constant in the HII region. In this critical flow, the radiation force due to ionizing photons exactly balances the gravity. We applied those solutions for a compact HII region forming around an accreting star. As the HII region grows in radius, the accreting flow evolves from an accelerating solution to a critical one. Soon after the critical flow is reached, accretion is halted by radiation force due to the ionizing radiation. However, even stars as massive as $`10^3M_{\mathrm{}}`$ are unable to emit ionizing photons sufficient to halt the accretion. The halting by gas pressure is even more difficult. Therefore, contrary to the previous expectation, the formation of the HII region does not impose a stringent mass upper limit (at least up to $`10^3M_{\mathrm{}}`$) on metal-free stars.

We have estimated the upper limit of the mass of Pop III stars in relation to the formation of an HII region. Note that other effects (e.g., stellar wind, mass outflow due to pulsational instability, etc.) could grow in importance and might decrease the upper mass limit of Pop III stars below the value obtained in this paper. We leave detailed studies on these topics for future work.

Also, in this paper, we have adopted very simplified assumptions, for example spherical symmetry, steady accretion, etc. Here, we discuss other complexities and possible deviations from our picture.

We have considered only photoionization and electron scattering as sources of radiation force. Here, we mention briefly the radiation force due to Ly$`\alpha `$ photons. Ly$`\alpha `$ photons are emitted from recombination in the HII region. The flow inside the HII region is exerted by the radiation force due to Ly$`\alpha `$ photons. We can extend our theory easily to include the Ly$`\alpha `$ pressure , since, in the HII region, the radiation force due to the Ly$`\alpha `$ photons is proportional to that due to ionization photons (equation 12 of Haehnelt 1995; see also Braun & Dekel 1989; Bithell 1990). However, according to Haehnelt(1995), the radiation force due to Ly$`\alpha `$ photons is at most of marginal importance relative to that due to the ionizing photons. Thus, for the sake of simplicity, we have not included it. The HII region is surrounded by an HI layer, which is very optically thick to Ly$`\alpha `$ photons. Without dust grains, which absorb Ly$`\alpha `$ photons and reemit them as infrared photons, the Ly$`\alpha `$ photons must diffuse out through the HI layer. In the course of this, the HI layer is pushed outward by those photons. Doroshkevich & Kolesnik (1976) argued that this mechanism expels the HI layer soon after the HII region is formed and thereby limits the mass of stars below 10 M. However, Harrington (1973) showed that, even without dust, the two-photon emission process decreases the Ly$`\alpha `$ photon density drastically, and the Ly$`\alpha `$ radiation force is not dynamically important in an HI layer surrounding an HII region. In our case, supersonic motion in the accreting envelope also decreases the Ly$`\alpha `$ photon density. Considering these facts, it is likely that the Ly$`\alpha `$ radiation force does not play a significant role in our case. Therefore, we chose to neglect it here and assumed the free-fall outside the HII region.

We have identified the base of an HII region as the stellar surface, and have taken the typical value of the stellar radius $`10R_{\mathrm{}}`$ as the inner boundary radius $`R_{\mathrm{in}}`$ of the HII region. However, when the accretion rate is high, a photosphere may be formed in the accreting flow (Wolfire & Cassinelli 1986; Stahler, Palla, & Salpeter 1986). In this case, $`R_{\mathrm{in}}`$ should be taken as the radius of the photosphere, since we assume the HII region is optically thin to continuum absorption except for photoionization outside $`R_{\mathrm{in}}`$. In spite of higher $`R_{\mathrm{in}}`$, this effect generally works towards the smaller HII region, since $`Q`$ drops (Wolfire & Cassinelli 1986). Hence, our conclusion of higher upper mass limit than the former estimate remains the same.

The spherical symmetry is clearly an oversimplification if the accreting matter has large angular momentum and an accretion disk is formed. Even in this case, the spherical symmetric flow is a good approximation outside the centrifugal radius. Suppose here that inside the centrifugal radius, there is a “cavity” around the disk. In this case, $`R_{\mathrm{in}}`$ should be taken as the radius of the cavity. If $`Q`$ is unattenuated inside the cavity, the halting of the flow outside becomes easier because of higher $`R_{\mathrm{in}}`$ and then lower $`Q_{\mathrm{crit}}`$ (see eq. ). The accretion through the disk might continue, however (Nakano 1989). Those issues are still too speculative and beyond the scope of this paper.

Although the accelerating solution is Rayleigh-Taylor stable, it might be unstable to perturbations: if a portion of the flow becomes slightly slower than the rest, it becomes slower and slower relative to the average flow owing to the increased radiation force. In this case, density irregularities or blobs could be formed in the flow, and our spherically symmetric solution might be regarded as an approximate description of the average flow. Further study of the stability and dynamical evolution of the flow will be interesting.

We thank H.Sato for helpful comments and V.Bromm for providing the data of Pop III ZAMS. This work is supported in part by Research Fellowships of the Japan Society for the Promotion of Science for Young Scientists, grant 6819.


text:
# Acknowledgements

## Acknowledgements

We are grateful to Gary Shiu for helpful discussions. G.W.G., H.L. and C.N.P. are grateful to UPenn for support and hospitality during the course of this work.


text:
# Scanning Tunneling Microscopy of Defect States in the Semiconductor Bi2Se3

## Abstract

Scanning tunneling spectroscopy images of Bi2Se3 doped with excess Bi reveal electronic defect states with a striking shape resembling clover leaves. With a simple tight-binding model we show that the geometry of the defect states in Bi2Se3 can be directly related to the position of the originating impurities. Only the Bi defects at the Se sites five atomic layers below the surface are experimentally observed. We show that this effect can be explained by the interplay of defect and surface electronic structure.

Understanding the electronic properties of defects and the ability to control them will be crucial for the performance of the future microelectronic devices kane . Scanning Tunneling Microscopy (STM) represents a unique tool for the studies of defects as it combines atomic scale resolution with local spectroscopic capability. However, STM observation and analysis of defect states in semiconductors are complicated by surface effects such as in-gap surface states and reconstruction. These effects are avoided at the (110) surfaces of a number of III-V semiconducting systems gaasbands , attracting extensive research sigaas charge . A number of point defect types have been observed. However, positions of these defects with respect to the surface plane could be inferred only from indirect observations. The interpretation of such observations is complicated by the drastic effect the surface proximity may have on the defect states el2 .

Modeling STM measurements of defects in semiconductors is not straightforward: Approximation of the STM images by maps of the local surface electronic density of states tersoff  is justified only if the charge relaxation rates of defect states significantly exceed the tunneling rate of electrons between the tip and the sample theory . Tip-induced effects also need to be taken into account. These may include both local band bending sigaas , and charging of the defect states by the tunneling current, resulting in bias voltage-dependant lattice relaxation in the vicinity of the defect atoms charge . Careful analysis is necessary to clearly separate these effects from the intrinsic defect properties, and the bulk features of the observed defect states from the surface effects.

In this paper we present cryogenic STM and scanning tunneling spectroscopy (STS) studies of the layered narrow gap semiconductor Bi2Se3, which can be viewed as a model system for the STM study of near-surface defect states. The bonding scheme of Bi2Se3 allows a direct determination of the position of a subsurface defect atom with respect to the surface: Atomic planes consisting of either Bi or Se hexagonally arranged atoms are stacked in a close-packing fcc fashion; 5 atomic planes with atomic order Se1-Bi-Se2-Bi-Se1 (Fig. 1(a)) form a layer. The layers are weakly bound to each other by Se1-Se1 bonds. Both valence and conduction bands are formed almost exclusively by the $`4p`$ and $`6p`$ orbitals of Se and Bi respectively bisebands . For each atom, the closest neighbors from the adjacent atomic planes form almost a regular octahedron, so the bonding can be roughly approximated by strongly interacting $`pp\sigma `$ chains of atomic $`p`$-orbitals (Fig. 1(b)), with a weaker $`pp\pi `$-type interaction between adjacent chains. A substitutional defect is therefore likely to produce a perturbation in the electronic local density of states (DOS) predominantly along the three $`pp\sigma `$ chains passing through the defect atom. Hence, the defect state should be observed at the surface as three spots of modified DOS around the atoms terminating these chains at the surface.

We have performed STM and STS measurements using a custom built low-temperature microscope with direct immersion in liquid He-4 tessmer . Stoichiometric Bi2Se3 single crystal samples, as well samples doped with 2–5% excess Bi or Se were grown by a directional solidification technique. The stoichiometric as-grown Bi2Se3 samples are n-type with carrier concentration of about 1019 cm-3. Doping samples with excess Bi introduces substitutional Bi defects at the Se sites (BiSe antisites), which are shallow acceptors. However, because of the low solubility of Bi in Bi2Se3 lovett , the Bi-doped samples are n-type due to the compensating defects. Doping samples with excess Se introduces shallow donor-type Se substitutional defects at the Bi sites (SeBi antisites). The samples are sufficiently inert to obtain atomic resolution in air. However, to minimize the surface contamination, in the experiments reported here, the samples were cleaved in situ or in a glove box directly attached to the STM setup in ultrapure He gas prior to transfer to the STM with subsequent cooling to T=4.2 K.

To map out the defect states, we have performed differential current imaging tunneling spectroscopy (CITS) measurements strosciobook . The CITS maps were acquired by fixing the tip at each point during topographic imaging and measuring the differential conductance at various bias voltages with the feedback loop disabled. Only Bi-doped samples exhibited an appreciable density of observable defect states, as shown in Fig. 2(a) for a sample doped with 5% excess Bi (Bi2Se2.85). The map was acquired with sample bias Vb=$`0.2`$ V in the feedback mode (with tunneling current set to $`50`$ pA), and Vb=$`0.45`$ V for the conductivity measurement. The defects appear as regular clover-shaped bright features, indicating areas with locally enhanced conductance at Vb=$`0.45`$ V. Irregular spots in this image (mostly in the upper left corner) resulted from topographic defects. The data discussed below were obtained on a more weakly doped Bi2Se2.95 sample, where the defect density was reduced.

Fig. 2(b),(c) present topographic maps of the sample area encompassing an isolated clover-shaped defect state. Topographic image Fig. 2(b), acquired at a sample bias voltage Vb=$`0.3`$ V, shows a periodic atomic structure, indicating no significant structural variation associated with the defect. The height of the atomic corrugations in Fig. 2(b) is about $`30`$ pm. Fig. 2(c) shows a topographic image of the same area acquired at Vb=$`0.6`$ V, where the largest corrugations, locally enhanced by tunneling through the defect state, are about $`100`$ pm high. The highest amplitude of the defect state correlates with positions of three surface Se atoms (marked with larger black circles in schematic Fig. 2(d)), forming a regular triangle. These atoms terminate three $`pp\sigma `$-bonded chains passing through the Se1 site five atomic layers below the surface, for which we also use notation Se3. Since the observed defects appear only in Bi-doped samples and they originate from Se sites, we attribute them to the BiSe3 antisites.

Fig. 3(a) presents a series of differential conductance spectra acquired in the vicinity of an isolated defect. The spectra were obtained by numerical differentiation of 60 I-V curves with setpoint parameters Vb=$`0.3`$ V, I=$`0.8`$ nA. At the measurement temperature of 4.2 K, thermal broadening is negligible on the displayed bias voltage scale. The defect state appears as a broad resonance (indicated by an arrow) in the energy range where the differential conductance is suppressed away from the defect. The established semiconducting gap value is about $`0.3`$ eV bisebands lovett , therefore the defect levels appear *inside* the valence band. Theoretical modeling is necessary to understand this spectroscopic feature, as well as why only BiSe antisites five atomic planes below the surface are observed.

First, we performed ab initio calculations in the full potential relativistic LAPW formalism lapw  within LDA approximation. To model the surface, a supercell geometry was used, with distance between slabs (consisting of 15 atomic planes, or 3 layers, each) increased by $`0.5`$ to $`1.5`$ nm as compared to the bulk crystal structure. The calculated band structure did not exhibit significant variation for the slab separation larger than $`0.4`$$`0.5`$ nm, therefore we found the slab separation of $`0.7`$ nm, used for the calculations presented below, sufficient for modeling the surface. The differential conductance spectra were approximated by the local DOS in the center of the gap between the slabs tersoff . The calculation presented in Fig. 3(b) was performed for a position above a Se1 atom, although we found the variation of the calculated spectra with position respective to the surface atoms to be insignificant. The calculation reproduces both the finite conductance in the bulk semiconducting gap, and the suppressed conductance just below the gap. In Fig. 3(c), the total DOS calculated in the slab geometry is compared to the calculation of the bulk DOS, which reproduces the accepted semiconducting gap value of $`0.3`$ eV LDAcomment . Band structure analysis indicates that the highest valence band (HVB) states are predominantly Se1-Se1 antibonding type. As the Se1-Se1 bonds are broken at the surface, the splitting of these states is reduced, resulting in the observed suppression of the differential conductance in the HVB energy range (Fig. 3(b)). The states that appear in the gap have high dispersion along the surface. They originate from the rehybridization of the surface Se1 valence with Bi conduction states, bringing the latter down below the bulk conduction band minimum.

ab initio calculations of a single defect state are complicated by the large cluster or supercell size necessary to model the impurity states without introducing artificial interaction between defects. Instead we use a LCAO approximation as a simple model of the system bisebands . This model presents just a qualitative argument and is not capable of reproducing the detailed electronic structure or the semiconducting gap value. However, it gives a surprisingly good qualitative agreement with the experimental observations and first principles calculations. In a tight binding formalism

$$H\psi =H\underset{i}{}u_i\varphi _i=\underset{i}{}E_iu_i\varphi _i+\underset{ij}{}V_{ij}u_i\varphi _j,$$
(1)
where $`\varphi _i`$ are atomic wave functions and $`V_{ij}`$ are off-diagonal matrix elements of $`H`$. We approximate $`E_i`$ by the atomic term values harrison  of Bi and Se, and take into account only $`pp\sigma `$ interaction between the closest neighbors, as shown in Fig. 1(b). Thus the problem is reduced to a system of noninteracting one-dimensional chains, with three matrix elements $`V_1`$,$`V_2`$, and $`V_3`$, corresponding to Se1-Bi, Se2-Bi, and Se1-Se1 $`pp\sigma `$ bonds. Consider first a 5-atom chain Se1-Bi-Se2-Bi-Se3 (we also call it a unit), representing a single layer in our approximation. Here Se3 position is equivalent to Se1. The highest energy filled state is nonbonding

$$\psi _0=\frac{1}{\sqrt{2+(V_1/V_2)^2}}(\varphi _{Se1}+\varphi _{Se3}V_1/V_2\varphi _{Se2}).$$
(2)
We model the bulk by a long chain of Se1-Bi-Se2-Bi-Se1 units, with the interaction between layers expressed by $`V_3`$. Due to the Se1-Se1 interaction, the nonbonding level is split, and the HVB states become antibonding in the sense of Se1-Se1 bond character. The splitting is large, because the state Eq.(2) has a significant weight on the interlayer Se1 atoms. As the interlayer bonds are broken at the surface, the HVB states of a chain decay at the surface (the end of the chain), as shown in Fig. 4(a), where the antibonding character of the valence band maximum (VBM) state can also be seen. In semiclassical terms, as illustrated in Fig. 4(b), the surface gap is larger than the bulk value. This effect is in agreement with the more accurate first principles calculations and spectroscopic measurements, Fig. 3(b), (c).

To highlight the importance of these surface effects for the observation of the defect states, in Fig. 4(c) we plot the calculated dependence of the near-gap energy levels on the position of the BiSe antisite. The second layer (positions 4–6) is only weakly affected by the proximity of the surface, so the defect level is split from the VBM and its energy is only weakly dependant on the position. As the antisite position approaches the surface (positions 1–3), the defect level energy is reduced as the surface gap opens up, so that the defect state merges with the bulk valence states, forming a resonance *inside* the bulk valence band. This behavior is supported by the semiclassical picture shown in Fig. 4(b). Only the BiSe3 state is observed in the experiment, because, as our model suggests, BiSe1 and BiSe2 states are so much lowered in energy by the proximity of the surface that they form small amplitude broad resonances in the valence band. STM images do not exhibit defect features associated with BiSe antisite in the second layer. Our model suggests, that they form bound states in the bulk gap, which cannot sustain STM current. Surface effects thus provide a mechanism for the charge relaxation of near-surface defect states through the bulk valence band.

Similar surface effects should be observable in other semiconductors, e.g. at (110) surface of GaAs, where the valence band is, like in Bi2Se3, suppressed at the surface gaasbands . As a result, in-gap impurity states may become resonances in the valence band, if the originating impurities are sufficiently close to the surface el2 . This suggests an alternative explanation for some of the published results sigaas2 . It may also be possible to induce the resonant behavior of near-surface defects by careful control of the surface band bending with doping and/or external field.

Resonances induced by near-surface defects can be contrasted to the bulk-like in-gap states. The origin of the spectroscopic broadening of the latter feenstra , and the mechanisms of their charge relaxation, allowing their observation with STM, need further theoretical and experimental studies theory . Variable temperature studies of the influence of the local defect distribution on the spectroscopic features of defect states may provide insight into these issues, and Bi2Se3 represents a convenient model system for such studies.

In summary, we have observed clover-shaped defect states in Bi2Se3 doped with excess Bi, which appear as resonances in the high valence band, and can be attributed to BiSe antisites in fifth atomic layer from the surface. In the analysis of these defect states, we have demonstrated the importance of the surface effects for the electronic structure of near-surface defects. While BiSe defects in the bulk Bi2Se3 form shallow acceptor levels, the near-surface defects produce resonances in the energy range of the bulk valence states, which are suppressed at the surface. We suggest that similar surface effects are likely to be observable in other semiconducting systems.

We thank J. Nogami, Norman O. Birge, M. I. Dykman, T. Hogan, and S. Lal for helpful discussions. This work was supported in part by NSF (DMR-0075230) and ONR/DARPA (N00014-01-1-0728). SHT acknowledges support of the Alfred P. Sloan Foundation.