text:
# The porous atmosphere of 𝜂-Carinae

## 1 Introduction

$`\eta `$-Carinae is probably one of the most remarkable stellar object to have ever been documented. About 150 years ago, the star began a 20 year long giant eruption during which it radiated a supernova-like energy of roughly $`3\times 10^{49}ergs`$ (Davidson & Humphreys 1997). Throughout the eruption it also shed some $`12_{}`$ of material carrying approximately $`6\times 10^{48}ergs`$ as kinetic energy (Davidson & Humphreys 1997), while expanding at a velocity of $`650km/sec`$ (Hiller & Allen 1992, Currie et al. 1996). $`\eta `$-Carinae can therefore serve as a good laboratory for the study of atmospheres at extreme luminosity conditions.

At first glance, it appears that the star shed a large amount of material. Indeed, the inferred mass loss rate during the great eruption of $`0.1_{}/yr`$ is significantly larger than the mass loss rate inferred for the star today ($`\begin{array}{c}<\\ \end{array}10^3_{}/yr`$, Davidson & Humphreys 1997 and references therein). However, considering that the luminosity during the great eruption is estimated to be significantly above the Eddington limit, we shall show that the star should have had a much higher mass loss rate. In fact, it should have lost during the 20 year eruption more mass than its total mass, giving rise to an obvious discrepancy.

A review of our current knowledge of $`\eta `$ Car can be found in Davidson & Humphreys (1997). In section 2 we summarize how a wind solution for the star $`\eta `$ Car should be constructed. Since the luminosity is very high, the effects of convection must be taken into account. In section 3 we integrate the wind equations to show that no consistent solution for $`\eta `$ Car exists within the possible range of observed parameters. Section 4 is devoted to possible classical solutions to the discrepancy, showing that no such possibility exists. In section 5, we show that a porous atmosphere is a simple and viable solution to the wind discrepancy.

## 2 Solving for the Wind

Since the mass of $`\eta `$-Car is estimated to be of order $`100120_{}`$ (Davidson & Humphreys 1997), the average luminosity in the great eruption was clearly super-Eddington (of the order of 5 times the Eddington limit). That is to say, the radiative force upwards, assuming the smallest possible opacity (for ionized matter) given by Thomson scattering, was significantly larger than the gravitational pull downwards. Optically thin winds formally diverge at the Eddington limit (e.g., Kudritzki et al. 1989 and references therein). Consequently, a consistent wind solution requires an optically thick wind. We thus look for a wind in which the sonic point (which is the point at which the local speed of the outflow equals the speed of sound) is below the photosphere. Moreover, since the duration of the eruption is longer than the sound crossing time of the star by about a factor of 50, a stationary wind appears to be a good approximation.

In practically all super-sonic wind theories which describe super-sonic outflows from an object at rest, a consistent stationary solution is obtained only when the net driving force of the wind (excluding the pressure gradient) vanishes at the sonic point111 The exception is line driven winds in which the force is explicitly a function of $`dv/dr`$ which is actually an approximation to the line transfer equations. If we had written the proper radiation transfer equations for this case which only implicitly depend on $`dv/dr`$, we would have recovered that the sonic point coincides with the point at which the total force vanishes (cf Mihalas & Weibel Mihalas 1984 §107). Moreover, line driven winds are important only under optically thin conditions while we describe the optically thick part of the wind. . Thus, material experiencing a super-Eddington flux necessarily has to be above the sonic point. If most of the envelope carries a super Eddington flux, then no consistent stationary wind solution can be obtained and in fact, the object will evaporate on a dynamical time scale. In most systems however, this need not be the case. For example, in very hot systems (e.g., hot neutron stars during strong X-ray bursts, Quinn & Paczynski 1985), the opacity in the deep layers is reduced due the reduced Klein-Nishina opacity for Compton scattering at high temperatures. Thus, the sonic point in these objects is found where the temperature is high enough to reduce the opacity to the point where the flux corresponds to the local Eddington limit.

Another important effect, which should be taken into account, is convection. Deep inside the atmosphere, convection can carry a significant part (or almost all) of the energy flux, thus reducing the radiative pressure to a sub-Eddington value. In fact, as the radiative flux approaches the Eddington limit, convection generally arises and carries the lion share of the total energy flux (if it can) to keep the system at a sub-Eddington level (Joss et al. 1973). Although the total flux in the entire envelope (or almost all of it) can be equivalent to a super-Eddington flux, up to some depth below the photosphere, convection carries most of the flux so as to reduce the radiative flux alone into a sub-Eddington value. A consistent wind solution should therefore, have its sonic point at the location where the most efficient convection cannot carry enough flux any more. As we shall soon see, the problem in $`\eta `$-Car is that this point is relatively deep within the atmosphere, where the density is so high that the expected mass loss is significantly higher than the observed one.

To see this in a robust way we integrated numerically the wind equations starting from the photosphere inwards. The equations are those that describe optically thick spherically symmetric winds (Quinn & Paczynski 1985; Żytkow 1972; Kato & Hachisu 1994). The equations of mass conservation, momentum conservation and temperature gradient are

$$4\pi r^2\rho v=\dot{M}=\mathrm{const}$$
(1)
$$v\frac{dv}{dr}+\frac{GM}{r^2}+\frac{1}{\rho }\frac{dP_g}{dr}\frac{\chi _r}{4\pi r^2c}=0$$
(2)
$$\frac{dT}{dr}=\frac{3\chi \rho _r}{16\pi r^2caT^3}(1+\frac{2}{3\chi \rho r}),$$
(3)
with standard notation. The parenthesized term in the last equation is a simple approximate interpolation that has the correct asymptotic limits for optical depths much larger and much smaller than unity (Quinn & Paczynski 1985). The last equation is the integrated form of the energy conservation equation. Unlike the aforementioned references, we specifically include advection by a maximally efficient convection. Thus, the integrated form of the energy conservation equation becomes

$`_r+\dot{M}\left({\displaystyle \frac{v^2}{2}}+w{\displaystyle \frac{GM}{r}}\right)+_{conv}`$ $``$ (4)
$`_r{\displaystyle \frac{GM\dot{M}}{r}}+_{adv}+_{conv}`$ $`=`$ $`\mathrm{\Lambda }_{tot}=_{obs}+_{kin,\mathrm{}},`$ (5)
were $`\mathrm{\Lambda }_{tot}`$, $`\dot{M}`$, $`_{obs}`$ and $`_{kin,\mathrm{}}`$ are the total energy output of the star, the wind mass loss rate, the observed luminosity at infinity and the kinetic energy flux at infinity. On the other hand, $`_r,v,w,_{conv},_{adv}`$ are respectively, the local radiative luminosity, velocity, enthalpy, convective flux and advected flux (as internal and kinetic energies). The expression adopted for $`_{conv}`$ is $`4\pi r^2uv_s`$ where $`u`$ is the internal energy per unit volume and $`v_s`$ is the speed of sound. By no means can convection be more efficient than this expression since highly dissipative shocks are unavoidable at higher speeds. It is likely that the maximally efficient convection is somewhat less efficient than this expression, but this will only aggravate the problem that we shall soon expose. Detailed calculations of the wind were carried out. The calculations include the latest version of the OPAL opacities (Iglesias & Rogers 1996). It is found that the total opacity below the photosphere has comparable contributions from Thomson scattering and absorption processes. This implies that the modified Eddington limit, in which the Thomson opacity is replaced by the local total opacity, is somewhat lower than the classical Eddington limit.

