text:
# Semi-hard scattering of partons at SPS and RHIC: A study in contrast

## Abstract

We analyze the contribution of pQCD based semi-hard parton (re)scattering to the reaction dynamics of relativistic heavy-ion collisions at SPS and RHIC. While such processes are able to account for the measured yield of high momentum direct photons at SPS energies, the conditions necessary for jet-quenching are not fulfilled. The situation changes dramatically at RHIC energies.

preprint: DUKE-TH-02-230

The detection of a quark gluon plasma would confirm one of the important predictions of QCD, namely the de-confinement of quarks and gluons when nuclear matter is heated to a temperature exceeding $`T_c170`$ MeV. Several of the proposed signatures of the formation of a quark-gluon plasma qgp\_rev  in relativistic heavy ion collisions have been observed in collisions of lead nuclei at the CERN SPS. Even though explanations in terms of solely hadronic interactions for many of these observations have been proposed from time to time, the simplest and most consistent interpretation of the observations is in terms of the formation of a de-confined state signs . Microscopic calculations confirm that the energy density reached in these collisions is, indeed, high enough for quark de-confinement to occur weber98 .

The success xnw  of perturbative QCD in describing the transverse momentum spectrum of high-$`p_T`$ pions measured in Pb+Pb collisions at the SPS pions  suggests that the perturbative scattering of nuclear partons contributes to the energy deposition in these reactions. On the other hand, the apparent absence of a visible energy loss by these partons, in stark contrast to the observations at RHIC jet-quench , indicates a dearth of multiple scattering that could lead to local equilibration of the partons before hadrons are formed.

Prompt photons are a good probe of these phenomena, since they escape from the strongly interacting matter without major distortion by final state interactions. We have recently shown that the measured photon yield can be related to the number of semi-hard scatterings among the partons, because a substantial fraction of the photons are emitted in secondary interactions of scattered partons fms ; bms\_phot .

In the present study we report the results of a parton cascade model calculation of photon emission in central collision of Pb nuclei at CERN-SPS energies. We address the following questions:

* How large is the contribution of pQCD based parton scattering to the reaction dynamics of nuclear collisions at the SPS?
* Can this mechanism account for the direct photon yield as measured by the WA98 experiment wa98 ?
* Does the perturbative scattering of partons create conditions necessary for jet quenching, and how do these conditions compare to those predicted by the same model for collisions at RHIC energies?

The parton cascade model (PCM) provides a detailed space-time description of nuclear collisions at high energy, from the onset of hard interactions among the partons of the colliding nuclei up to the moment of hadronization GM . It is based on the relativistic Boltzmann equation for the time evolution of the parton density due to perturbative QCD interactions:

$$p^\mu \frac{}{x^\mu }F_i(x,p)=𝒞_i[F].$$
(1)
The PCM collision term $`𝒞_i`$ is a nonlinear functional of the phase-space distribution function $`F(x,p)`$, containing the matrix elements which account for the following processes:

$$\begin{array}{ccc}gggg\hfill & ggq\overline{q}\hfill & qgqg\hfill \\ qq^{}qq^{}\hfill & qqqq\hfill & q\overline{q}q^{}\overline{q}^{}\hfill \\ q\overline{q}q\overline{q}\hfill & q\overline{q}gg\hfill & \\ qgq\gamma \hfill & q\overline{q}\gamma \gamma \hfill & q\overline{q}g\gamma \hfill \end{array}$$
(2)
with $`q`$ and $`q^{}`$ indicating different quark flavors. The corresponding scattering cross sections are expressed in terms of spin- and color-averaged amplitudes $`||^2`$ Cutler.78 :

$$\left(\frac{\mathrm{d}\widehat{\sigma }}{\mathrm{d}Q^2}\right)_{abcd}=\frac{1}{16\pi \widehat{s}^2}||^2$$
(3)
The total cross section, necessary for the transport calculations is obtained from (3):

$$\widehat{\sigma }_{ab}(\widehat{s})=\underset{c,d}{}\underset{(p_T^{\mathrm{min}})^2}{\overset{\widehat{s}}{}}\left(\frac{\mathrm{d}\widehat{\sigma }}{\mathrm{d}Q^2}\right)_{abcd}dQ^2.$$
(4)
The low momentum-transfer cut-off $`p_T^{\text{min}}`$ regularizes the infrared divergence of the parton-parton cross section. We have used a value of $`Q_0^2=(p_T^{\text{min}})^2`$. We choose the GRV-HO parameterization grv  and sample the distribution functions at the initialization scale $`Q_0^2`$.

Additionally, we include the branchings $`qqg`$, $`qq\gamma `$, $`ggg`$ and $`gq\overline{q}`$ frag . The soft and collinear singularities in the showers are avoided by terminating the branchings when the virtuality of the time-like partons drops below $`\mu _0=1`$ GeV. Some of these aspects were discussed in SG  and explored within framework of the PCM implementation VNI vni . The present work is based on our thoroughly revised, corrected, and extensively tested implementation of the parton cascade model, called VNI/BMS bms\_big1 .

The starting point of our investigation is the comparison of the number and average properties of perturbative parton-parton collisions between nuclear reactions at the SPS and RHIC, displayed in table 1. Following eskola2k  we choose a low-momentum cut-off of 0.7 GeV for our SPS calculation and 1.0 GeV for the RHIC calculation. The number of perturbative parton-parton scatterings at SPS is on average 225 per event – clearly too small to provide a significant effect on the dynamics and characteristics of the bulk of the matter. The number of perturbative scatterings increases by a factor of 14 from SPS to RHIC – at RHIC these interactions are as important for the collision dynamics as the soft, non-perturbative processes. The number of fragmentations increases likewise by more than two orders of magnitude – at SPS it is virtually negligible and has no influence on the parton collision cascade. The final entries of table 1 refer to the average momentum transfer $`Q^2`$ in semi-hard collisions and their average center of mass energy $`\sqrt{\widehat{s}}`$. For the SPS, both values are clearly cut-off dominated. The lack of fragmentations at the SPS and the small values for the average c.m. energy and momentum transfer are a clear indication that the conditions for jet-quenching to occur are not fulfilled at this energy.

A more detailed comparison between the regimes at SPS and RHIC can be seen in Fig. 1. The upper frame shows the $`\sqrt{\widehat{s}}`$ distribution for binary parton-parton collisions at SPS and RHIC. Obviously the collisional regime at SPS is much softer than at RHIC reflecting the strong suppression in the number of semi-hard scatterings observed at the SPS. The two distributions differ strongly not only in magnitude, but also in shape – the distribution at RHIC cannot be obtained by simply rescaling the the distribution for the SPS energy and exhibits strong enhancement towards higher values of $`\sqrt{\widehat{s}}`$. The bottom frame of Fig. 1 shows a comparison of the collision number distributions for SPS and RHIC. Although multiple parton-parton scatterings occur at SPS energy, such re-scattering is strongly suppressed compared to RHIC. E.g., the probability for a parton to scatter five times is more than two orders of magnitude larger at RHIC than at the SPS, and the difference between the two curves grows with the number of multiple scatterings.

In the following, we shall focus on photon production at the SPS, since high-$`k_T`$ photons constitute a direct probe of semi-hard parton scattering bms\_phot . As a first step, we discuss invariant transverse momentum spectra of photons obtained when performing the PCM calculation within the eikonal approximation, as well as with the restriction to primary interactions. These should correspond to a lowest order pQCD estimate (see upper frame of Fig. 2). We consider near central collision of Pb nuclei with impact parameters $`0<b<4.5`$ fm, corresponding to the centrality covered in the WA98 experiment. All the results are for the region of the central rapidity. The close agreement between these three different calculations confirms their insensitivity to the assumptions as well as the accuracy of our numerical implementations. Deviations from these results would then provide evidence for multiple scattering.