Since $`T_{\mathrm{𝑒𝑓𝑓}}9000^{}`$K222This is the typical observed effective temperature for LBV’s during outbursts (Humphreys & Davidson 1994). If the temperature is higher than this value, the inferred bolometric magnitude of $`\eta `$ Car during the eruption would be more negative, increasing the Eddington factor. If the temperature is lower than $`7000^{}`$K, the opacity at the photosphere and outwards rises abruptly (Davidson 1987), thus reducing the modified Eddington limit. In both cases, the discrepancy will be aggravated., the average luminosity implies a photospheric radius of $`10^{14}cm`$. Note that since it is a thick wind, the exact definition of the photosphere is ambiguous. Nevertheless, different definitions do not change the results by more than $`1020\%`$. Knowing that the observed mass loss rate is roughly $`0.1_{}/yr`$ (which gives the observed $`2_{}`$ of shed material in 20 years, Davidson & Humphreys 1997), a specified flow speed at the photosphere can be translated to a required density. We can therefore integrate the wind equations inwards. If a consistent wind solution can be obtained for some value of the imposed velocity in the photosphere (which has to be between $`v_s`$ and $`v_{\mathrm{}}`$), then the integration inwards should reach a sonic point at which the total force on the gas vanishes. This will be attained if the convective and advective fluxes can carry a significant amount of the total flux so as to reduce the residual radiative flux to a sub-Eddington one.

## 3 The Discrepancy

We define the luminosity needed to be carried by convection and advection in order to bring about the vanishing of the total local force as $`_{crit}`$. If enough energy is advected and convected then $`_r`$ will be reduced to the local modified Eddington flux:

$$_{Edd,mod}=\frac{4\pi cGM}{\chi }$$
(6)
with $`\chi `$ the local opacity which can be larger than the Thomson opacity. Thus, from eq. (5), the critical advective+convective flux can be written as

$$_{crit}=\mathrm{\Lambda }_{tot}+\frac{GM\dot{M}}{r}_{Edd,mod}.$$
(7)
Figure 1 shows the fraction $`\eta (_{adv}+_{conv})/_{crit}`$ at the sonic point. A consistent solution can be found only if $`\eta =1`$ at the sonic point.

Inspection of the figure clearly shows that the space of possible observed values does not contain a viable and consistent solution. This is of course irrespective of whether a solution from the photosphere outward can or cannot be obtained. The discrepancy arises because a wind corresponding to the observed low mass loss rate necessarily has a sonic point that is not deep enough to have either convection or advection as an efficient mean of transporting energy. This can be seen from the optical depth at which the sonic point is obtained. In all cases, $`1\begin{array}{c}<\\ \end{array}\tau <300`$. However, convection is efficient only up to an optical depth of $`\tau c/v_s300`$ for $`p_{rad}p_{gas}`$ (Shaviv 2000b), or even deeper for larger radiation pressures (i.e., when close to the Eddington limit).

## 4 Unfeasible Solutions to the Discrepancy

Can the discrepancy be resolve with a classical assumption? Since the discrepancy is rather large, assuming the wind to emerge from an angular fraction $`f`$ from the star does not relax the problem (it actually aggravates the problem because more material will be blown away from the higher luminosity regions). Another possibility that fails is having a higher velocity in the photosphere than the one observed today for the shed material. This might be the case if the wind collides with previously ejected slow moving material. Even if such material did exist, the necessarily reduced mass loss rate inferred from the present day observed momentum aggravates the problem.

The problem is not mitigated if we relax the assumption that the mass loss rate and the luminosity are assumed to be constant in time throughout the eruption.

If one wishes to solve the problem using magnetic fields, then a solution can be found only if the magnetic energy density at the photosphere is significantly larger (by several orders of magnitude) then the equipartition value with the gas pressure. This of course seems unlikely.

Another option is to have the distance estimate to $`\eta `$-Car be three times smaller than $`2300pc`$. A shorter distance will remove $`\eta `$-Car out of the cluster Tr16 of massive stars inside which it is observed and leave it instead roaming the inter galactic arm space. Considering the short lifetime of the star, just a few million years, this possibility appears as very unlikely.

The problem can be solved if the mass of the star corresponds to a sub-Eddington luminosity. This proposed solution requires $`\eta `$-Car to be at least a $`1000_{}`$ star. However, this suggestion is at variance with much lower estimates (see for instance Davidson & Humphreys 1997 and references therein). Nevertheless, having such a massive star is in fact not completely unrealistic and would have far reaching consequences if found to be true.

## 5 A Viable Solution to the Discrepancy

As the title suggests, there is a clear solution to the discrepancy. As the results show, the sonic point appears to be between the optical depths of $`1`$ and $`300`$. The exact value cannot be obtained since it requires the integration outward from the photosphere, which owing to the relatively inaccurate effective temperature and therefore opacity, yields a wide range of results. If the mean radiative force between the point $`\eta =1`$ and the above found optical depth, is smaller than classically estimated, then a solution to the discrepancy can be found. Such a reduction in the mean radiative force is a natural result if the atmosphere is inhomogeneous.

Shaviv (1998) has shown that in an inhomogeneous atmosphere, the effective opacity used to calculate the average force is reduced relative to the effective opacity used for the radiation transfer in a homogeneous medium. The effective opacity used for the average force should be a volume flux weighted average of the opacity per unit volume333When the flux is frequency dependant, a similar average should be taken in order to find the radiative force. However, one then takes a flux weighted mean over frequency space. $`\chi _v\chi \rho `$. Namely,

$$\chi _{\mathrm{𝑒𝑓𝑓}}=\frac{\chi \rho F}{F\rho }.$$
(8)
The effect is universal and arises in inhomogeneous systems that conduct heat or electricity. Extensive discussions exist in the literature under a different terminology (Isichenko 1992). The only requirement is therefore, that close to the Eddington limit the star develop inhomogeneities. The transformation from an homogeneous to an inhomogeneous atmosphere at luminosities close to but below Eddington luminosities, was recently found to take place generically even in Thomson scattering atmospheres (Shaviv 1999; Spiegel & Tao 1999; Shaviv 2000a).

It was found that two different types of instabilities arise naturally when the luminosity approaches the Eddington limit (Shaviv 2000a). One instability is of a phase transition into a stationary nonlinear pattern of “fingers” that facilitate the escape of the radiation. The second type of instability allows the growth of a propagating wave, from which one expects a propagating nonlinear pattern to form. The two possibilities are summarized in figure 2. Both instabilities bring about a reduction of the average radiative force on the matter and a significant reduction of the mass loss rate since the sonic surface can sit near (or not much below) the photosphere. In both cases, the nonlinear pattern is necessarily expected to form in the region between the radius $`r_{conv}`$ at which $`\eta =1`$, or in other words, that $`_{conv}+_{adv}`$ is large enough to have $`_r\begin{array}{c}<\\ \end{array}_{Edd,mod}`$, and the photosphere. When the pattern is stationary, the rarefied regions have a larger than Eddington flux and the sonic surface in these regions is near $`r_{conv}`$. On the other hand, if the pattern is propagating, the flux may be larger locally than the Eddington limit but the time average of the force on a mass element is less than Eddington. Since the instability does not occur above the photosphere, it should be homogeneous and hence super-Eddington with a super-sonic flow.

Further analysis of the instabilities is needed to know which instability will dominate though it is more likely to be the phase transition since it is dynamically more important.

## 6 Summary

To summarize, the super-Eddington luminosity emitted by $`\eta `$-Car should have generated a much thicker wind with a sonic point placed significantly deeper than what can be directly inferred from the observations. A solution which lives in harmony with observations and theoretical modeling is a porous atmosphere, which allows more radiation to escape while exerting a smaller average force. It also means that the Eddington limit is not as destructive as one would a priori think it must be, even in a globally spherically symmetric case. Namely, all astrophysical analyses that employ the Eddington limit as a strict limit should be reconsidered carefully, even if they involve only unmagnetized Thomson scattering material. If $`\eta `$-Carinae could have been super-Eddington for such a long duration without “evaporating”, other systems can display a similar behavior.


format: text

text:
# Pulsation in two Herbig Ae stars: HD 35929 and V351 OriBased on observations carried out at the European Southern Observatory, La Silla, Chile under proposals number 62-I-0533, 63-I-0053

## 1 Introduction

Young pre–main-sequence (PMS) stars of mass greater than $``$1.5 M cross the region of pulsational instability during their contraction toward the main-sequence. The location of the instability strip in the H-R diagram has been recently identified by means of nonlinear models for the first three radial modes (Marconi & Palla 1998). The time spent by intermediate-mass stars within the boundaries of the strip represents a small fraction of the Kelvin-Helmoltz time scale, varying between $``$10% for a star of 1.5 M and $``$5% for a 4 M star. Despite the brevity of this phase, however, a number of known Herbig Ae stars have the appropriate combination of luminosity and effective temperature to become pulsationally unstable. In our previous study, we suggested to look for $`\delta `$ Scuti-type photometric variations with periods of minutes to several hours and amplitudes less than few tenths of magnitudes in a sample of Herbig Ae stars whose position in the H-R diagram coincides with the instability strip.
The identification of a few PMS objects pulsating as $`\delta `$ Scuti stars (Breger 1972), the prototype being the star HR 5999 (Kurtz & Marang 1995), has provided some support to the connection between variability and stellar pulsation. Of particular interest is the Ae star HD 104237 that shows both short- and long-term velocity changes of spectral lines (Donati et al. 1997). These variations indicate that the star is undergoing radial pulsations with a period of approximately 40 minutes and an amplitude of about 1 km s-1 (see also Böhm et al. 2000, in preparation). Interestingly, HD 104237 is the first intermediate mass PMS star with a measured magnetic field (Donati et al. 1997). A few new PMS candidate pulsators have been recently identified by Pigulski et al. (2000).
Stimulated by these initial results, we have started a photometric investigation of a sample of seven Herbig Ae stars with spectral types in the range A5 to F5, located within or near the boundaries of the instability strip. For some of them, large time scale variations have been observed during the long term monitoring program of variable stars conducted at ESO (LTPV project: Sterken et al. 1995 and references therein). However, no information is available on their variability on time scales shorter than 2 or 3 days.
The main goal of our study is to detect and characterize the pulsation properties of young stars. This way, we can improve our knowledge of their internal structure and obtain unique constraints on the theoretical predictions of the models. Ultimately, the analysis of the pulsation characteristics can yield an indirect estimate of the stellar mass. This represents a powerful method for stars that are not part of the restricted group of spectroscopic binary systems. In this Letter, we report the discovery of two additional Herbig Ae stars which show evidence for variability of the $`\delta `$ Scuti-type.

## 2 Selection of the sample, observations and data reduction

The selection of the stars was based on their spectral type. As an initial choice, we have adopted values published in the literature. To have an independent check on the effective temperature of the selected stars, we have used the Strömgren photometry provided by the ESO catalogues of the LTPV project. Then, we have placed the stars in the dereddened color-color \[m1\]–\[c1\] diagram (see Strömgren 1966), and compared their position with the theoretical colors from model atmospheres (Kurucz 1992). As a result, we derived a sample of seven stars with spectral types between F5 and A5: V346 Ori, V351 Ori, BF Ori, BN Ori, NX Pup, HD 35929, and AK Sco.

The photometric observations were performed in two runs with different telescopes at ESO, La Silla. In the first period (Dec. 23–29, 1998), we used the 0.5k$`\times `$0.5k CCD (Tektronix TK512CB) attached to the 0.9m DUTCH telescope. Five of the seven nights were photometric. In the second period (June 30 to July 4, 1999) we used the 2k$`\times `$2k CCD (LORAL/LESSER) attached to the 1.5 m. DANISH telescope. Unfortunately, weather conditions were excellent for accurate photometric observations only during one night. Due to their brightness, all of the sources are suitable for observations in the Strömgren $`uvby`$-system. Two comparisons stars were used for each candidate variable; whenever possible, we observed the same comparison stars as in the LTPV program. Data reduction and analysis to derive instrumental aperture magnitudes were performed following the standard procedures under the MIDAS environment. The typical intrinsic instrumental photometric error in each Strömgren filter is of the order of a few thousands of a magnitude.

## 3 Results: $`\delta `$ Scuti candidates

The stars V346 Ori, V351 Ori, BF Ori, HD 35929, and AK Sco have been found to show photometric variability. For V351 Ori and HD 35929, we could derive reliable pulsational periods that fall in the expected range of $`\delta `$ Scuti type variability111As for V346 Ori and AK Sco, the variation seems $`\delta `$ Scuti-like, with $`u`$ amplitudes of about 0.035 and 0.08 mag respectively but no period could be derived due to the poor phase coverage. As for BF Ori, its variation is clear ($`\mathrm{\Delta }u`$ $``$ 0.11 mag), but seems rather monotonic. In Table 1 we report their main properties.
HD 35929 was observed during two nights (Dec. 28 and 29, 1998). The stars HD 37399 and HD 37210 were used for comparison. Unfortunately, HD 37210 turned out to be variable with a maximum amplitude in the $`u`$-filter of $``$0.03 mag and a rather monotonic variation. Therefore, the differential light curve of HD 35929 was derived using the data points relative to HD 37399 only. The instrumental phased light curves in $`\mathrm{\Delta }u`$ and $`\mathrm{\Delta }b`$ are shown in the two upper panels of Fig. 1. The frequency spectrum of HD 35929 obtained from the $`u`$-filter is shown in the lower panel of Fig. 1. The frequency spectrum peaks at 5.1 d-1, i.e. 0.196d. We note that no significant differences in the calculated period are found using light curves in other filters.

V351 Ori was observed during the nights of Dec. 23 and 25, 1998. HD 38155 and HD 35298 were used as comparison stars. The latter turned out to be a variable star, with a maximum amplitude in $`u`$ of $``$0.04 mag and a quasi-periodic behaviour. Thus, we use only HD 38155 for comparison. The differential magnitudes in each filter ($`\mathrm{\Delta }u,\mathrm{\Delta }v,\mathrm{\Delta }b,\mathrm{\Delta }y`$, computed as variable$``$comparison) are shown in the upper panels of Fig. 2. The light curve in the $`u`$-band, where the largest amplitude is observed, was used to determine the observed period of V351 Ori. The frequency spectrum of V351 Ori obtained from the $`u`$-filter is shown in the lower panel of Fig. 2. The main maximum in this spectrum corresponds to a period of 0.058 d. Again, no significant differences in the derived periods occur when using light curves in other filters.

## 4 Evolutionary state of HD 35929 and V351 Ori

Although more systematic observations are needed to define the precise pulsational periods of HD 35929 and V351 Ori, the available data already provide some preliminary constraints on their position in the H-R diagram, stellar mass and evolutionary state. As remarked by the referee, the discussion presented in this section is useful as an illustration of how well-defined periods for these stars could indeed tightly constrain their evolutionary states.
Fig. 3 shows the location of HD 35929 in the H-R diagram. The dotted box accounts for the uncertainty in the spectral type (A5 to F0: Malfait et al. 1998, Miroshnichenko et al. 1997, van den Ancker et al. 1998) and distance ($`d=`$360 to 430 pc: van den Ancker et al. 1998). The instability strip for the first three radial modes, as predicted by Marconi & Palla (1998), is also shown together with the PMS evolutionary tracks computed by Palla & Stahler (1993) and the post-MS evolutionary tracks for 2.5 and 3.0 $`M_{}`$ of Castellani et al. (1999). The two circles indicate the best combination of the stellar parameters ($`M`$, $`L`$, $`T_{\mathrm{eff}}`$) that yield a period equal to the observed one, $`P=0.196\pm 0.005`$ d. These values are listed in Table 2. The two solutions indicate a mass of 3.4 or 3.8 M, pulsating in the first overtone (FO) and second overtone (SO) respectively: in both cases, HD 35929 can be considered a PMS pulsator, as expected.