The lower frame of Fig. 2 shows calculations when (a) only multiple scattering among the partons, and (b) multiple scatterings as well as fragmentations of the partons are included. Multiple scattering of partons leads to a broadening of the $`k_T`$ distribution for intermediate $`k_T`$ in the range 1.0 GeV $`k_T`$ 2.5 GeV. We find that at SPS the fragmentations do not result in an increasing number of collisions. The increase in photon production, all due to the fragmentation $`qq\gamma `$, is of the of the order of 30%. This result contrasts sharply with the results bms\_phot  obtained for collisions at RHIC where fragmentation processes enhance multiple scatterings by a factor of almost 2.4 (in the absence of LPM suppression).

The effect of higher order corrections via an effective $`K`$-factor or inclusion of an intrinsic $`k_T`$ in the momenta of the partons is investigated in Fig. 3. The intrinsic $`k_T`$ is sampled from a distribution $`\mathrm{exp}(k_T/p_0)`$, with $`p_0=`$ 0.44 GeV vni . While the inclusion of $`K`$-factor increases the yield of photons roughly in proportion (here $`K=2`$), the inclusion of intrinsic $`k_T`$ affects the results in a more complex manner: it reduces the yield of photons with $`p_T<1.2`$ GeV, but increases the yield of photons with larger $`p_T`$ by up to a factor of 2 (we do not consider the results below $`p_T^{\mathrm{min}}`$ which are strongly affected by the cut-off). A similar observation has been made in dumi\_phot , where the importance of nuclear broadening effects for the understanding of the WA98 photon data at the SPS was emphasized.

Figure 4 shows the comparison of our most complete scenario (parton re-scattering, fragmentation and intrinsic $`k_T`$) with data from the WA98 experiment wa98 . The PCM is able to account for the measured direct photon yield for $`k_T2.7`$ GeV. For smaller $`k_T`$ the PCM falls short of the data, which is not surprising since hadronic processes not included in the PCM dominate the photon yield in this range photons .

In summary, our calculations are able to describe the observed photon yield at $`k_T2.7`$ GeV measured by WA98 if higher-order corrections in the form of intrinsic $`k_T`$ are taken into account. Only a small increase in the photon yield is seen when the scattered partons are allowed to fragment (mostly due to $`qq\gamma `$ reactions), since the number of secondary collisions does not change as the fragmentations are rather few. We conclude that the parton medium is dilute and does not form a dense parton plasma, at SPS. This is in stark contrast to the results obtained at RHIC energies bms\_phot  where the multiple scatterings and fragmentations increase the production of photons manifold over the prompt photon yield. Our results suggest that the observed absence of jet quenching at SPS is a consequence of the modest amount of cascading and fragmentation of the scattered partons, while these processes occur in abundance at RHIC energies, resulting in a significant suppression of high-momentum hadrons due to multiple scattering.

###### Acknowledgements.

This work was supported in part by RIKEN, Brookhaven National Laboratory and DOE grants DE-FG02-96ER40945 as well as DE-AC02-98CH10886 and the Natural Sciences and Engineering Research Council of Canada. We thank Charles Gale for useful comments.


format: text

text:
# The effect of bandpass uncertainties on component separation

## 1. Introduction

Small scale anisotropies in the cosmic microwave background can arise from a number of sources. In addition to the ‘primary’ anisotropies, generated at the surface of last scattering, secondary anisotropies and foregrounds can contribute to the observed brightness fluctuations. Separating the components from multi-frequency observations is an important part of the data reduction and interpretation. Here we consider the effect of uncertainties in the frequency response of the instrument within the observational bands and how that impacts our ability to perform component separation.

## 2. Model of the sky

We assume for simplicity that the intensity in any given direction of the sky is the sum of 5 components. The first is the CMB itself with specific intensity

$$I_\nu =\frac{dB_\nu }{d\nu }\frac{x^4e^x}{(e^x1)^2}$$
(1)
where $`B_\nu `$ is a blackbody spectrum and $`x=h\nu /k_BT_{CMB}\nu /56.84`$GHz is the dimensionless frequency. The kinetic Sunyaev-Zel’dovich effect (SZ; Sunyaev & Zel’dovich SZ72 (1972, 1980); for recent reviews see Birkinshaw Bir (1999) and Rephaeli Rep (1995)), arising from the motion of ionized gas with respect to the rest-frame of the CMB, has the same frequency dependence as the CMB signal. The second component is the thermal SZ effect – one of the primary sources of secondary anisotropies in the CMB on small angular scales. Ignoring relativistic corrections, the change in the (thermodynamic) temperature of the CMB resulting from scattering off non-relativistic electrons is

$`{\displaystyle \frac{\mathrm{\Delta }T}{T}}`$ $`=`$ $`y\left(x{\displaystyle \frac{\mathrm{e}^x+1}{\mathrm{e}^x1}}4\right)`$ (2)
$``$ $`2y\text{for }x1,`$ (3)
where the second expression is valid in the Rayleigh-Jeans limit and $`y`$ is the Comptonization parameter which is proportional to the integrated electron pressure along the line of sight.

We also include dust, bremsstrahlung (or free-free) emission and synchrotron radiation following Knox (knox99 (1999), see below). Specifically for the dust we assume the spectral dependence of a modified blackbody with an emissivity index of 2 and a temperature of $`18`$K. For the bremsstrahlung and synchrotron we assume power-laws in frequency with indices $`0.16`$ and $`0.8`$ respectively (Bennett et al. Ben (1992)). Our model does not include a contribution from extragalactic point sources. For the purposes of this analysis we hold the spectral indices of the components fixed (and known). Though we do not expect this to be true of real astrophysical foregrounds, the impact of this on our analysis is negligible.

## 3. Method

To understand the effect of an uncertainty in the response of the instrument to a particular signal let us postulate the following situation. We imagine that the sky is observed at $`N`$ frequencies, with measurements $`\theta _i`$ ($`i=1,\mathrm{},N`$). The signal is the sum of $`M`$ components with amplitudes $`s_\alpha `$ ($`\alpha =1,\mathrm{},M`$) such that

$$\theta _i=\underset{\alpha }{}f_{i\alpha }s_\alpha +n_i$$
(4)
where $`n_i`$ is the noise in channel $`i`$. We write our linear estimator of $`s_\alpha `$ as

$$\widehat{s}_\alpha =\underset{i}{}W_{i\alpha }\theta _i.$$
(5)
Requiring $`\widehat{s}`$ to be an unbiased estimator and minimizing the rms residual, $`(\widehat{s}s)^2`$, gives

$$W_{i\alpha }=\underset{\beta j}{}F_{\alpha \beta }^1f_{j\beta }N_{ji}^1$$
(6)
where $`N_{ij}n_in_j`$ and $`F_{\alpha \beta }=_{ij}f_{i\alpha }N_{ij}^1f_{j\beta }`$ is the Fisher matrix for the $`s_\alpha `$. In the case of diagonal noise, $`N_{ij}=\sigma _i^2\delta _{ij}`$, the Fisher matrix simplifies to

$$F_{\alpha \beta }=\underset{i}{}f_{i\alpha }f_{i\beta }/\sigma _i^2$$
(7)
and the estimator can be written as

$$\widehat{s}_\alpha =\underset{i\beta }{}F_{\alpha \beta }^1\frac{f_{i\beta }}{\sigma _i^2}\theta _i.$$
(8)
The covariance matrix of the statistical errors for this estimator is $`(\widehat{s}_\alpha s_\alpha )(\widehat{s}_\beta s_\beta )=F_{\alpha \beta }^1`$.

We are interested in the effects on our component separation of a band error which causes the actual bandpass to deviate from the design bandpass by $`\delta f`$. Using Eq. (8) with the design bandpass, but replacing $`\theta _i`$ with the right-hand side of Eq. (4) evaluated with the actual bandpass, and subtracting off the true signal we find