A combination of parameters for a post-MS stellar mass can also reproduce the pulsational period observed in HD 35929. In this case, the best choice would be a post-MS model of a 2.7 $`M_{}`$ star with $`L=83L_{}`$ and $`T_{\mathrm{eff}}=6900`$ K, pulsating in the SO mode. The location of this solution is quite close to the lower circle shown in Fig. 3. Although only few specific studies exist on HD 35929, some evidence supports the fact that HD 35929 is a young star associated with the Ori OB-1c association. For example, Malfait et al. (1998) have discussed the infrared excess observed toward this star, whereas Miroshnichenko et al. (1997) have suggested that the star might be in a transition phase between a PMS Herbig Ae star and a $`\beta `$ Pictoris-type object. We also note that the pulsational character (period and H-R diagram location) of HD 35929 is similar in many respects to that of the well known Herbig Ae star HR 5999. From these considerations, we favor the conclusion that HD 35929 is a PMS star with a mass in the narrow range 3.4-3.8 $`M_{}`$, pulsating in the FO or SO.
Finally we note that if the Strömgren indices \[c1\] and \[m1\] measured for this star are taken into account, together with published values for the $`H\beta `$ index, one derives an effective temperature $`\mathrm{log}T_{\mathrm{eff}}=3.84`$ and a luminosity varying between 26 and 36 $`L_{}`$. This means that according to the Strömgren photometry, HD35929 should be located on a $`2.5M_{}`$ evolutionary track and the period based on present data would be too long even for a pulsation in the fundamental mode. A possible explanation for this inconsistency could be that HD35929 is also a rapid rotator ($``$150 Km s-1) so that the assumption of radial pulsation could be not completely correct. Future observations and numerical simulations are needed in order to properly address this problem.