$$\delta s_\alpha \widehat{s}_\alpha s_\alpha =\underset{\beta }{}M_{\alpha \beta }s_\beta $$
(9)
where $`M`$ is the ‘component mixing matrix’ given by

$$M_{\alpha \beta }\underset{i\gamma }{}F_{\alpha \gamma }^1\frac{f_{i\gamma }\delta f_{i\beta }}{\sigma _i^2}$$
(10)
and $`\delta f_{i\beta }`$ is the difference between the design and actual bandpasses. The bandpass uncertainty then induces a relative rms error on component $`\alpha `$ from component $`\beta `$ characterized by

$$W_{\alpha \beta }\frac{M_{\alpha \beta }s_\beta ^2^{1/2}}{s_\alpha ^2^{1/2}}.$$
(11)
How is $`\delta f_{i\alpha }`$ related to the bandpass error? If component $`\alpha `$ has frequency dependence $`g_\alpha (\nu )`$ then

$$f_{i\alpha }=𝑑\nu g_\alpha (\nu )\left[h_i(\nu )+\delta h_i(\nu )\right]$$
(12)
where $`h_i(\nu )+\delta h_i(\nu )`$ is the total bandpass, with the latter term the error. It is easy to see that an error in the amplitude of the bandpass, $`\delta h_i(\nu )h_i(\nu )`$, will have no effect on component separation since it will ‘calibrate out’. In general we shall model this calibration process by demanding that

$$f_{i0}=𝑑\nu g_0(\nu )\left[h_i(\nu )+\delta h_i(\nu )\right]=1i$$
(13)
where component $`0`$ is the CMB. We normalize $`g_\alpha (\nu =30\mathrm{GHz})=1`$ for all components $`\alpha `$. Note that this means $`s_\alpha `$ is the amplitude of component $`\alpha `$ at 30 GHz.

Note that the above treatment applies in both pixel space and spherical harmonic space. In pixel space the $`\sigma _i`$ should all be calculated for the same pixel size. In spherical harmonic space, the $`\sigma _i`$ (interpreted as errors on the beam–deconvolved maps) are $`\mathrm{}`$–dependent:

$$\sigma _i(\mathrm{})=\widehat{\sigma }_i\vartheta _i\mathrm{exp}\left[\frac{1}{2}\left(\frac{\mathrm{}\vartheta _i}{2.355}\right)^2\right]$$
(14)
where $`\vartheta `$ and $`\widehat{\sigma }`$ are the beam and noise defined in Table 1 with $`\vartheta `$ converted to radians.

## 4. Results

As an example let us consider simple shifts in the bandpass, which we shall model as a top-hat of width $`\mathrm{\Delta }\nu `$ about the central frequency as given in Table 1. To gain intuition let us restrict ourselves to a $`2\times 2`$ subspace consisting only of CMB and tSZ signals with a $`1`$GHz shift in the $`217`$GHz channel to higher frequency. In this case (at angular scales larger than the largest beam)

$$M_{\alpha \beta }=\left(\begin{array}{cc}0.0& 0.00391\\ 0.0& 0.00522\end{array}\right)$$
(15)
which indicates that the CMB channel does not ‘contaminate’ the any of the signals but there is leakage from the tSZ channel. That $`M_{\alpha 0}=0`$ is a direct result of our calibration procedure which enforces $`\delta f_{i0}=0`$. Note that for other calibration sources typically used by ground–based experiments, e.g. planets, the CMB can, via bandpass errors, contaminate other components.

A more realistic scenario is that the bandpass frequency be quite well determined, but the amplitude of the response as a function of frequency be somewhat uncertain. This holds for both HEMTs (M. Seiffert, private communication) and Bolometers (P. Ade, private communication). As a simple model of this effect we introduce a linear ‘ramp’ into our otherwise top-hat bandpasses with a change in amplitude of 10% across the band. The results of including a $`+10\%`$ ramp in all of the channels with our 5 component model is given in Table 2.

The elements of the mixing matrix in Table 2 tell us, for example, that our CMB map, $`\widehat{s}_0`$ will have a contribution from dust of $`0.25s_2`$, where $`s_2`$ is the dust amplitude at 30 GHz. Fortunately the dust amplitude is very low at 30 GHz!

The importance of this contamination is easier to read from the normalized mixing matrix, $`W_{\alpha \beta }`$. To calculate $`W_{\alpha \beta }`$ we need to know the rms sky fluctuations of the various components. We use the model described in Knox (knox99 (1999)) and references therein, most notably Bouchet & Gispert (bouchet99 (1999)). The exception is our SZ power spectrum which we take from White, Hernquist & Springel (WHS (2002)) and extend to $`\mathrm{}<400`$ by assuming $`C_{\mathrm{}}=C_{400}`$ appropriate for the Poisson dominated regime (see Fig. 1).

In Table 3 we assume $`s_\alpha `$ are the amplitudes of an $`\mathrm{}=500`$ spherical harmonic, so that $`s_\alpha ^2=C_{\mathrm{}}^\alpha `$. We find that all the elements are quite small; the largest contribution to the CMB map comes from SZ and the rms of this contaminant is .05% of the CMB signal rms. The galactic contaminants of the CMB map, $`W_{0\beta }`$ for $`\beta =2,3,4`$ greatly increase with increasing angular scale. Even so, they are very small at all angular scales; at $`\mathrm{}=2`$ $`W_{0\beta }`$ = -0.006, 0.0002 and $`2\times 10^5`$ for dust, synchrotron and bremstrahlung respectively.

In Fig. 1 we show how the SZ contamination of the CMB map affects the CMB power spectrum and the CMB-SZ cross–correlation power spectrum. Denoting the CMB contaminant from SZ as $`a_{lm}^{\mathrm{CMB}(\mathrm{SZ})}`$ we have

$$C_{\mathrm{}}^{\mathrm{CMB}(\mathrm{SZ})\mathrm{SZ}}a_\mathrm{}m^{\mathrm{CMB}(\mathrm{SZ})}(a_\mathrm{}m^{\mathrm{SZ}})^{}=M_{01}(\mathrm{})C_{\mathrm{}}^{\mathrm{SZ}\mathrm{SZ}}.$$
(16)
The $`\mathrm{}`$-dependence of the mixing matrix arises from the $`\mathrm{}`$-dependence of the noise in the beam–deconvolved maps but is quite mild: $`M(0,1)`$ monotonically decreases from $`4\times 10^3`$ to its high $`\mathrm{}`$ asymptote of $`5.3\times 10^3`$. In general the cross–correlations are the most affected since they are first order in the component mixing matrix. Fortunately, for the 10% ramp error, the bandpass error–induced cross–correlation is well below the level of the Planck noise.

Our results are specific to the conservative component separation procedure we have assumed. Other methods can reduce the statistical noise by including assumptions about the statistical properties of the various components (Tegmark & Efstathiou tegmark96 (1996); Bouchet, Gispert & Puget bouchet95 (1995); Hobson et al. hobson98 (1998)). The results will also differ in detail if one adopts a more realistic (and more complicated) model of the foregrounds. However, we do not expect the component mixing matrix to be qualitatively different for these different procedures or for more realistic foreground modeling. The fact that the two–component model and the five–component model give similar results for both $`M(0,1)`$ and $`M(1,1)`$ we take as evidence of this robustness.

The bandpass uncertainties can lead to systematic errors in the velocities of clusters inferred from the kinetic SZ effect. The mixing between thermal SZ and primary CMB on a cluster of optical depth $`\tau `$ and temperature $`T_e`$ induces an error

$$\delta T8\times 10^3\tau \left(\frac{kT_e}{m_ec^2}\right)$$
(17)
in the CMB signal. If the contamination were all erroneously attributed to kinetic SZ from the moving cluster, it would bias the inferred velocity by $`v=8\times 10^3(kT_e/m_ec^2)c50\mathrm{k}\mathrm{m}\mathrm{s}^1`$ for a rich cluster. This bias is negligible compared to other sources of uncertainty for individual cluster velocities (Haehnelt & Tegmark (1996); Nagai et al. (2002); Holder (2002)) but is comparable to errors that might be achievable by Planck on bulk flows in $`10^6h^3`$ Mpc3 volumes (Aghanim et al. (2001)). This systematic contaminant of the bulk flows will appear as a $`T_e`$–weighted monopole. Such a pattern is not expected cosmologically, would be evidence of bandpass errors and could be removed from the data with negligible residuals. Note that while we have ignored relativistic corrections to the SZ distortion for the purposes of estimating the magnitude of the effect, they will be important in the actual analysis of the data (Diego et al. (2002)).

## 5. Conclusions

We can conclude that bandpass errors at the $`𝒪(10\%)`$ level are acceptably small for Planck. More sensitive experiments though may have more stringent requirements on the quality of the bandpass measurements because an overall scaling of the sensitivity of each channel leaves the component mixing matrix unchanged.

We thank P. Ade, G. Holder, C. Lawrence and M. Seiffert for useful conversations. M.W. was supported by a NASA Astrophysical Theory Grant, the NSF and a Sloan Fellowship. L.K. was supported by NASA NAG5-11098. We thank the KITP (supported by NSF PHY99-07949) for their hospitality.


text:
# Dynamical Brittle Fractures of Nanocrystalline Silicon using Large-Scale Electronic Structure Calculations

## Abstract

A hybrid scheme between large-scale electronic structure calculations is developed and applied to nanocrystalline silicon with more than 105 atoms. Dynamical fracture processes are simulated under external loads in the  direction. We shows that the fracture propagates anisotropically on the (001) plane and reconstructed surfaces appear with asymmetric dimers. Step structures are formed in larger systems, which is understood as the beginning of a crossover between nanoscale and macroscale samples.

Silicon is one of ideally brittle materials and its fracture behavior is studied intensively, because we can obtain essentially dislocation-free single crystals. A pioneering theory of brittle fracture was given, in 1920’s, by Griffith GRIFFITH  within a continuum theory, which is the foundation of the present understanding of brittle fractures BRITTLE-TEXTBOOK . The fracture in single crystals should be investigated also in atomistic pictures, on the point of how and why the fracture path is formed and propagates in the crystalline geometry. This point includes surface reconstruction processes. Since fracture is a thermal non-equilibrium process, the atomic structure on a cleavage surface can be different from that on equilibrium clean surfaces. For example, the easiest cleavage plane in macroscale samples of silicon is the (111) plane, in which the surface structure is not the ground state ($`7\times 7`$) structure but a metastable $`2\times 1`$ structure PANDEY1 ; SI111-21-PARRINELLO-MD .

This letter is devoted to the atomistic fracture behaviors in nanocrystalline silicon, especially, its possible difference from macroscale samples. Such a difference can be expected, as explained below; Now a typical atomistic length in silicon is defined as $`d_0^3\sqrt{v_0}3`$ Å, where $`v_0`$ gives the volume per atom. The essence of the Griffith theory GRIFFITH  is the energy competition between the energy gain of the strain relaxation and the loss of the surface formation energy. The former energy is a volume term proportional to (length)3, while the latter energy is a surface term proportional to (length)2. As analogous to the theory of nucleation LANDAU , the dimensional analysis gives the critical crack length for the spontaneous fracture propagation. The critical crack length $`c_\mathrm{G}`$ is given as GRIFFITH ; BRITTLE-TEXTBOOK

$`c_G{\displaystyle \frac{\gamma E}{\sigma ^2}}`$ (1)
with the stress $`\sigma `$, the Young modulus $`E(10^2\mathrm{GPa})`$ and the surface energy per unit area $`\gamma `$. The value of $`\gamma `$ was estimated to be in the order of $`1\mathrm{J}/\mathrm{m}^2`$ SPENCE93 ; GUMBSCH2 , which can be reduced to the bond breaking energy $`(\gamma d_0^21\mathrm{e}\mathrm{V})`$ in the atomistic picture. In a recent experiment with macroscale samples CRAMER , the stress is $`\sigma 10^1`$ MPa and Eq. (1) gives a macroscale length ($`c_G1`$ mm). Since the length $`c_\mathrm{G}`$ is not dependent on the sample size $`L`$, the fracture behavior can be expected to be different from the above picture, in case that the sample size $`L`$ is smaller than the critical length $`c_\mathrm{G}`$ $`(L<c_\mathrm{G})`$. In this letter, we will see such a situation in nanocrystalline silicon, in which the numbers of atomic layers for these lengths $`(c_\mathrm{G}/d_0,L/d_0)`$ are not macroscale numbers.

For atomistic fracture simulations of silicon crystals, a recent work of classical modelings HOLLAND-MARDER  was done with $`10^5`$ atoms. A more recent workGUMBSCH2 , however, pointed out the limited applicability of classical modelings and the importance of electronic structure calculations. On the other hand, there are several ab initio calculations with $`10^2`$ atoms SPENCE93 ; GUMBSCH2 . Due to the system size of simulations, these investigations are limited in situations, such as the preparation of the initial cleavage plane in which the reconstructed surface structure is assumed. Therefore, large-scale electronic structure calculations are essential.

So far, we have developed several order-N methods for large-scale electronic structure calculations HOSHI2000A . The order-N method is the general name of methods in which the computational cost is proportional to the system size ($`N`$). We have developed the variational and perturbative order-N methods based on generalized Wannier states HOSHI2000A ; HOSHI2001A . The Wannier states $`\{\varphi _i\}`$ are localized and the index $`i`$ denotes its localization center. The equation for wave functions is common between the two methods and is given by the one-body density matrix $`(\rho _i^{\mathrm{occ}.}|\varphi _i\varphi _i|)`$ HOSHI2000A . In the computational algorithm, the perturbative method is simpler than the variational method. Figure 1 demonstrates our large-scale calculation among $`10^210^6`$ atoms.

With the above two methods, we now construct a novel hybrid scheme, in which the variational method is used only for the wave functions whose centers locate near fracture regions. The regions contain, typically, $`4\times 10^4`$ electrons. Some of such wave functions change their character dynamically from the bulk (sp3 bonding) states to surface ones, as discussed later. The other wave functions, in bulk regions, keep the character of the bulk bonding states and can be obtained by the perturbative method. The wave functions $`\{\varphi _i\}`$ calculated by the two methods are commonly used for constructing the one-body density matrix $`\rho `$. Any physical quantity is expressed by the density matrix HOSHI2000A  and is well defined in the present hybrid scheme.

The present work is based on a transferable tight-binding Hamiltonian with s and p orbitals KWON . It is used for several crystalline phases and non-crystalline phases, such as liquid KWON  and surfaces SURFACE-TBMD . Since the fracture is the formation of surfaces in a bulk region, the theory should reproduce the atomic structures both in bulk and surface phases, which is satisfied in the present Hamiltonian. The essence of the quantum mechanical freedoms is the fact that the sp3-hybridized bonds are formed in the bulk region, while are not on surfaces. To analyze the hybridization freedom, a parameter $`f_\mathrm{s}^{(j)}`$ is defined, for a wave function $`\varphi _j`$, as

$`f_\mathrm{s}^{(j)}{\displaystyle \underset{I}{}}|\varphi _j|Is|^2,`$ (2)
where $`|Is`$ is the s orbital at the $`I`$-th atom. For example, $`f_\mathrm{s}=1/4`$ in an ideal sp3 hybridized state.