Fig. 4 shows the location of V351 Ori in the H-R diagram. Here, the uncertainty on the luminosity is larger than for HD 35929, because of the distance ambiguity. The lower value corresponds to the minimum distance of 260 pc given in the Hipparcos catalogue (van den Ancker et al. 1998). The upper limit (inverted triangle) assumes that V351 Ori is located in the Orion molecular cloud at a distance of 460 pc. The dashed box corresponds to an uncertainty of $`\pm 0.01`$ dex in $`\mathrm{log}T_{\mathrm{eff}}`$. Finally, the filled circle marks the position estimated by van den Ancker et al. (1996) with the associated error bar. As already pointed out, present data for V351 Ori suggest a pulsation period of $`0.058\pm 0.002`$.
Using the constraints provided by this preliminary period and by the topology of the instability strip, we have computed linear nonadiabatic models to find the best set of stellar parameters that reproduce the pulsation of V351 Ori. The results are shown in Fig. 4 and the stellar parameters are given in Table 2. The solutions yield a stellar mass of 1.85 M or 2.15 M, respectively, pulsating in the SO (open square in Fig. 4), or the third overtone (TO) mode (open circle in Fig. 4. For lower modes, the luminosity of the model would be lower than the estimated lower limit for V351 Ori, whereas higher modes are probably excluded by the closeness of the observational box to the second overtone blue boundary and we did not consider their occurrence. These solutions would tend to favor a distance of V351 Ori, smaller than that of the young stellar population of the Orion complex.
Recently, Koval’chuk & Pugach (1998) have argued on the basis of several peculiar photospheric properties that V351 Ori is in fact more evolved than previously thought and conclude that this star does not belong to the group of Herbig stars. From Fig. 4, we see that the the post-MS track of a 2 M star intersects the corresponding PMS track at about the position of the best pulsational models. The degeneracy of the tracks does not allow to use the pulsational analysis to discriminate between the two evolutionary phases. However, this uncertainty would disappear if the distance of V351 Ori were the same as that of the Orion population stars. Then, the difference between the pre- and post-MS tracks would be large enough that the period analysis would rule out one of the two solutions. Future studies of V351 Ori should address this important aspect.
As in the case of HD35929, we used the measured Strömgren indices \[c1\] and \[m1\] to provide an independent evaluation of the location in the HR diagram of this star. The results, $`\mathrm{log}T_{\mathrm{eff}}=3.89`$, $`L30L_{}`$, now suggest a stellar mass of $`2.25M_{}`$ for this pulsator, in which case the observed periods would be consistent with an oscillation in the TO pulsation mode.
In conclusion, the present observations yield compelling evidence for the occurrence of $`\delta `$ Scuti-type pulsation in two Herbig Ae stars, HD 35929 and V351 Ori. The comparison with evolutionary and pulsational models provides independent, even if quite preliminary, constraints on the mass and evolutionary state of these stars.

###### Acknowledgements.

We thank the referee, J. Matthews, for his very helpful report. We also thank Dr. R. Silvotti for useful discussions. This research has made use of the Simbad database, operated by CDS, Strasbourg, France.


text:
# Depinning and dynamic phases in driven three-dimensional vortex lattices in anisotropic superconductors

## Abstract

We use three-dimensional molecular dynamics simulations of magnetically interacting pancake vortices to study the dynamic phases of vortex lattices in highly anisotropic materials such as BSCCO. Our model treats the magnetic interactions of the pancakes exactly, with long-range logarithmic interactions both within and between planes. The pancake vortices decouple at low drives and show two-dimensional plastic flow. The vortex lattice both recouples and reorders as the driving current is increased, eventually forming a recoupled crystalline-like state at high drives. We construct a phase diagram as a function of interlayer coupling and show the relationship between the recoupling transition and the single-layer reordering transitions.

In highly anisotropic superconductors such as BSCCO, the vortex lattice is composed of individual pancake vortices that may be either coupled or decoupled between layers depending on such factors as the material stoichiometry or the magnitude and angle of the applied magnetic field. Of particular interest is the possible relationship between a coupling/decoupling transition and the widely studied second peak or fishtail effect .

As a function of interlayer coupling strength $`s`$, there are two limits of vortex behavior in a system containing pointlike disorder. For zero interlayer coupling $`s=0`$, each plane behaves as an independent two-dimensional (2D) system. For infinite interlayer coupling $`s=\mathrm{}`$, the vortices form perfectly straight three-dimensional (3D) lines, and all of the planes move in unison. At finite coupling strength $`s0`$, a transition between these types of behaviors should occur with coupling strength, but it is unclear whether this transition is sharp or if an intermediate state of the lattice exists. Furthermore, it is known that 2D systems with pointlike pinning can exhibit dynamic reordering under the influence of an applied driving current, passing from a liquid-like state at zero drive to a recrystallized state at high current . Thus, in a 3D system, a dynamically driven recoupling transition could be expected, but it is unclear where this transition falls in relation to the 2D reordering transitions already seen.

To study the coupling transitions, we have developed a simulation containing the correct magnetic interactions between pancakes . This interaction is long range both in and between planes, and is treated according to Ref. .

The overdamped equation of motion, $`T=0`$, for vortex $`i`$ is given by $`𝐟_i=_{j=1}^{N_v}𝐔(\rho _{i,j},z_{i,j})+𝐟_i^{vp}+𝐟_d=𝐯_i`$, where $`N_v`$ is the number of vortices, $`\rho `$ and $`z`$ are the distance between pancakes in cylindrical coordinates. The magnetic energy between pancakes is

$`𝐔(\rho _{i,j},0)=2dϵ_0\left((1{\displaystyle \frac{d}{2\lambda }})\mathrm{ln}{\displaystyle \frac{R}{\rho }}+{\displaystyle \frac{d}{2\lambda }}E_1(\rho )\right)`$
$`𝐔(\rho _{i,j},z)=s{\displaystyle \frac{d^2ϵ_0}{\lambda }}\left(\mathrm{exp}(z/\lambda )\mathrm{ln}{\displaystyle \frac{R}{\rho }}E_1(R)\right)`$
where $`R=\sqrt{z^2+\rho ^2}`$, $`E_1(x)=_\rho ^{\mathrm{}}\mathrm{exp}(x/\lambda )/\rho ^{}`$ and $`ϵ_0=\mathrm{\Phi }_0^2/(4\pi \xi )^2`$. The pointlike pins are randomly distributed in each layer and modeled by parabolic traps. We vary the relative strength of the interlayer coupling using the prefactor $`s`$. We have simulated a $`16\lambda \times 16\lambda `$ system containing 89 vortices and 4 layers, with a total of 356 pancake vortices. Further work on systems containing up to 16 layers will be reported elsewhere .

In Fig. 1(a) we present a phase diagram as a function of interlayer coupling strength $`s`$ and driving force $`f_d`$. At zero drive, we find a recoupling transition for a coupling strength of $`s>4.5`$. In samples with $`s5`$, the pancakes remain coupled into lines at all drives and show the same transitions seen in previous work , exhibiting plastic flow of stiff lines above depinning, and reordering into first a smectic state and then a recrystallized state of stiff lines at higher drives.

For samples with weaker interlayer coupling, $`s<5`$, the vortex lattice is broken into decoupled planes at zero drive. Upon application of a driving current, the samples exhibit 2D plastic flow in which each layer moves independently of the others. Once the individual layers reach the driving force at which a transition to a smectic state occurs, the vortices simultaneously form the smectic state and recouple, as can be seen from the measure shown in Fig. 1(c) of the $`z`$-axis correlation $`C_z=1(𝐫_{i,L}𝐫_{i,L+1})\mathrm{\Theta }(a_0/2|(𝐫_{i,L}𝐫_{i,L+1})|)`$, where $`a_0`$ is the vortex lattice constant. The dynamic recoupling transition line follows the smectic transition line down to $`s=2`$ and is associated with a peak in the $`dV/dI`$ curve seen in Fig. 1(b). Both the static and dynamic transition lines between decoupled 2D and recoupled 3D behavior are sharp.

As a function of the number of layers, we observe the same behavior, but the depinning current in the 3D stiff state drops with the number of layers. The depinning current in the 2D decoupled phase is not affected since in this case the individual planes behave as isolated entities. The recoupling transition sharpens as the number of layers is increased .

We acknowledge helpful discussions with L. N. Bulaevskii, A. Kolton, C. Reichhardt, R.T. Scalettar, G. T. Zimányi. This work was supported by CLC and CULAR (LANL/UC) and by the Director, Office of Adv. Scientific Comp. Res., Div. of Math., Information, and Comp. Sciences, U.S. DoE contract DE-AC03-76SF00098.


text:
# Phase synchronization in coupled nonidentical excitable systems and array enhanced coherence resonance

## Abstract

We study the dynamics of a lattice of coupled nonidentical Fitz Hugh-Nagumo system subject to independent external noise. It is shown that these stochastic oscillators can lead to global synchronization behavior without an external signal. With the increase of the noise intensity, the system exhibits coherence resonance behavior. Coupling can enhance greatly the noise-induced coherence in the system.

The study of coupled oscillators is one of the fundamental problems with applications in various fields . Mutual synchronization of the oscillators is of great interest and importance among the collective dynamics of the coupled oscillators. The notation of synchronization has been extended to include a variety of phenomena in the context of interating chaotic oscillators, such as complete synchronization , generalized synchronization , phase synchronization  and lag synchronization . Also, synchronization phenomenon has been studied in stochastic systems. Stochastic resonance can be understood from the view point of frequency locking and phase synchronization of the noise-induced motion to the external signal . Due to this synchronization, stochastic bistable elements subject to the same periodic signal, when appropriately coupled, can display global synchronization to the external periodic signal. This synchronization has the effect of enhancing the stochastic resonance in the array .

Noise-induced coherent motion has been a topic of great interest recently. For a system at a Hopf bifurcation point, an optimal amount of noise can induce the most coherent motion in the system . More generally, noise-induced coherent motion has been demonstrated in a variety of excitable systems  and system with delay . Although described by different terms such as stochastic resonance without external periodic force, autonomous stochastic resonance or coherence resonance, a common feature of this type of systems is the increased coherence of the motion and resonant like behavior of the coherence induced purely by noise without an external signal. An interesting question about this type of systems is how these stochastic elements behave when coupled. Do the elements subjected to independent noise display synchronization? and can the synchronization enhance the coherence in the system motion? A recent study demonstrated noise-induced global oscillation and resonant behavior in a subexcitable media . In Ref. , it is shown that two interacting coherence oscillators can be synchronized. However, whether the coherence of the motion can been enhanced by the interaction between the elements are not clear.

This paper studies these problems by the following $`N`$ coupled nonidentical Fitz Hugh-Nagumo (FHN) neurons, a simple but representative model of excitable system and nerve pulses ,

$`ϵ\dot{x_i}`$ $`=`$ $`x_i{\displaystyle \frac{x_i^3}{3}}y_i+g(x_{i+1}+x_{i1}2x_i),`$ (1)
$`\dot{y_i}`$ $`=`$ $`x_i+a_i+D\xi _i(t),`$ (2)
where $`a_i`$ is a parameter of the $`i`$th element. For a single FHN model, if $`|a|>1`$, the system has only a stable fixed point, while $`|a|<1`$ a limit cycle. The system with fixed point dynamics ($`|a|`$ slightly larger than one) is excitable because it will return to the fixed point only after a large excursion (“near limit cycle”) when perturbed away from the fixed point. To make the study more general, we suppose $`a_i`$ is not the same for the elements, but has a uniform distribution $`a_i(1,1.1)`$. This implementation of the model is physically significant because many physical systems are diffusively coupled, and nonidentity is more natural in physical situations. With nonidentical $`a_i`$, the uncoupled elements will have different response to the same level of noise $`\xi _i(t)`$, e.g., different average firing frequency. The Gaussian noise is uncorrelated in different elements, i.e., $`\xi _i(t)\xi _j(t^{})=\delta _{ij}\delta (tt^{})`$. Periodic boundary condition is employed in our studies. The parameter of the time scale is fixed at $`ϵ=0.01`$.

To characterize synchronization behavior in the lattice of nonidentical excitable systems, we introduce the phase of the elements

$$\varphi _i(t)=2\pi \frac{t\tau _k}{\tau _{k+1}\tau _k}+2\pi k,$$
(3)
where $`\tau _k`$ is the time of the $`k`$th firing of the element defined in simulations by threshold crossing of $`x_i(t)`$ at $`x=1.0`$. The quantity

$$s_i=\mathrm{sin}^2\left(\frac{\varphi _i\varphi _{i+1}}{2}\right)$$
(4)
measures the phase synchronization effect of neighbouring elements . A spatiotemporal average of $`s_i`$, i.e.,

$$S=\underset{T\mathrm{}}{lim}\frac{1}{T}_0^T\left(\frac{1}{N}\underset{i=1}{\overset{N}{}}s_i\right)𝑑t$$
(5)
gives a measure of the degree of phase synchronization in the coupled system. For completely unsynchronized motion $`S0.5`$, while for globally synchronized system $`S0`$.

To measure the temporal coherence of the noise-induced motion, we examine the distribution of the pulse duration $`T_k=\tau _{k+1}\tau _k`$. For a single element subject to noise, the distribution $`P(T)`$ has a peak at a certain value of $`T_k`$ and an exponential tail at large values . A measure of the sharpness of the distribution, for example,

$$R=T_k/\sqrt{\text{Var}(T_k)},$$
(6)
which can be viewed as signal-to-noise ratio, provides an indication of the coherence of the firing event. Biologically, this quantity is of importance because it is related to the timing precision of the information processing in neural systems . For a single element, it has been shown that $`R`$ possesses an optimal value at a certain level of noise . Here we adopt the same measure of coherence, with the distribution $`P(T)`$ constructed by pulse duration of all the $`N`$ elements during a long enough period of time.

Typical behaviors observed in numerical simulations of the system are in the following:

(1) Without external noise, each element comes to a fixed point state and no firing takes place. The fixed points are slightly different due to the nonidentity.

(2) When strong enough noise is added to the originally quiescent system, the system is excited. The firing process of the elements may be in phase synchronization if the coupling is strong enough.

(3) The coherence of the temporal and spatial pattern of the firing process is greatly enhanced in a region of noise and coupling strength.

The main features of the system are illustrated in Fig. 1, where the results of $`R`$ and $`S`$ in the parameter space $`(\mathrm{log}_{10}(g),\mathrm{log}_{10}(D))`$ are shown for a lattice of $`N=100`$ elements. For a fixed coupling $`g`$, $`R`$ increases first with the noise level $`D`$, reaches a maximum, then decreases again, showing the typical resonant behavior without an external signal. In general, for a stronger coupling, a higher level of noise is needed to excite the system. Similarly, for a fixed noise level $`D`$, $`R`$ increases with increasing $`g`$ until it reaches an optimal value; after that, it decrease again.

From Figure 1, one can observe several dynamical regimes in the systems. For very weak coupling ($`g<10^2`$), the firing of the elements are essentially independent, because a noise-induced firing of an element cannot excite its neighbouring lattices. Due to the independent firing, the phase difference has a uniform random distribution on $`(0,2\pi )`$, resulting $`S0.5`$. In this region, with the increase of $`D`$, the temporal behavior of the system display coherence resonance similar to that observed in a single element. The enhancement of the coherence by the noise is not very pronounced. A maximal $`R4.5`$ is found around $`D10^1`$. A typical spatiotemporal pattern of $`x_i`$ and $`s_i`$ for weakly coupled elements subject to relatively weak noise is shown in Fig. 2(a). Both the spatial and temporal behaviors are quite irregular. Each element has its individual average firing frequency due to the nonidentity (smaller frequency for larger $`a_i`$).

With the increase of the coupling strength, the system becomes sensitive to weak noise because the firing events induced by noise now become the source of excitation of the neighbouring elements. This mutual excitation enhances the coherence of the motion in the coupled system, as indicated by increasing $`R`$ and decreasing $`S`$. The lattice displays clusters of synchronization. The clusters breaks and reunites during the evolution, so that each element has slightly different firing frequency. Typical behavior of this partial synchronization regime is shown in Fig. 2(b). However, if the external noise is strong enough so that the noise dominates over the coupling, firing of each element is governed mainly by its individual noise, and synchronization clusters cannot survive, as seen in Fig. 2(c). The coherence of the temporal behavior is relatively low because quite large noise deforms greatly the “near limit” cycle. Note that firing frequency is not locked, but the difference is not as pronounced as in weak noise region (Fig. 2(a)). The next regime where $`R`$ takes large values ($`R18`$) while $`S`$ is very close to zero, is the most interesting, because the system performs quite regular motion globally, as seen in Fig. 2(d). All elements are locked to a relatively large firing frequency, and the distribution of the pulse duration becomes very sharp. After that, with stronger coupling, the system keeps global synchronization, however, the temporal behavior becomes irregular again, as indicated by decreasing $`R`$. This can be understood qualitatively from the global dynamics $`X=x_i_N`$ and $`Y=Y_i_N`$, with $`_N`$ denoting average over the lattice. From Eqs. (1-2), we get approximately

$`ϵ\dot{X}`$ $`=`$ $`X{\displaystyle \frac{X^3}{3}}r^2XY,`$ (7)
$`\dot{Y}`$ $`=`$ $`X+a_0+{\displaystyle \frac{D}{\sqrt{N}}}\xi (t),`$ (8)
where $`r^2=(x_iX)^2_N`$ is the fluctuation level of local dynamics and higher order terms of this fluctuation are ignored; $`a_0=a_i_N`$. The summation of the independent noise is still a Gaussian noise $`\xi (t)`$, but with a weaker strength $`D/\sqrt{N}`$. For strong enough coupling, the system achieves global synchronization, $`x_iX`$, and $`r^20`$, and the coupled system can be viewed as a single element subject to a white noise with strength $`D/\sqrt{N}`$. For a single element, a rather weak noise may not excite the system, or the excited motion is quite irregular . This explains the globally synchronized but irregular motion of the lattice, as shown in Fig. 2(e).

The locations of the above five representative dynamical regimes shown in Fig. 2 have been indicated by the black dots in Fig. 1. These regimes are typical for different size of the lattice. For a larger lattice, the regime of global regular motion ($`R`$ large and $`S0`$) is wider in the parameter space.

Now let us discuss the phenomenon of array enhanced coherence resonance. For a fixed size of lattice and a certain value of coupling $`g`$, there is an optimal level of noise at which the system has a maximal value $`R_{max}`$. Figure 3 shows $`R_{max}`$ as a function of $`g`$ for different size of the lattice. The dashed line represents $`R_{max}`$ for a single element with $`a=1.05`$. As seen from this figure, in the weak coupling region, $`R_{max}`$ of the coupled lattice is lower than that of a single element, because the nonidentity makes the distribution of the pulse duration broader. However, in some intermediate coupling region, even only two coupled elements can enhance the coherence of the temporal behavior, even though they are nonidentical. The enhancement is larger for larger $`N`$. For large enough lattice, $`R_{max}`$ seems to be saturated, but the region of coupling in which $`R_{max}`$ takes large values is broader for larger lattice. As pointed above, for strong enough coupling, the coupled system tends to act as a single element, and $`R_{max}`$ converges to that of the single element at large value of $`g`$. Clearly, for a larger lattice, stronger coupling is needed to make the whole lattice act as a single element due to the local diffusive coupling, resulting in a slower convergence of $`R_{max}`$ to that of the single element.

The properties observed here have some similarity to that in coupled bistable stochastic resonance oscillators subject to the same periodic signal . The difference is that in Refs.  and , global synchronization is observed only at an optimal level of noise, while in the present system, global synchronization occurs for strong enough coupling if the noise level is not too high. The global motion can be regular or irregular. Here we should emphasize that there is no periodic signal in our system, and the coherent motion is purely induced by noise and enhanced by coupling.

To conclude, we demonstrate global synchronization and array enhanced coherence in a lattice of locally coupled, nonidentical FHN model neurons. The results are similar for identical neurons. Interaction between the elements not only renders synchronization of the firing process induced by independent noise, but also improves the coherence of the noise-induced motion greatly. The phenomena demonstrated may be of importance in neurophysiology, where background noise may play an active role to increase the order and timing precision of a large ensemble of interacting neurons in biological information processing.

Zhou Thanks Prof. Gang Hu and Dr. Jinghua Xiao for helpful discussions. This work is supported in part by grants from the Hong Kong Research Grants Council (RGC) and the Hong Kong Baptist University Faculty Research Grant (FRG).


text:
# The superbubble model for LiBeB production and Galactic evolution

## 1. Introduction

Galactic nucleosynthesis and chemical evolution are about how a given element is produced in the universe and how its abundance evolved from the primordial universe on. In the case of the light elements, it is widely agreed that the nucleosynthesis occurs through spallation reactions induced by energetic particles (EPs) interacting with the interstellar medium (ISM). In these reactions, an heavier nucleus (most significantly C, N or O) is ‘broken into pieces’ and transmuted into one of the lighter 6Li, 7Li, 9Be, 10B or 11B nuclei. Except for 7Li, this spallative nucleosynthesis is thought to be the main (if not the only) light element production mechanism. The case of 11B is slightly more complicated, as neutrino-induced spallation in supernovae (the so-called $`\nu `$-process) is sometimes invoked to increase the B/Be and $`{}_{}{}^{11}\mathrm{B}/^{10}\mathrm{B}`$ ratios which one would expect should the light elements be produced by nucleo-spallation alone.

Concerning the Galactic evolution of light element abundances, Fields et al. (2000) have recently re-analyzed the available data as a function of O/H, discussing the uncertainties associated with the methods used to derive the O abundance, the stellar parameters and the incompleteness of the samples. According to their results, Be and B evolution can be described in terms of two distinct production processes: i) a primary process dominating at low metallicity and leading to a linear increase of the Be and B abundances with respect to O – ‘slope 1’ – followed by ii) a secondary process compatible with the standard expectations of the Galactic cosmic ray nucleosynthesis scenario (GCRN) – ‘slope 2’. This behaviour is characterized by a transition metallicity, $`Z_\mathrm{t}(\mathrm{O}/\mathrm{H})_\mathrm{t}`$, below which the Be/O and B/O ratios are constant and above which they are proportional to O/H. Although the value of $`Z_\mathrm{t}`$ is rather uncertain because very few data points have yet been reported at $`Z<Z_\mathrm{t}`$, energetics arguments show that a primary process *is* indeed required below, say, $`10^2Z_{}`$ (Parizot & Drury 1999a,b,2000b, Ramaty et al. 2000), and therefore the very existence of a transition metallicity separating a primary from a secondary evolution scheme seems reasonably well established. In spite of the current large uncertainties on the exact value of $`Z_\mathrm{t}`$, Fields et al. find a range of possible values between $`10^{1.9}`$ and $`10^{1.4}(\mathrm{O}/\mathrm{H})_{}`$ (see also Olive, this conference).