In this letter, we focus on the Si(001) surface, a standard template of the modern silicon technology. A characteristic feature in the Si(001) surface is the formation of asymmetric dimers CHADI79-SI001 ; RAMSTAD . The asymmetric dimer is connected by a ‘$`\sigma `$’ bonding state. Another state is localized on the ‘up’ atom, the dimerized atom near the vacuum region. This localized state is called ‘$`\pi `$’ state, because the direction of its p components is nearly perpendicular to the dimer bond. Here an energy quantity is defined as

$`\mathrm{\Delta }\epsilon _i^{(\mathrm{cov})}\varphi _i|H|\varphi _i\left[f_\mathrm{s}^{(i)}\epsilon _\mathrm{s}+(1f_\mathrm{s}^{(i)})\epsilon _\mathrm{p}\right].`$ (3)
A negative value of $`\mathrm{\Delta }\epsilon _i^{(\mathrm{cov})}`$ corresponds to the energy gain of a covalent bonding. The ‘$`\sigma `$’ state has the gain of $`\mathrm{\Delta }\epsilon _i^{(\mathrm{cov})}2\mathrm{e}\mathrm{V}`$, which mainly contributes to the dimerization energy (about $`2\mathrm{e}\mathrm{V}`$)RAMSTAD . The ‘$`\pi `$’ state has much smaller $`\mathrm{\Delta }\epsilon _i^{(\mathrm{cov})}`$, which is comparable to the energy difference between the asymmetric and symmetric dimers (the order of $`0.1\mathrm{eV}`$)RAMSTAD .

The simulation details are as follows; The hybrid order-N scheme is used for systems with 104 atoms or more. In systems smaller than the above size, the variational method is used in the whole region. The samples are isolated tetragonal clusters, whose geometries are labeled with the number of atomic layers, such as $`n_{100}\times n_{010}\times n_{001}`$ or $`n_{110}\times n_{1\overline{1}0}\times n_{001}`$. As the boundary condition, the Wannier states at the sample surfaces are terminated by fixed sp3 bonding states and are not reconstructed. The time step of the molecular dynamics is 3 fs. The total kinetic energy is controlled to be that with 300 K by the Nose thermostat method NOSE . The numerical accuracy is checked among bulk, surface and fracture properties, such as the elastic constants, the dimer formations on the clean (001) surface, and the critical stress for fracture. The last quantity is checked, with smaller samples, in comparison with the standard diagonalization. The calculated fracture propagation velocity is always in the same order of, but less than, the Rayleigh surface wave velocity ($`4.5`$km/s), as expected from the continuum theory BRITTLE-TEXTBOOK .

For fracture propagations, external loads in the  direction are imposed. During the simulations, the external loads can be dynamically controlled by the atoms on the sample surfaces in the $`z`$ direction. These atoms are fixed or under artificial constant-velocity motions in the $`z`$ axis. The velocity, typically $`10^2`$ km/s, is much smaller than that of observed fracture propagation velocities (km/s). As a seed of fractures, a short range repulsive potential is imposed on one particular pair of atoms, as a defect bond. For smaller samples, the simulations begin without initial deformations. The fracture always occurs with the external loads in the order of $`\sigma 1`$ GPa, which corresponds to the strain energy of $`\sigma d_0^30.1\mathrm{eV}`$ per atom. For larger samples, the simulations begin with initial static deformations in the above magnitude of external loads. The length $`c_\mathrm{G}`$ in Eq. (1) is calculated as $`c_G100\mathrm{n}\mathrm{m}`$, which is longer than the present sample sizes ($`L20\mathrm{n}\mathrm{m}`$).

In results, a two-stage reconstruction process is commonly observed as the elementary process during successive bond breakings; In Fig. 2(a), we monitor the one-electron energy $`\epsilon _i\varphi _i|H|\varphi _i`$ and the hybridization freedom $`f_\mathrm{s}^{(i)}`$ of a Wannier state $`|\varphi _i`$. Before the bond breaking ($`t<0\mathrm{ps}`$), the wave function $`|\varphi _i`$ is a bonding state in the bulk region, deformed due to the external load. At $`t0\mathrm{ps}`$, a bond breaking occurs and the wave function $`|\varphi _i`$ loses the bonding character with rapid increase of the bond length. Then ($`0\mathrm{ps}<t<0.2\mathrm{ps}`$), a twofold coordinated surface atom appears, since another bond is broken almost simultaneously. The wave function $`|\varphi _i`$ forms a loan pair state that is stabilized by the increase of $`f_\mathrm{s}^{(i)}`$ $`(0.60.8)`$. The corresponding energy gain can be estimated to be $`0.2\times (\epsilon _\mathrm{p}\epsilon _\mathrm{s})1.3`$eV, which explains the energy gain in the figure ($`\epsilon _i=2.7\mathrm{eV}3.8\mathrm{eV}`$). In other words, the bond breaking process is caused by the local electronic instability, that is, the energy competition between the loss of the bonding (transfer) energy and the gain due to the increase of the weight on the s orbitals ($`f_\mathrm{s}`$). Finally, after the thermal motions with a finite time ($`t0.4\mathrm{ps}`$), a pair of twofold coordinated atoms forms an asymmetric dimer with a $`\sigma `$ bonding state $`|\varphi _i`$. The corresponding covalent-bonding energy, defined in Eq. (3), is $`\mathrm{\Delta }\epsilon _i^{(\mathrm{cov})}1.9`$ eV. This energy explains the gain in the figure ($`\epsilon _i=3.8\mathrm{eV}4.8\mathrm{eV}`$) and the energy loss (about 1.3eV) due to the decrease of $`f_\mathrm{s}^{(i)}`$ $`(0.80.6)`$. This asymmetric dimer is preserved until the end of the simulation, during a couple of pico seconds. Figure 2(b) is an example of observed asymmetric dimers.

Figure 3 shows the fracture process of a cubic sample with 4501 atoms. Each Wannier state is classified from its weight distribution into a bonding or atomic orbital, which is shown as a rod or a ball in the figures, respectively. A black rod or ball corresponds to one in the layer that contains the defect bond. One almost flat (001) surface is being created with many asymmetric dimers. The surface contains, however, many twofold coordinated atoms that have two back bonds (white rods) and a lone pair orbital (black ball). This is because a lone pair states are metastable, as discussed above. In Fig. 3, an anisotropic bond-breaking propagation is seen in the  and \[1$`\overline{1}`$0\] directions, especially in the early snapshots. In the $`[110]`$ direction, the successive bond breakings propagate along the nearest neighbor bond sites, which forms a zigzag path, as the black rods in Fig. 4(a). A bond breaking process drastically weakens the nearest neighbor bonds, due to the local electronic instability, observed in Fig. 2(a). Therefore, the successive bond breakings propagate easily in the $`[110]`$ direction. In the $`[1\overline{1}0]`$ direction, on the other hand, the bond-breaking paths are not connected, as the red rods in Fig. 4(a). In this direction, the bond breakings are propagated through the local strain relaxation, not by the local electronic instability. As results, the bond breaking propagation along the nearest neighbor bond sites (in the $`[110]`$ direction of the present surface) is faster than that in the perpendicular direction (in the $`[1\overline{1}0]`$ direction), due to the difference of the successive bond breaking mechanisms. Note that a flat (001) surface is also obtained by a similar simulation without the initial defect bond, in which the fracture begins at the sample edges.

Figures 4 (b) and (c) show larger samples with step formations NOTE-STEP1 . In the two cases, all the conditions are the same, except the sample sizes. To see the step structures clearly, the broken bond sites are shown as rods in the ideal crystalline geometry. The defect bond is located in the center of the drawn area. The anisotropic fracture propagation in one (001) plane increases the anisotropic strain energy NOTE-ELAS-ISO . The anisotropy originates from the inequivalence between the  and \[1$`\overline{1}`$0\] directions within one (001) layer. Since the above inequivalence does not appear within two successive layers, a step formation between them will release the anisotropic strain energy. In Fig. 4 (b), a step is formed between the layer of black rods and that of red rods. In the $`[110]`$ direction, the bond-breaking propagation reaches the sample surfaces without step formations. In the \[1$`\overline{1}`$0\] directions, the bond-breakings propagate slower and a step is formed in the central area at an early period of the crack propagation. After that, the fracture propagates among the two atomic layers. Since the two layers are symmetrically equivalent, the resultant step formation path, is almost a line in the  or  directions, as the boundary of the fractured areas between the two layers.

In Fig.4(c), the largest sample in the present letter, the above line structure does not reach the sample surfaces but is canceled with additional step formations in complicated paths. The sample size dependence of the step structures is understood by the beginning of the crossover between nanoscale and macroscale samples; If the sample contains so many atoms, the geometry of the resultant crack will be almost circular, as in Fig. 4 (c), so as to minimize the anisotropic strain energy NOTE-ELAS-ISO . If not, the strain energy is accumulated only within the confined bulk region due to the finite sample size. The resultant fracture behavior is directly related to the anisotropic atomic structure of the cleaved surface, as in Fig. 4 (b).

Since the above mechanism of step formations is two-dimensional, the present samples may be nanoscale ‘thin’ samples. In larger or thicker samples, an expected fracture behavior is the bending of the fracture plane into the (111) plane, the easiest cleavage plane in macroscale samples, which is the crossover in the present context. In an enough large sample, the fracture mode with the easiest cleavage plane will grow with no regard for sample shapes and details of conditions. Note that the dynamical simulation with 105 atoms is the practical limitation within a single CPU workstation and the program code with parallel computations is now being developed for simulations with larger samples.

This letter shows a possible difference in fracture behaviors among nanoscale and macroscale silicon crystals. Its origin is the size dependence of the energy competition between bulk and surface regions. The electronic structures between the two regions are essentially different and can be described by the present method with the well-defined total energy. This energy competition is also inherent in other phenomena, such as crystal growth and self organizations, which may be candidates for applications of the present method.

This work is supported by a Grant-in-Aid for COE Research ‘Spin-Charge-Photon’ and a Grant-in-Aid from the Japan Ministry of Education, Science, Sports and Culture. This work is also supported by ‘Research and Development for Applying Advanced Computational Science and Technology’ of Japan Science and Technology Corporation.


text:
# Stellar structure and mass loss during the early post asymptotic giant branch

## 1 INTRODUCTION

Models for the acceleration of winds from asymptotic giant branch (AGB) stars which are based on radiation pressure on dust agree well with observations (e.g., Elitzur, Ivezic & Vinkovic 2002). The exact mass loss rate depends on the process which initiates the outflow (Elitzur et al. 2002). One of the most popular mechanism to initiate the outflow, which is based both on observations and theoretical studies, is stellar pulsation (e.g., Wood 1979; Bowen 1988 Höfner et al. 1996; Höfner 1999; Winters et al. 2000; Andersen, Höfner & Gautschy-Loidl 2002). In such models, it is found that the mass loss rate substantially increases as the star evolves along the AGB (e.g., Bedijn 1988; Bowen & Willson 1991; Höfner & Dorfi 1997; Wachter et al. 2002, hereafter WA02), with strong dependence on the stellar mass, effective temperature, and luminosity (or radius). This behavior along the AGB is in agreement with observations which require the mass loss rate to substantially increase near the termination of the AGB, in what was termed superwind (Renzini 1981). The strong dependence on the temperature (e.g., WA02) implies that as the star evolves along the post-AGB track and becomes hotter, at constant luminosity, the mass loss rate steeply decreases with time. However, observations of planetary nebulae (PNs) and stellar evolution calculations along the post-AGB require the high mass loss rate to continue during the post-AGB phase (e.g., Tylenda & Stasinska 1994), until the star heats-up to $`T5000\mathrm{K}`$ (Blöcker 1995). Indeed, the post-AGB evolution time obtained by WA02 is much too long (R. Szczerba 2002, private communication). We elaborate on this in Section 2.

It seems that additional physical parameter(s) are needed to determine the mass loss rate (temperature still should be a factor as it influence dust formation). This is also supported by the observations that some elliptical PNs have spherical halo$``$which was presumably formed from the AGB ‘regular’ wind$``$while the inner dense region$``$which was probably formed from a superwind$``$is axisymmetrical. The change in mass loss geometry suggests that other factors are involved in determining the mass loss rate and geometry. Soker & Harpaz (1999, hereafter SH99) propose that the density and entropy profiles are the additional factors. Because of the high mass loss rate the envelope mass decreases very fast, such that even when the star starts its contraction on the post-AGB early phase, the average density of the stellar envelope decreases. This makes the density profile below the photosphere much shallower, and the entropy profile much steeper. SH99 suggest that one possibility by which these profiles enhance mass loss rate is that they enhance magnetic activity. Other possibilities are a switch in the pulsation mode, possibly to include non-radial pulsation modes, or a more violent convection, because of the steeper entropy profile. A more violent convective region may influence the pulsation and deformation of the surface of the AGB star, and cause shocks to propagate through its photosphere (Jacobs et al. 2000), possibly increasing the mass loss rate. These mechanisms require more investigation.

SH99 concentrate on the onset of the superwind until the stellar radius shrinks to $`R260R_{}`$ and the temperature reaches $`T3300\mathrm{K}`$. In the present paper we show that the same arguments used by SH99 for the onset of the final intensive wind$``$the superwind$``$imply that the superwind ends when the effective temperature is around $`T5000\mathrm{K}`$. We emphasize that we do not propose a new mass loss mechanism. We accept that pulsations coupled with radiation pressure on dust is the mechanism for mass loss, and that stellar mass, effective temperature, and radius play a significant role. However, we suggest (Sect. 3) that the envelope structure, in particular the entropy and density profiles, are important physical factors in determining the mass loss rate on the early post-AGB phase.

## 2 PROBLEMS WITH STRONG TEMPERATURE DEPENDANCE

In the physical models cited above for mass loss from upper AGB stars (e.g., Bowen & Willson 1991; WA02), the mass loss rate depends mainly on the stellar effective temperature, $`T`$, the stellar luminosity, $`L`$ (or stellar radius $`R`$) , and the total stellar mass $`M_{}`$. These models cannot reproduce the correct mass loss rate after the star evolves to the left on the HR diagram, i.e., when it becomes hotter at a constant luminosity. The reason is as follows. In order to get a very large increase in the mass loss rate as stars approach the tip of the AGB, as inferred from observations, the mass loss rate should strongly depend on the different variables, and in particular on the temperature. WA02 find the mass loss rate, in $`M_{}\mathrm{yr}^1`$, to be

$$\mathrm{log}\dot{M}_T=8.861.95\mathrm{log}M_{}/M_{}\beta \mathrm{log}T/K+2.47\mathrm{log}L/L_{},$$
(1)
with $`\beta =6.81`$. The subscript $`T`$ stands for winds where the stellar effective temperature is the main parameter determining the mass loss rate. An AGB star starts its blueward motion (to the left on the HR diagram) when its envelope mass is $`M_{\mathrm{env}}0.10.4M_{}`$ (Schönberner 1983; Soker 1992); stars with a more massive core turn blueward while retaining a more massive envelope. With the strong dependence on the temperature, the mass loss rate decreases very fast. The result is that the evolution time between the maximum mass loss rate epoch and the PN stage is much too long.

We now show that the post-AGB evolution time according to the mass loss formula of WA02 is too long. In the model calculated by Soker (1992) the star reaches the tip of the AGB when its envelope mass is $`M_{\mathrm{env0}}0.16M_{}`$, then stays at that radius until the envelope mass decreases to $`M_{\mathrm{env1}}0.08M_{}`$. After that point the envelope shrinks according to