This two-slope picture seems to reconcile the two competing theories for light element nucleosynthesis, namely GCRN which predicts a secondary behaviour for Be and B evolution (e.g. Vangioni-Flam et al. 1990, Fields & Olive 1999), and the superbubble model which predicts a primary behaviour at low metallicity (Parizot & Drury 1999b,2000b, Ramaty et al. 2000). However, we show here that when applied to the whole lifetime of the Galaxy (not only to the early Galaxy), the superbubble model *alone* actually predicts the entire two-slope behaviour inferred from observations, and accounts for all the qualitative and quantitative constraints currently available. Implications for particle acceleration inside a superbubble (SB) are also analyzed.

## 2. Description of the SB model

The superbubble model is based on the observation that most massive stars are born in associations, and evolve quickly enough to explode as SNe in the vicinity of their parent molecular cloud. The dynamical effect of repeated SN explosions in a small region of the Galaxy is to blow large bubbles – superbubbles – of hot, rarefied material, surrounded by shells of swept-up and compressed ISM. The interior of superbubbles consists of the ejecta and stellar winds of evolved massive stars *plus* a given amount of ambient ISM evaporated off the shell and dense clumps passing through the bubble. The exact fraction of the ejecta material inside the SB is not well known, and can be expected to vary with time and from one SB to another. However, this fraction, which we note $`x`$, is all we need to know in order to fully determine the mean composition of the matter inside SBs. Noting $`\alpha _{\mathrm{ej}}(\mathrm{X})`$ and $`\alpha _{\mathrm{ISM}}(\mathrm{X})`$ the abundances of element X among the SN ejecta and in the ISM, respectively, we can indeed write the abundance of X inside the SB as:

$$\alpha _{\mathrm{SB}}(\mathrm{X})=x\alpha _{\mathrm{ej}}(\mathrm{X})+(1x)\alpha _{\mathrm{ISM}}(\mathrm{X}).$$
(1)
The second assumption of the SB model is that the material inside superbubbles is efficiently accelerated by a combination of shocks produced by SN explosions and supersonic stellar winds, secondary shocks reflected by other shocks or clumps of denser material, and a strong magnetic turbulence created by the global activity of all the massive stars. Two different SB models have been proposed so far, assuming different EP compositions and energy spectra. In our model (Parizot & Drury 1999b,2000b), we follow Bykov & Fleishman (1992) and Bykov (1995,1999) and argue that the SB acceleration process produces a rather flat spectrum at low energy, namely in $`E^1`$, as expected from multiple shock acceleration theory (Markowith & Kirk 1999), up to a few hundreds of MeV/n, say. Above this value, the spectrum of the superbubble EPs (SBEPs) is either cut off through a steep power-law or turned into the standard cosmic ray source spectrum (CRS), in $`E^2`$. The exact behaviour of this so-called ‘SB spectrum’ at high energy is important in itself and should be derived from a detailed calculation of the particle acceleration, but we do no consider it here, as it is not relevant to our problem (most of the LiBeB production arises from the most numerous low-energy particles anyway). The other SB model proposed so far (Ramaty & Lingenfelter 1999, Ramaty et al. 2000) assumes that the SBEPs *are* actually the cosmic rays and thus their energy spectrum is the standard CRS spectrum ($`Q(p)p^2`$).

To summarize, the essence of the SB model is that repeated SN explosions occurring in OB associations lead to the acceleration of EPs having either the CRS spectrum or the SB spectrum, and a composition given by Eq. (1), where the only free parameter is the proportion of the ejecta inside the SB: $`x`$. In principle, $`x`$ can be derived from the study of SB evolution dynamics, coupled with a gas evaporation model. But we shall first study LiBeB evolution for itself with no external prejudice about the value of $`x`$, and therefore consider it as a free parameter which we vary from 0 (i.e. SBEPs have the ambient ISM composition) to 1 (i.e. SBEPs are made of pure SN ejecta). Later, we compare the value derived from the LiBeB constraints with the value expected from standard SB dynamical models.

## 3. Be and B Galactic evolution

### 3.1. Qualitative features

Having parameterized our problem as above, we can easily calculate the Be/O production ratio in the Galaxy as a function of $`Z_{\mathrm{ISM}}(\mathrm{O}/\mathrm{H})_{\mathrm{ISM}}`$. We consider SBs blown by 100 SNe exploding continuously over a lifetime of 30 Myr. We then integrate the Be production rates induced by the SBEPs and divide the result by the total O yield (added up assuming a Salpeter IMF and SN yields from Woosley & Weaver 1995). The result is plotted in Fig. 1 for various values of $`x`$ and the two investigated spectra. The main difference between the latter is the Be production efficiency, i.e. the number of Be produced per erg of SBEP. Apart from the SB spectrum being more efficient, both figures show distinctively the sought-for two-slope behaviour, with a transition metallicity $`Z_\mathrm{t}`$ depending on the actual value of $`x`$. This behaviour derives directly from Eq. (1). Replacing X by O there, we see that the abundance of O among the SBEPs is essentially $`x\alpha _{\mathrm{ej}}(\mathrm{O})`$ at low metallicity, and $`(1x)\alpha _{\mathrm{ISM}}(\mathrm{O})`$ above a transition metallicity $`Z_\mathrm{t}\frac{x}{1x}Z_{\mathrm{ej}}`$ (where $`Z_{\mathrm{ej}}10Z_{}`$). Therefore, remembering that O is the main progenitor of Be, we find that the SB model predicts a primary behaviour below $`Z_\mathrm{t}`$ (production efficiency independent of $`Z_{\mathrm{ISM}}`$), and a secondary behaviour above $`Z_\mathrm{t}`$, since the SB model is then essentially identical to the GCRN model (except maybe for the assumed energy spectrum).