$$R400\left(\frac{M_{\mathrm{con}}}{0.1M_{}}\right)^{0.2}R_{},\mathrm{for}M_{\mathrm{con}}10^5M_{},$$
(2)
where $`M_{\mathrm{con}}`$ is the mass in the convective region. To a good accuracy the equality $`M_{\mathrm{con}}M_{\mathrm{env}}`$ holds until the envelope mass reduces to $`M_{\mathrm{env}}0.002M_{}`$, the temperature rises to $`4000\mathrm{K}`$, and the stellar radius shrinks to $`R180R_{}`$. The stellar luminosity in that model is $`L=7500L_{}`$, hence the temperature is $`T=2700(R/400R_{})^{1/2}K`$. During the later phase, the stellar mass does not change much, and stays in the range of $`0.60.68M_{}`$. Substituting the temperature, luminosity, and mass at the lower end in the expression for the mass loss rate of WA02 (eq. 1), and using equation 2, gives for the mass loss rate

$$\dot{M}_T3\times 10^5\left(\frac{M_{\mathrm{env}}}{0.1M_{}}\right)^{0.1\beta }M_{}\mathrm{yr}^1,$$
(3)
where $`\beta `$ is defined in equation 1. The total time required to lose an envelope mass of $`M_{\mathrm{env1}}`$, during which approximation (2) holds, is given then by (assuming $`\beta <10`$)

$$\tau _{T1}=\frac{dM_{\mathrm{env}}}{\dot{M}_T}3300\frac{1}{10.1\beta }\left(\frac{M_{\mathrm{env1}}}{0.1M_{}}\right)^{10.1\beta }\mathrm{yr}.$$
(4)
Because $`\beta `$ is close to 10, the factor $`1/(10.1\beta )`$ is large, and the evolution time is too long. For the parameters cited above with $`\beta =6.81`$, it takes $`\tau _{T1}10^4\mathrm{yrs}`$ to lose the last $`0.08M_{}`$. This evolution is too long, as the dynamical age of the dense shell of many PNs is much below that value (Sect. 1). Also, the maximum mass loss rate according to equation (1), when applied to the model from Soker (1992), is reached already when $`M_{\mathrm{env0}}=0.16M_{}`$. The evolution time from that point to the point when $`M_{\mathrm{env1}}=0.08M_{}`$ is $`\tau _{T0}=(0.160.08)/3\times 10^5=2700\mathrm{yr}`$. Overall, the PN stage will start $`12000\mathrm{yrs}`$ after the superwind has started. This is a too long time when compared to most PNs. Accurate calculations by WA02 find this evolution time for the more massive cores, and a much longer evolution time from the AGB tip to the PN stage for lower mass stars. This is in contradiction with observations (R. Szczerba 2002, private communication). Note that Bowen & Willson (1991) get an ever increasing mass loss rate because they don’t treat properly the AGB and post-AGB evolution, as they don’t include the blueward stellar evolution part.

Although crude, the treatment above has a robust conclusion: Models which attribute the final intensive wind on the AGB and and post-AGB (the superwind) mainly to the decreasing effective temperature of the star (or its increasing radius), are in contradiction with observations of PNs and proto-PNs, because these models predict a too long post-AGB evolution.

## 3 THE DEPENDANCE OF MASS LOSS RATE ON THE ENVELOPE STRUCTURE

### 3.1 The density profile

Following SH99, we examine the ratio of the density at the photosphere, $`\rho _p`$, to the average envelope density $`\rho _a`$. 111Note that the density scale in Figs. 1-5 of SH99 is too low by a factor of 10; the correct scale is displayed in their Fig. 6.. The photospheric density is given by (Kippenhahn & Weigert 1990)

$$\rho _p=\frac{2}{3}\frac{\mu m_H}{k_B}\frac{GM_{}}{R^2\kappa T},$$
(5)
where $`\mu m_H`$ is the mean mass per particle, $`k_B`$ is the Boltzmann constant, and $`\kappa `$ is the opacity. Substituting typical values for a post-AGB star in the last equation gives

$$\rho _p=9.3\times 10^{10}\left(\frac{\kappa }{10^3\mathrm{cm}^2\mathrm{g}^1}\right)^1\left(\frac{T}{5000\mathrm{K}}\right)^3\left(\frac{L}{10^4L_{}}\right)^1\left(\frac{M_{}}{0.6M_{}}\right)\mathrm{g}\mathrm{cm}^3.$$
(6)
During the post-AGB phase, most of the envelope volume is convective. We are interested in the density in that region, hence we consider the mass of the convective region, $`M_{\mathrm{con}}`$. The average density in the convective region is given by

$$\rho _a=\frac{M_{\mathrm{con}}}{4\pi R^3/3}=6.0\times 10^8\left(\frac{T}{5000\mathrm{K}}\right)^6\left(\frac{L}{10^4L_{}}\right)^{3/2}\left(\frac{M_{\mathrm{con}}}{0.1M_{}}\right)\mathrm{g}\mathrm{cm}^3.$$
(7)
The ratio of photospheric to average density is then

$$\frac{\rho _p}{\rho _a}0.016\left(\frac{\kappa }{10^3\mathrm{cm}^2\mathrm{g}^1}\right)^1\left(\frac{T}{5000\mathrm{K}}\right)^3\left(\frac{L}{10^4L_{}}\right)^{1/2}\left(\frac{M_{}}{0.6M_{}}\right)\left(\frac{M_{\mathrm{con}}}{0.1M_{}}\right)^1.$$
(8)
At an effective temperature of $`T2800\mathrm{K}`$, appropriate for the tip of the AGB, the opacity is $`\kappa 3\times 10^4`$. We find from equation (8) that at the tip of the AGB, when most of the envelope is convective and $`M_{\mathrm{env}}=M_{\mathrm{con}}`$, that $`\rho _p/\rho _a0.3(M_{\mathrm{env}}/0.1M_{})^1`$. Of course, the average density must be larger than the photospheric density, hence the star starts to contracts, i.e. becomes hotter, before this stage is reached, i.e., when the envelope still has a significant amount of mass.

During the early post AGB phase, when $`T4000\mathrm{K}`$, the envelope shrink according to equation (2) and the opacity (e.g., Alexander & Ferguson 1994) can be approximated by $`\kappa 4\times 10^4(T/3000K)^4\mathrm{cm}^2\mathrm{g}^1`$. Using these in equation (8) with $`L=7500L_{}`$, gives for the density ratio (SH99)

$$\frac{\rho _p}{\rho _a}0.3\left(\frac{M_{\mathrm{con}}}{0.1M_{}}\right)^{0.3}0.4\left(\frac{T}{3000\mathrm{K}}\right)^3\mathrm{for}3000\mathrm{K}\mathrm{T}4000\mathrm{K}.$$
(9)
This ratio increases with mass loss, and implies that the density profile becomes shallower, almost flat, as indeed found in the numerical stellar models of SH99 for the early post-AGB phase. In SH99, the evolution of the envelope was followed at the end of the AGB phase, where the envelope mass decreased down to $`0.015M_{}`$ (where the mass of the stellar core was kept constant at $`0.6M_{}`$).

In the present work we follow the evolution during the early post-AGB. The numerical methods used here are the same as those used by SH99, but the present model has somewhat higher luminosity, $`L=9300L_{}`$ instead of $`L=8400L_{}`$ in SH99 or $`L=7500L_{}`$ in Soker (1992). To demonstrate the same behavior as found for late AGB stars for these later stages, we plot in Figure 1 the density profiles at several post-AGB evolutionary points, where the envelope mass dropped down to the values of $`0.0002M_{}`$. The envelope mass at the points displayed in the figure, has the values of $`0.0015,0.001,0.0005,0.0003`$, and $`0.0002M_{}`$. Similar figure is plotted for the entropy profile in the same models in Figure 2 (see next section).