Incidentally, it is interesting to note that the SB model can be considered as a correction of the GCRN scenario, taking into account the chemical inhomogeneity of the early Galaxy. Indeed, since particle acceleration occurs precisely in those places where metals are released (i.e. superbubbles), the SBEP composition is considerably richer in O than the average ISM, as long as the SN ejecta dominate the O content of the SBs. Afterwards, it makes little difference, as far as composition is concerned, whether the EPs producing LiBeB are accelerated inside SBs or in the regular ISM.

Finally, we see from Fig. 1 that the predicted slope 1 and slope 2 correlations between Be and O are limit behaviours for very low and very high metallicity respectively. Depending on the value of $`x`$, any intermediate value for the Be-O slope is reached over a given range of stellar metallicity. This is in contrast with what would arise if the two-slope behaviour were to be explained in terms of two different mechanisms (e.g. the SB model at low $`Z`$ and GCRN at high $`Z`$). In that case, indeed, one would have a sharp change of slope at the precise metallicity where the secondary process becomes dominant, with no intermediate values. Of course, expected physical fluctuations of the parameters would weaken this effect, and current observational error bars prevent us from distinguishing conclusively between the two pictures. But we argue that the observed ‘slope 1.45’ behaviour reported by Boesgaard & Ryan (this conference) can be explained (in principle) only if there is a continuous *transition* from slope 1 to slope 2 within a *unique* model (as in the SB model above ), rather than two unrelated models with a slope 2 eventually superseding a slope 1.

### 3.2. Quantitative features

Quantitatively, the Be/O ratio at low metallicity derived from the observations is about $`\mathrm{4\hspace{0.17em}10}^9`$ (Parizot & Drury 2000b). This can be achieve either by the CRS spectrum model, provided $`x50\%`$, or by the SB spectrum model, provided $`x2or3\%`$ (see Fig. 1). So far, both models are equally acceptable since we chose not to accept any prejudice about the value of $`x`$ from outside the restricted field of LiBeB evolution. But when considering the transition metallicity associated with the two possible models, we see that the CRS spectrum implies $`Z_\mathrm{t}>10^1Z_{}`$, well outside the range derived by Fields et al. On the other hand, the value of $`Z_\mathrm{t}`$ predicted by SB spectrum model falls exactly in the required range. As a conclusion, the SB model is fully consistent with the observations provided that i) the SBEP spectrum is flattened at low energy (in $`E^1`$), and ii) the SN ejecta amount to a few percent of all the matter present inside SBs.

Now let us extend the scope of our study. Quite remarkably, the first condition above is exactly what is expected from the SB acceleration model developed by Bykov et al. As for the second condition, it is in perfect agreement with the dynamical model for SB evolution worked out by Mac Low & Mc Cray (1988). In other words, had we looked beforehand for a theoretically preferred value of $`x`$, we would have chosen just the particular value which turns out to account for the various constraints of Be Galactic evolution. Therefore, our results actually bring support not only to the SB model as the natural framework for Be and B evolution studies, but also to the SB acceleration model and standard SB dynamics.

Concerning B, unfortunately, only qualitative constraints can be checked against the SB model (successfully in this instance), since either a significant $`\nu `$-process or a LECR component is required anyway to account for the observed B/Be and $`{}_{}{}^{11}\mathrm{B}/^{10}\mathrm{B}`$ ratios. However, Li does provide additional quantitative constraints. First, in order not to break the Spite plateau, the Li/Be production ratio must be lower than $`100`$. This is shown to be satisfied for any value of $`x`$ greater than about 1% in Fig. (2a). Second, the measurement of the 6Li abundance in two halo stars of metallicity $`Z10^{2.3}Z_{}`$ indicates that the 6Li/9Be ratio in these stars should be in the range 20–80 (see Vangioni-Flam, Cassé, & Audouze 2000 and references therein), in contrast with the solar value of $`6`$. This could not be explained if the proportion of SN ejecta inside SBs were of the order of 50% (CRS spectrum model). However, it is quite remarkable again that the value of a few percent derived from the SB spectrum model is totally consistent with the observed value of the 6Li/9Be ratio, both a low metallicity and at solar metallicity.

## 4. Conclusion

The SB model described above has been shown to be fully consistent with the qualitative and quantitative constraints of LiBeB Galactic evolution: 1) it explains the inferred two-slope behaviour in the framework of one sole model; 2) it provides the correct value of Be/O at low metallicity; 3) it predicts the correct value of the transition metallicity; 4) it does not break the Spite plateau; 5) it is consistent with the 6Li/9Be ratio at any metallicity. Most importantly, these successes rely on the value of only one free parameter, namely the proportion of SN ejecta inside a SB. The value which we find is of the order of a few percent, i.e. exactly in the range derived from standard SB dynamical evolution. Likewise, the SB model is found to be successful only if the SBEPs have the SB spectrum, i.e a flattened shape at low energy (in $`E^1`$). But this is exactly what is predicted by the SB acceleration model of Bykov et al. As a conclusion, the SB model appears to account for all the available constraints about LiBeB evolution by making only the most standard assumptions about the involved models relating to other fields of astrophysics. This may be considered as lending support to these models as well.

#### Acknowledgments.

This work were supported by the TMR programme of the European Union under contract FMRX-CT98-0168.

## References

Bykov A. M. 1995, Space Sci. Rev., 74, 397

Bykov A. M. 1999, in “LiBeB, cosmic rays and gamma-ray line astronomy”, ASP Conference Series, eds. R. Ramaty, E. Vangioni-Flam, M. Casse, K. Olive

Bykov A. M., & Fleishman G. D. 1992, MNRAS, 255, 269

Fields B. D., & Olive K. A., 1999, ApJ, 516, 797

Fields B. D., Olive K. A., Vangioni-Flam E., & Cassé M. 2000, submitted to ApJ (astro-ph/9911320)

Mac Low M. M., & McCray R., 1988, ApJ, 324, 776

Markowith A., & Kirk J. G. 1999, A&A, 347, 391

Parizot E., & Drury L. 1999a, A&A, 346, 686

Parizot E., & Drury L. 1999b, A&A, 349, 673

Parizot E., & Drury L. 2000a, A&A, submitted

Parizot E., & Drury L. 2000b, A&A, submitted

Ramaty R., & Lingenfelter R. E. 1999, in “LiBeB, cosmic rays and gamma-ray line astronomy”, ASP Conference Series, eds. R. Ramaty, E. Vangioni-Flam, M. Casse, K. Olive

Ramaty R., Scully S. T., Lingenfelter R. E., & Kozlovsky B. 2000, ApJ, in press (astro-ph/9909021)

Vangioni-Flam, E., Cassé, M., Audouze, J., & Oberto, Y. 1990, ApJ 364, 586

Vangioni-Flam, E., Cassé, M., & Audouze, J. 2000, submitted to Physics Reports

Woosley S. E., & Weaver T. A. 1995, ApJSS, 101, 181