An important point confronted in such evolutionary calculations, concerns the fact that due to the very low density in the outer envelope, the opacity in this region is very low, and the effective temperature of the star, calculated at the point where the optical depth, $`\tau `$, reaches the value of $`2/3`$, is deep in the outer envelope. In the figures, we display the run of the density and of the entropy out to the point where they actually do not play any role in the construction of the envelope structure. This point is about $`515R_{}`$ further than the point where the effective temperature was calculated.

An interesting feature that appears in the density run of all the models is the density inversion, that takes place in the outer part of the envelope. This phenomenon was discussed by Harpaz (1984). Close to the outer region where convection takes place, where the density is very low, the convective flux becomes very unefficient, because it depends on the density (see eq. 10 below). A very steep temperature gradient develops in order to increase the convective flux. However, this steep temperature gradient, causes an increase of the pressure, and since the pressure gradient have to satisfy the hydrostatic equation, a density inversion is developed, in order to compensate for the increase of the pressure due to the increase of the temperature. Certainly, the demand on the increased temperature gradient is removed when the opacity falls steeply at the point where the recombination of the hydrogen is completed, and there, the density inversion vanishes.

At the temperature range of $`5000\mathrm{K}\mathrm{T}6000\mathrm{K}`$, a steep increase of the opacity implies that the photospheric density steeply decreases with increasing effective temperature. Using the opacity table of Alexander & Ferguson (1994), we find from equation (6), with the same luminosity and stellar mass used there, that at $`T=4000\mathrm{K}`$, $`\kappa =10^3\mathrm{cm}^2\mathrm{g}^1`$ and $`\rho _p9\times 10^{10}\mathrm{g}\mathrm{cm}^3`$, at $`T=5000\mathrm{K}`$, $`\kappa =2.2\times 10^3\mathrm{cm}^2\mathrm{g}^1`$ and $`\rho _p4\times 10^{10}\mathrm{g}\mathrm{cm}^3`$, and at $`T=6000\mathrm{K}`$, $`\kappa =0.03\mathrm{cm}^2\mathrm{g}^1`$ and $`\rho _p=3\times 10^{11}\mathrm{g}\mathrm{cm}^3`$. This demonstrates the sharp drop in the photospheric density occurring as the star moves blueward across $`T5000\mathrm{K}`$. Comparing the change in photospheric density with the average density in the convective region (Fig. 1), shows that the photospheric density drops faster than the average density in the convective region, as the post-AGB star heats up in the temperature range of $`T50006000\mathrm{K}`$. Over all, in evolving from the AGB tip blueward to a temperature of $`T5000\mathrm{K}`$ the general envelope density profile becomes shallower, as evident from the increase of the ratio $`\rho _p/\rho _a`$. In the evolution from $`T5000\mathrm{K}`$ to $`T_f70009000\mathrm{K}`$, with higher temperatures for more luminous post-AGB stars, the trend of this behavior changes, but the density profile in the convective region stays shallow (Fig. 1). Becoming hotter than $`T_f`$, the star contracts very fast and the envelope density becomes very steep (Soker 1992). The range relevant to the present paper is the post-AGB track cooler than $`T_f`$.

If the shallow density profile is behind the high mass loss rate during the late AGB and early post-AGB phase (SH99; sec. 1 above), then it is expected that the mass loss rate will decrease significantly only as the star gets hotter than $`T=T_f70009000\mathrm{K}`$. Including the effect of temperature–low temperature is required to form dust–the sharp drop in the mass loss rate can occur somewhat earlier, at $`T6000\mathrm{K}`$.

### 3.2 The entropy profile

SH99 show that during the late-AGB and early post-AGB phases the entropy gradient becomes steeper with mass loss. In Figure 2 here we show that the entropy profile stays steep for most of the relevant post-AGB phase.

For stars on the upper AGB and early post-AGB phases the convective velocity, $`v_c`$, in the mixing length prescription is limited to not exceed the isothermal sound speed, and the pressure profile is relatively shallow (e.g., figs. 1-5 in SH99). Under these circumstances, the convective luminosity at each radius $`r`$ in the envelope can be written as

$$L_{\mathrm{conv}}(r)4\pi \alpha r^2T_e(r)v_c(r)\rho _e(r)H_P(r)|dS(r)/dr|,$$
(10)
where $`T_e`$ and $`\rho _e`$ are the temperature and density, respectively, inside the envelope (rather than on the photosphere), $`H_P`$ is the pressure scale height, and $`\alpha `$ is a parameter of order unity which is basically the ratio of the mixing length to $`H_P`$. The entropy $`S`$ is in units of $`\mathrm{erg}\mathrm{K}^1\mathrm{g}^1`$. The envelope temperature profile, $`T_e(r)`$, during the relevant post-AGB evolution segment and in the relevant envelope parts, evolves very slowly (e.g., Soker 1992). The same holds for $`H_P`$ and $`v_c`$ (SH99). Since the density in the envelope is very low, equation (10) implies that the entropy profile must be steep, as indeed we find from accurate numerical calculations (Fig. 2).

Again, If the steep entropy gradient is an important factor in creating the high mass loss rate during the late AGB and early post-AGB phase (SH99; sec. 1 above), then it is expected that the mass loss rate will decrease significantly only as the star gets hotter than $`T=T_f70009000\mathrm{K}`$ (higher temperatures for more luminous stars). As mentioned in the previous subsection, effects of dust formation due to effective temperature will cause a sharp drop in the mass loss rate already at somewhat lower temperatures.

## 4 SUMMARY

The main goal of the present paper is to argue that the decreasing temperature and increasing radius can’t be the only main parameters which determine the final intensive wind$``$ the superwind$``$on the tip of the AGB and the early post-AGB phases. The decreasing temperature and increasing radius, as well as the decrease in total stellar mass, lead to an increase in the mass loss rate of red giant branch and AGB stars. We accept this general trends. However, we showed (Sect. 2) that models with mass loss rates which depend mainly on these variables predict that the mass loss rate steeply decreases as the star leaves the AGB tip and becomes hotter (at constant luminosity). Because at that stage the envelope still contains a significant mass, the post-AGB evolution predicted by such models (e.g., WA02) is much too long. From observations and stellar evolution models fitted to observations, it is known that the final intensive wind should continue well into the post-AGB phase (Tylenda & Stasinska 1994), until the stellar effective temperature becomes $`T50006000\mathrm{K}`$, (higher temperature for more massive stars, e.g., Blöcker 1995).

Following SH99, we showed in this paper (Sect. 3) that the density profile stays shallow, and the entropy profile stays steep until the post-AGB stars heats up above $`T=T_f70009000`$, when the star shrinks fast (Soker 1992). Therefore, we suggest that the envelope structure, in particular the entropy and density gradients, should be among the main parameters which determines the mass loss rate on the tip of the AGB and the early post AGB. Some mechanisms which may influence the mass loss rate, and which are influenced by the envelope structure were mentioned in Section 1.

The sharp drop in the mass loss rate is likely to occur earlier, before the star heats-up to $`T=T_f`$, since it becomes more difficult to form dust as the star gets hotter. Over all, we suggest a general form for the mass loss rate on the tip of the AGB and early post-AGB phases

$$\dot{M}=\dot{M}_T[1+G(S^{},\rho ^{})],$$
(11)
where $`\dot{M}_T`$ is the mass loss rate mainly determined from the effects of the effective temperature, stellar radius and mass (e.g., eq. 1 above taken from WA02). $`G(S^{},\rho ^{})`$ is a function of the density and entropy gradients in the convective envelope, with $`G1`$ only when the density profile is shallow and the entropy profile is steep in the envelope convective region.

We thank Ryszard Szczerba for helpful discussions. This research was conducted in part while N.S. was visiting the Copernicus Astronomical Centre at Torun. This research was supported in part by grants from the US-Israel Binational Science Foundation and the Israel Science Foundation.