text:
# Phonon dispersion and lifetimes in MgB₂

## Abstract

We measure phonon dispersion and linewidth in a single crystal of $`\mathrm{MgB}_2`$ along the $`\mathrm{\Gamma }`$-A, $`\mathrm{\Gamma }`$-M and A-L directions using inelastic X-Ray scattering. We use Density Functional Theory to compute the effect of both electron-phonon coupling and anharmonicity on the linewidth, obtaining excellent agreement with experiment. Anomalous broadening of the $`E_{2g}`$ phonon mode is found all along $`\mathrm{\Gamma }`$-A. The dominant contribution to the linewidth is always the electron-phonon coupling.

The discovery of 39 K superconductivity in $`\mathrm{MgB}_2`$ Nagamatsu  has led to in-depth study of the material and a picture has emerged of a phonon-mediated superconductor with multiple gaps Shulga ; Liu ; Choi , moderate electron-phonon coupling (EPC) An ; Kortus ; Kong ; Liu , and anharmonicity Liu ; Yildirim ; Choi ; Boeri . However no measurements exist concerning either phonon dispersion or the evolution of phonon lifetimes over the Brillouin Zone (BZ), due to the absence of large single crystals. Neutron scattering on powder samples Osborn ; Yildirim  has been limited to the determination of phonon density of states. Raman spectroscopy, which is restricted to the BZ center, has shown that the optical mode with $`E_{2g}`$ symmetry, corresponding to in-plane distortions of the B hexagons, is strongly damped Bohnen ; Goncharov ; Postorino ; Hlinka .

Phonon damping can be caused by (i) EPC mediated phonon decay into electron-hole pairs Allen , or (ii) phonon-phonon interaction due to anharmonicity Cardona . The linewidth (the inverse of the lifetime) of a given phonon is the sum of both contributions. Direct determination of the contribution of each phonon mode to EPC Allen  from the measured linewidth is only possible if the anharmonic contribution is negligible Grimvall  and is seemingly questionable for $`\mathrm{MgB}_2`$ where many calculations Yildirim ; Choi ; Liu  suggest strong anharmonic effects.

In this work we present the first measured phonon dispersion curves and linewidths (where possible) along three directions in the BZ, $`\mathrm{\Gamma }`$-A, $`\mathrm{\Gamma }`$-M and A-L. We circumvented the problem of sample size by using high resolution inelastic scattering of a focused and intense X-Ray beam at the European Synchrotron Radiation Facility (beamline ID28), a technique Ruf  successfully used in single crystalline samples and in particular for the measurement of high energy optical modes Dastuto . To understand the mechanisms governing the measured phenomena, we calculated phonon dispersion, the contributions of EPC and anharmonicity to the linewidth and the structure factors using Density Functional Theory (DFT).

Small single crystals of $`\mathrm{MgB}_2`$, suitable for inelastic X-ray scattering experiments have recently become available. The crystal used in our experiment was grown at a pressure of 30-35 kbar. A mixture of Mg and B was put into a BN container in a cubic anvil device. The temperature was increased during one hour up to $`17001800^{}`$C, kept stable for 1-3 hours and decreased during 1-2 hours. As a result plate-like $`\mathrm{MgB}_2`$ crystals were formed of which we used a sample of about 400 x 470 x 40 $`\mu m^3,`$ with a measured in-plane mosaicity of $`0.007^{}`$. The beam incident on the sample was obtained from a high-resolution Silicon backscattering monochromator using the (8 8 8) reflection at an incident energy of 15.816 keV. The X-ray beam was focused onto the sample by a toroidal mirror into a spot of 270 x 90 $`\mu m^2`$(horizontal $`\times `$ vertical), full width at half maximum (FWHM). Slits before the sample further limited the vertical beam size to 30 $`\mu m`$. The scattered photons were analyzed in energy by five spherical silicon crystal analyzers operating at the same reflection order and mounted in pseudo Rowland circle geometry. The total energy resolution was 6.1 meV FWHM, as determined by a fit to a Lorentzian lineshape. The momentum transfer Q was selected by rotating the 7 m long analyzer arm around the sample position, in the horizontal plane, which also contained the linear x-ray polarization vector of the incident beam. The momentum resolution was set to $`0.04\mathrm{\AA }^1`$ in the horizontal direction and $`0.07\mathrm{\AA }^1`$ in the vertical direction. The following measurements were performed at a temperature of 300K: i) Q=(1 2 $`\xi `$), in almost transverse configuration along the $`\mathrm{\Gamma }`$-A direction, i.e. with $`𝐐𝐪0`$, q=(0 0 $`\xi `$) being the phonon wavevector, ii) Q=(1-$`\delta `$ 1+$`\xi `$ 0), while nearly following the $`\mathrm{\Gamma }`$-M direction ($`0\delta 0.05`$). iii) Q=(1-$`\delta _1`$ 2+$`\xi `$ 0.5+$`\delta _2`$), while nearly following the A-L direction ($`0\delta _1,\delta _20.04`$). The choice of BZ points measured was dependent on a series of conditions including the need to optimize the structure factors, single ($`\mathrm{\Gamma }`$-A) or multi-analyzer ($`\mathrm{\Gamma }`$-M, A-L) measurement mode and spectrometer and time limitations.

In Fig. 1 we show data taken at the 0.6 $`\mathrm{\Gamma }`$-A point in the BZ. The acoustic mode as well as the lower energy optical mode ($`E_{1u}`$) are visible as resolution-limited peaks. Most importantly, a broad peak is observed at higher energy loss, corresponding to the $`E_{2g}`$ mode. We performed least square fits to sums of Lorentzian functions with FWHM corresponding to the experimental resolution for the resolution limited peaks and a free parameter for the strongly damped phonon. These yield the dispersion as well as the linewidth variation over the BZ. Despite statistical limitations (3-6 counts per minute on this peak along $`\mathrm{\Gamma }`$-A) and tails of the peaks from the stronger, low energy phonons, the peak energy as well as the linewidth can be estimated with reasonable confidence.

Fig. 2 shows a similar energy loss scan at (0.97 2.29 0.54) close to the 0.58A-L point. A strong acoustic mode is seen at 30 meV. The peak at 50 meV corresponds to the $`A_{2u}`$ branch and the one at 65 meV to another acoustic branch. Finally two resolved features are seen at 85 and 97 meV. These are the two $`E_{2g}`$ modes which in this region of reciprocal space are well separated from other modes. Though, given the statistics, it would be hazardous to estimate a linewidth, the comparison between the experimental and ab-initio spectra suggests that the linewidths of both the $`E_{2g}`$ modes are resolution limited and so the damping is much less than that along $`\mathrm{\Gamma }`$-A. As for the measurement nearly along $`\mathrm{\Gamma }`$-M, the structure factor for the optical modes strengthens only near the zone boundary. At the point measured nearest to M ((1.05 1.45 0), not shown) the $`E_{2g}`$ and $`E_{1u}`$ modes are comparable in intensity but only separated by about 1 meV according to our calculation and we do measure a single peak only somewhat broader (FWHM $`10meV`$) than the experimental resolution, indicating reduced $`E_{2g}`$ linewidth. The proximity of the $`E_{1u}`$ mode however prevents a firm conclusion in this regard. We mention that the calculated structure factors and energies show excellent quantitative agreement with our measured data of which we have shown only two examples.

Similar analysis was done for several points along the three directions in order to experimentally determine the phonon dispersion and the linewidths. The difference in calculated phonon energies between the measured points and corresponding points exactly along $`\mathrm{\Gamma }`$-M and A-L ($`\delta ,\delta _1,\delta _2=0`$) is less than half a meV in all cases. We can thus compare the experimental phonon dispersion with the theoretical calculation along the high symmetry lines, as shown in the bottom panel of fig. 3 (circles).

The measured intrinsic linewidth of the $`E_{2g}`$ branch, shown in the top panel of fig. 3, is strongly anisotropic in the BZ. Along $`\mathrm{\Gamma }`$-A it is particularly large (ranging from 20 to 28 meV), while near L and probably near M it is below the experimental resolution.

Electronic structure calculations PWSCF  were performed using DFT in the generalized gradient approximation PBE . We used norm conserving pseudo-potentials Troullier . For $`\mathrm{Mg}`$ we used non-linear core corrections NLCC  and we treated the $`2s`$, $`2p`$ levels as core states. The wave-functions were expanded in plane waves using a $`35`$ Ry cutoff. The calculations were performed with the experimental crystal structure, namely $`a=3.083\AA `$ and $`c/a=1.142`$. The harmonic phonon frequencies were computed in the linear response DegironcPRB . We used a $`16\times 16\times 16`$ Monkhorst-Pack grid for the electronic BZ integration and first order Hermite-Gaussian smearing Degironcsmear  of $`0.025`$ Ry. The dynamical matrix at a given point of the BZ was obtained from a Fourier interpolation of the dynamical matrices computed on a $`6\times 6\times 4`$ phonon mesh. The resulting phonon frequencies are shown in fig. 3 and are in good agreement with a recent calculation Bohnen . The agreement with experiment is remarkable.

The contribution to the FWHM linewidth $`\gamma _{𝐪\nu }`$ at momentum $`𝐪`$ for the $`\nu `$ phonon mode due to the electron-phonon interaction can be written as Allen :

$$\gamma _{𝐪\nu }=\frac{4\pi \omega _𝐪}{N_k}\underset{𝐤,n,m}{}|g_{𝐤n,𝐤+𝐪m}^\nu |^2\delta (\epsilon _{𝐤n})\delta (\epsilon _{𝐤+𝐪m})$$
(1)
where the sum is extended over the BZ, $`N_k`$ is the number of $`k`$-points in the sum, and $`\epsilon _{𝐤n}`$ are the energy bands measured with respect to the Fermi level at point $`𝐤`$. The matrix element is $`g_{𝐤n,𝐤+𝐪m}^\nu =𝐤n|\delta V/\delta u_{𝐪\nu }|𝐤+𝐪m/\sqrt{2\omega _{𝐪\nu }}`$, where $`u_{𝐪\nu }`$ is the amplitude of the displacement of the phonon $`\nu `$ of wavevector $`𝐪`$, $`\omega _{𝐪\nu }`$ is the phonon frequency and $`V`$ is the Kohn-Sham potential.

In the calculations we used $`N_k=30^3`$ inequivalent $`k`$-points and, in Eqs. (1) and (3), we substituted the $`\delta `$ functions with Gaussians. The electron phonon coupling $`\lambda _{𝐪\nu }`$ is obtained from the linewidth Allen  as:

$$\lambda _{𝐪\nu }=\frac{\gamma _{𝐪\nu }}{2\pi N(0)\omega _{𝐪\nu }^2},$$
(2)
$`N(0)=0.354`$ states/($`\mathrm{MgB}_2`$ eV spin) being the density of states at the Fermi level.

The second contribution to the linewidth is given by the anharmonicity in the crystal potential. At lowest order for the mode $`\nu `$ of a zone center phonon the FWHM linewidth is Cardona ; Debernardi ; Lang :

$$\mathrm{\Gamma }_{\mathrm{𝟎}\nu }=\frac{\pi \mathrm{}}{8N_q}\underset{𝐪,\mu \eta }{}|\frac{^3E}{u_{\mathrm{𝟎}\nu }u_{𝐪\mu }u_{𝐪\eta }}|^2\frac{I_{𝐪\mu \eta \nu }^\mathrm{D}+I_{𝐪\mu \eta \nu }^\mathrm{A}}{\omega _{\mathrm{𝟎}\nu }\omega _{𝐪\mu }\omega _{𝐪\eta }}$$
(3)
$`E`$ being the total energy, $`n_{𝐪\mu }`$ the Bose occupation for mode $`\mu `$ at wavevector $`𝐪`$, $`I_{𝐪\mu \eta \nu }^\mathrm{D}=(n_{𝐪\mu }+n_{𝐪\eta }+1)\delta \left(\omega _{\mathrm{𝟎}\nu }\omega _{𝐪\mu }\omega _{𝐪\eta }\right)`$ describes the decay in the two phonons $`\mu `$ and $`\eta `$, and $`I_{𝐪\mu \eta \nu }^\mathrm{A}=2(n_{𝐪\mu }n_{𝐪\eta })\delta \left(\omega _{\mathrm{𝟎}\nu }\omega _{𝐪\mu }+\omega _{𝐪\eta }\right)`$ describes the $`\eta `$-phonon absorption and the $`\mu `$-phonon emission.

We computed the anharmonic linewidth at the high-symmetry points $`\mathrm{\Gamma }`$, A, M. For the calculation at A we consider a $`1\times 1\times 2`$ supercell with 6 atoms, while for the M point we use a $`2\times 2\times 1`$ cell with 12 atoms. The third order matrices were evaluated using linear response theory and the $`2n+1`$ theorem for metals Lazzeri . The anharmonic contribution was evaluated at $`0K`$ and $`300K`$.

At $`\mathrm{\Gamma }`$ the anharmonic linewidth is largest for the $`E_{2g}`$ mode and equal to $`0.16`$ meV at $`T=0K`$ and $`1.21`$ meV at $`T=300K`$. Both the values are negligible if compared with the experimental Raman linewidth of roughly $`40`$ meV Postorino , suggesting that the main source of broadening is the electron-phonon interaction.

The results of the calculation of the two contributions to the linewidth at A and M are shown in Table 1. At the A point the linewidth of the $`E_{2g}`$ mode due to electron-phonon scattering is very large while the anharmonic contribution is more than an order of magnitude smaller. For sizable values of the EPC and at low temperature the anharmonicity is negligible showing that a measurement of the linewidth at $`𝐪`$ for the mode $`\nu `$ is essentially equivalent to determining $`\lambda _{𝐪\nu }`$. This is unexpected since earlier theoretical work Yildirim ; Choi ; Liu  estimates anharmonicity to be important for a different but related quantity, the frequency shift.

Since the anharmonic contribution is negligible at $`\mathrm{\Gamma }`$, A, M, and computationally demanding, we only evaluated the EPC linewidth for the other points along the three directions. Using the calculated phonon frequencies, displacements and linewidths we compute the structure factors for one phonon processes using the X-ray form factors, obtaining good agreement with experiment as shown in Figs. 1,2. In the top panel of Fig. 3 we show the theoretical results for the $`E_{2g}`$ branch. The result is consistent with the experimental value wherever the $`E_{2g}`$ branch is clearly visible namely along $`\mathrm{\Gamma }`$-A and near the L point. Finally, the anharmonicity being negligible, we used Eq. (2) to extract $`\lambda _{𝐪\nu }`$ of the $`E_{2g}`$ mode along $`\mathrm{\Gamma }`$-A using the measured linewidths and frequencies together with the calculated electronic density of states. In Table 2 the experimental values are compared with the theoretical predictions. The anomalously large EPC along $`\mathrm{\Gamma }`$-A is due to the nesting factor of the B bonding $`p_{x,y}`$ Fermi surfaces, which are concentric cylinders centered on $`\mathrm{\Gamma }`$-A An ; Kortus . The $`E_{2g}`$ modes, which modify the B-B distances, are the only ones with a sizeable matrix element, $`g_{𝐤n,𝐤+𝐪m}^\nu `$, between electrons on these surfaces.

In conclusion, we have measured phonon dispersion and linewidths in a sub-mm sized $`\mathrm{MgB}_2`$ crystal with inelastic X-ray scattering confirming the power and versatility of this technique. Both acoustic and optical modes are detected and we find that the $`E_{2g}`$ mode is anomalously broadened along $`\mathrm{\Gamma }`$-A but that this broadening is not generalized over the Brillouin Zone. Our Density Functional Theory calculations of the dispersion and linewidth are in excellent agreement with experiment. They show that the dominant contribution to the broadening for all modes is the electron-phonon coupling, the anharmonic contribution being much smaller. Thus phonon linewidth in $`\mathrm{MgB}_2`$ is a direct measure of electron-phonon coupling and could, with the availability of larger samples, be measured for all modes over the whole Brillouin zone so as to extract the anisotropic Eliashberg coupling function.

We acknowledge illuminating discussions with R. S. Gonnelli, P. Giannozzi, M. Xu, and F. Sette. The calculations were performed at the IDRIS supercomputing center. M.C. was supported by a Marie Curie Fellowship of the European Commission, contract No. IHP-HPMF-CT-2001-01185.


format: text

text:
# Light quark masses from UKQCD’s dynamical simulations with O(a)-improved Wilson fermions

## 1 Introduction

The determination of the masses of the light quarks represents an important challenge in particle physics today. They are among the least well determined parameters of the Standard Model, having uncertainties of the same order as their values as given in the particle data book.

Lattice QCD is in a position to improve this situation because it allows one to study the quark mass dependence of physical observables such as hadron masses. The physical values of these observables can then be used to fix the quark masses. Unfortunately current simulations cannot probe the up/down quark mass range directly and so some form of extrapolation is necessary. Chiral perturbation theory is at present the best procedure for providing a link to the low mass regime. The simulations studied in this paper however are still too heavy for $`\chi ^{PT}`$ and a simpler extrapolation ansatz is used instead.

In this paper I present results for the average up/down quark mass, $`m_l`$, and the strange quark mass, $`m_s`$, from a partially quenched analysis of pseudoscalar and vector meson masses for each of UKQCD’s dynamical fermion datasets. The suitably averaged values of the hadron masses used to fix the quark masses are as follows: $`M_\pi =137.3\text{MeV}`$ ; $`M_K=495.7\text{MeV}`$ ; $`M_K^{}=893.9\text{MeV}`$ and $`M_\varphi =1020\text{MeV}`$.

## 2 Simulation Details

All of the lattice simulations studied here have been carried out on a $`16^3\times 32`$ lattice with a non-perturbatively O(a)-improved Wilson fermion action and standard Wilson glue. Details of these simulations have been discussed previously in . Additionally the set of $`\kappa _{\mathrm{val}}`$ values for the datasets at $`\beta =5.2,\kappa _{\mathrm{sea}}=0.13565`$ and $`\beta =5.2,\kappa _{\mathrm{sea}}=0.1358`$ have been increased to $`\{0.1342,0.1350,0.1355,0.13565,0.13580\}`$. Pseudoscalar and vector meson correlation functions have been generated with all non-degenerate combinations of these $`\kappa _{\mathrm{val}}`$ values. For the $`\kappa _{\mathrm{sea}}=0.1358`$ dataset these correlators have also been calculated at $`\kappa _{\mathrm{val}}=0.13595`$.

## 3 Hadron masses

The procedures used here to extract hadron masses from correlation functions have been described in detail in . In brief the masses of the pseudoscalar and vector mesons were determined from double cosh fits to their corresponding correlation functions. Local and fuzzed correlators were fitted simultaneously to constrain fit parameters and a sliding window analysis used to find an appropriate time window for each fit.

## 4 Chiral Extrapolations

The pseudoscalar meson masses for each dataset were extrapolated in $`1/\kappa _{\mathrm{val}}`$ using the simple ansatz

$`M_{PS}^2(\kappa _1,\kappa _2)`$ $`=`$ $`A_{PS}(m_1+m_2)`$ (1)
$`=`$ $`A_{PS}({\displaystyle \frac{1}{2}}({\displaystyle \frac{1}{\kappa _1}}+{\displaystyle \frac{1}{\kappa _2}}){\displaystyle \frac{1}{\kappa _{\mathrm{crit}}}})`$ (2)
The value of $`\kappa _{\mathrm{crit}}`$ was determined as a free parameter in this fit.

For the vector meson the ansatz

$$M_V(m_1,m_2)=A_V(m_1+m_2)+B_V$$
(3)
was found to fit the data well.

An example plot of these chiral extrapolations is given in Figure 1 for the $`\kappa _{\mathrm{sea}}=0.1358`$ dataset. The lightest pseudoscalar and vector points at $`\kappa _{\mathrm{val}}=0.13595`$ are included in the plot but were not included in the extrapolation. The lightest vector point clearly deviates from the linear behaviour suggesting large finite size effects. For more discussion of finite size effects in these datasets I refer the reader to .

## 5 Fixing the bare light quark masses

Determination of the light quark masses from the results of the chiral extrapolations requires setting the scale by some physical quantity. In this work the Sommer scale parameter, $`r_0=0.49\text{fm}`$, and $`M_\rho =769.3\text{MeV}`$ were used.

The average up/down quark mass, $`m_l`$, was determined by requiring that the pseudoscalar mass equal the physical pion mass at that point i.e.

$$M_{PS}^2(m_l,m_l)=aM_\pi $$
(4)
The strange quark mass, $`m_s`$, was then determined from each of $`M_K`$, $`M_K^{}`$ and $`M_\varphi `$ by solving

$$M_{PS}^2(m_l,m_s)=aM_K$$
(5)
$$M_V(m_l,m_s)=aM_K^{}$$
(6)
$$M_V(m_s,m_s)=aM_\varphi $$
(7)
When the lattice spacing was fixed by $`r_0`$ however only the $`M_K`$ method gave a sensible value for $`m_s`$ and only these results are given.

## 6 Renormalisation and matching

In this work the values for the renormalised quark masses are given in the $`\overline{MS}`$ scheme at a reference scale of $`2\text{GeV}`$. The relationship between the renormalised and bare quark masses is given by the relation

$$m_R^{\overline{MS}}(\mu )=Z_m^{\overline{MS}}(\mu )(1+b_mam_q)m_q$$
(8)
where the matching has been performed at a scale, $`\mu `$. Running the masses from the scale, $`\mu `$, to a scale of $`2\text{GeV}`$ was achieved using the Mathematica package, RunDec , which implements 4-loop running from perturbation theory.

One-loop tadpole-improved perturbative results  were used for $`Z_m`$ and $`b_m`$ reorganised in terms of $`\alpha ^{\overline{MS}}(\mu )`$

$$Z_m^{\overline{MS}}(\mu )=[1\frac{\alpha ^{\overline{MS}}(\mu )}{4\pi }(8\mathrm{ln}(a\mu )15.085)]u_0^1$$
(9)
$$b_m=[1/20.7850\alpha ^{\overline{MS}}(\mu )]u_0^1$$
(10)
The tadpole factor, $`u_0`$, was taken to be $`1/8\kappa _c`$ here.

An estimation of the best value for the matching scale, $`\mu =q^{}`$ , has not been determined for this work. Instead the cases $`\mu =1/a`$ and $`\mu =\pi /a`$ have been investigated and the variation treated as a systematic effect.

The values of $`\alpha ^{\overline{MS}}(\mu )`$ were also determined with the RunDec package using values for $`\mathrm{\Lambda }_{\overline{MS}}`$ given in .

## 7 Results

Results for $`m_l`$ and $`m_s`$ have been plotted in Figures 2 and 3 respectively. Datasets matched in $`r_0`$ have been distinguished from the others and a matched quenched dataset included for comparison. Only results for $`\mu =1/a`$ are plotted as choice of matching scale was found to have around a $`1\text{\%}`$ effect. The results for $`m_s(M_\varphi )`$ are omitted as they were almost identical to those for $`m_s(M_K^{})`$.

A large dependence on the choice of scale setting parameter can be seen for both $`m_l`$ and $`m_s`$. It should be noted however that a smaller physical value for $`r_0`$ would bring the results into better agreement. There is also a clear difference between $`m_s(M_K)`$ and $`m_s(M_K^{})`$ though this would be expected to vanish in the continuum limit.

There is some indication of a trend in the matched datasets as the sea quark mass is decreased, as expected. When compared to the matched quenched point however there is no clear signal of unquenching from either quark mass though this may improve with statistics.

The ratio, $`m_s/m_l`$, was calculated for each scheme and found to be between $`25.0`$ and $`25.5`$ for $`m_s(M_K)`$ and around $`31.0`$ for $`m_s(M_K^{})`$. This can be compared with the result from $`\chi ^{PT}`$  of $`24.4(1.5)`$.

## 8 Acknowledgments

I would like to thank the Carnegie Trust for the Universities of Scotland, the European Community’s Human Potential programme under HPRN-CT-2000-00145 Hadrons/LatticeQCD and PPARC under grant PP/G/S/1998/00777 for their financial support of this work. I also thank Craig McNeile for his invaluable help.


text:
# Efficient radiative transfer in dust grain mixtures

## 1 Introduction

The simulation of spectral energy distributions (SEDs), images, and polarization maps of young stellar objects has become a profound basis for the analysis and interpretation of observing results. Many techniques and approximations for the solution of the radiative transfer (RT) problem in different model geometries (1D–3D), considering more and more special physical processes, such as the stochastic and photo-electric heating of small grains (see, e.g, Draine & Li 2001; Bakes & Tielens 1994, Siebenmorgen, Krügel, & Mathis 1992), scattering by spheroidal grains (see, e.g., Wolf, Voshchinnikov, & Henning 2002; Gledhill & McCall 2000), or the coupling of line and continuum RT (see, e.g., Rybicki & Hummer 1992), have been developed. The simulation of the temperature structure in simple-structured circumstellar shells or disks (see, e.g., Malbet, Lachaume, & Monin 2001; Chiang et al. 2001), the estimation of the properties (luminosity, temperature, mass) of heavily embedded stars (see, e.g., Kraemer et al. 2001), or the determination of the inclination of a circumstellar disk (see, e.g., Chiang & Goldreich 1999; Men’shchikov, Henning, & Fischer 1999; Wood et al. 1998) represented modest, first attempts of the application of the existing sophisticated numerical techniques. More recent efforts are directed to derive the dust grain size distribution in a circumstellar disk from its SED (Wood et al. 2002; D’Alessio, Calvet & Hartman 2001). Thus, it is clear that beside strong observational constrains, the model parameters and considered physical processes have to be questioned in depth in order to derive such detailed conclusions. However, looking behind the scenes, many of the RT models are based on simplifying assumptions of very basic processes such as isotropic instead of anisotropic scattering, mean dust parameters representing dust grain ensembles (different radii and chemical compositions), or the flux-limited diffusion approximation \- approximations which are well-suited for handling the energy transfer in hydrodynamic simulations or a rough data analysis but which may not necessarily guarantee the desired accuracy for a detailed SED and image/polarization data analysis.

In the advent of (space) observatories such as SIRTF111Space Infrared Telescope Facility, which will be able to obtain SEDs of evolved debris disks around young stars, there exists a strong need for adequate numerical RT techniques in order to allow to trace dust grain growth (Meyer et al. 2001) and the influence of other physical effects and processes such as the Poynting-Robertson effect (see, e.g., Srikanth 1999) and dust settling (Dubrulle, Morfill, & Sterzik 1995, Miyake & Nakagawa 1995) on the dust grain size distribution in these disks. Therefore, we clearly have to understand which influence the different approximations (as long as they are required) in the RT simulations may have on the resulting observables. Based on two different grain size distributions consisting of astronomical silicate, the differences in the resulting dust grain temperature distributions and the resulting SEDs between simulations of the RT in (A) a “real” dust grain mixture and (B) under the assumption (approximation) of mean dust grain parameters will be discussed. This investigation is therefore mainly focused on two questions: (1) Of what order of magnitude are the differences (?) and (2) How many grain sizes have to be considered to represent the properties of a real grain size distribution? The spatial temperature distribution of the considered dust configurations is calculated on the basis of local thermal equilibrium. Stochastic heating processes which are expected in case of very small grains consisting of tens to hundreds of atoms (see, e.g., Draine & Li 2000 and references therein) are not subject of this investigation.

In Sect. 2, the RT and the dust grain model are briefly introduced. In Sect. 3, the definition of the mean dust grain parameters is given and the expected deviations of the RT results – once based on the mean dust parameters and once on a real grain size distribution – are outlined. In Sect. 4, the RT in a spherical shell with variable optical depth and density distribution (see Sect. 4.1) is considered, while the temperature structure in a model of the HH 30 circumstellar disk is investigated in Sect. 4.2. The SED resulting from a models with a dust grain mixture is compared to the mean particle approximation in Sect. 5.

## 2 Radiative transfer and Dust model

The RT simulations presented in this article have been performed with the three-dimensional continuum RT code MC3D which has been described by Wolf & Henning 2000 (see also Wolf, Henning, & Stecklum 1999; Wolf 2002). It is based on the Monte-Carlo method and solves the RT problem self-consistently. Instead of the iterative procedure of estimating the dust temperature distribution (as described in Wolf et al. 1999), the concept presented by Bjorkman & Wood (2001), which is based on the immediate correction of the dust grain temperature after absorption of a photon package, has been used. Furthermore, the method described by Lucy (1999), which takes into account the absorption not only at the points of interaction of the photons with the dust grains but also in-between, has been applied. While the first method allows to simulate the RT for optical depths222Throughout the whole article the quantity $`\tau `$ refers to the optical depth based on the extinction cross section of the dust grains at a wavelength of 550 nm. $`\tau 10^3`$, the second was used in order to increase the efficiency of the simulation for low optical depths. In addition to the test of MC3D described by Wolf et al. (1999) considering a single grain size, it has been successfully tested for the case of dust grain mixtures against the one-dimensional code of Chini et al. (1986).

The RT is simulated at separate wavelengths within the wavelength range \[$`\lambda _{\mathrm{min}}`$, $`\lambda _{\mathrm{max}}`$\]. For this reason, the radiation energy of the emitting source(s) is partitioned into so-called weighted photons (Fischer, Henning, & Yorke 1999) each of which is characterized by its wavelength $`\lambda `$ and Stokes vector $`\widehat{I}=(I,Q,U,V)^\mathrm{T}`$. In order to consider dust grain ensembles instead of a single grain species, the RT concept as described by Wolf et al. (1999), Lucy (1999), and Bjorkman & Wood (2001) has to be extended in the following manner:

1. The mean free path length $`l`$ between the point of emission and the point of the first photon-dust interaction (and - subsequently - between two points of interaction) is given by

$$\tau _{\mathrm{ext},\mathrm{l}}=\mathrm{ln}(1Z),$$
(1)
$$\stackrel{~}{\tau }_{\mathrm{ext},\mathrm{l}}=\underset{i=1}{\overset{n_\mathrm{P}}{}}\left[\underset{j=1}{\overset{n_\mathrm{D}}{}}\underset{k=1}{\overset{n_\mathrm{S}}{}}\rho _{j,k}(\stackrel{}{r}_{\mathrm{l}_\mathrm{i}})C_{\mathrm{ext}_{\mathrm{j},\mathrm{k}}}\right]\mathrm{\Delta }l_i,l=\underset{i=1}{\overset{n_\mathrm{P}}{}}\mathrm{\Delta }l_i.$$
(2)
Here, $`\tau _{\mathrm{ext},\mathrm{l}}`$ is the optical depth along the path of the length $`l`$, $`\stackrel{~}{\tau }_{\mathrm{ext},\mathrm{l}}`$ is the corresponding numerical accumulated optical depth ($`\tau _{\mathrm{ext},\mathrm{l}}\stackrel{~}{\tau }_{\mathrm{ext},\mathrm{l}}0`$), $`n_\mathrm{D}`$ is the number of chemically different dust grain species, $`n_\mathrm{S}`$ is the number of particle radii being considered, and $`\stackrel{}{r}_{\mathrm{l}_\mathrm{i}}`$ is the spatial coordinate corresponding to the $`i`$th integration point along the path length $`l`$. The quantity $`\rho (\stackrel{}{r})`$ is the number density at the spatial coordinate $`\stackrel{}{r}`$, the quantity $`C_{\mathrm{ext}}`$ is the extinction cross section, and $`n_\mathrm{P}`$ is the number of integration points along the path. Furthermore, $`Z`$ is a random number uniformly distributed in the interval $`[0,1]`$.
2. According to the concept of immediate reemission described by Bjorkman & Wood (2001), at each point of interaction of a photon with a dust grain either scattering or absorption occurs. The probability for a photon to undergo the one or the other interaction process with a grain with the chemical composition #$`j`$ and the particle radius #$`k`$ is given by

$$\mathrm{\Pi }_{x_{j,k}}(\stackrel{}{r})=\frac{\rho _{j,k}(\stackrel{}{r})C_{x_{j,k}}}{_{m=1}^{n_\mathrm{D}}_{p=1}^{n_\mathrm{S}}\left[\rho _{m,p}(\stackrel{}{r})C_{ext_{m,p}}\right]},$$
(3)
where ’$`x`$’ stands either for absorption or scattering. In case of absorption, the immediate reemission occurs from the same dust grain (species).
3. In Lucy’s concept (1999), which considers the absorption of the electromagnetic radiation field not only at the end points of the photon path (points of interaction) but also in-between, simply the absorption due to all dust grain species (instead of only one) has to be taken into account:

$$\tau _{\mathrm{abs}}=\underset{j=1}{\overset{n_\mathrm{D}}{}}\underset{k=1}{\overset{n_\mathrm{S}}{}}_{\mathrm{point}_1}^{\mathrm{point}_2}C_{\mathrm{abs}_{\mathrm{j},\mathrm{k}}}\rho _{j,k}(\stackrel{}{r})𝑑r.$$
(4)

The investigations presented in this paper are based on dust grain ensembles with different grain size distributions but the same chemical composition. Since the formalism (and therefore at least the qualitative conclusions) are the same for ensembles consisting of grains with different chemical composition or size distribution or even both, this restriction is justified. The optical properties of astronomical silicate (Draine & Lee 1984) for \[A\] a size distribution of small grains with radii $`a=(0.0050.25)\mu `$m and \[B\] a size distribution of larger grains with radii $`a=(0.11)\mu `$m have been used in Sect. 3, 4.1, and 5.

The grains are assumed to be spherical with a size distribution described by the widely applied power law $`n(a)a^{3.5}`$ (Mathis, Rumpl, & Nordsieck 1977). The optical properties such as the extinction and scattering cross sections as well as the scattering distribution function have been derived on the basis of Mie scattering, using the Mie scattering algorithm published by Bohren & Huffman (1983). The correct (anisotropic) scattering distribution function (for each dust grain size) has been considered in the RT process.

## 3 Mean dust grain parameters

In case of non-self-consistent RT in a dust grain mixture, the numerical effort (finally, the computing time and RAM requirement) may be substantially decreased by assuming weighted mean values of those parameters which describe the interaction of the electromagnetic field with the dust grains: the extinction, absorption, and scattering cross section, the Stokes parameters, and – as as function of the extinction and scattering cross section – the Albedo.

The weight $`w_j(a)`$, which represents the contribution of the $`j`$th component of the dust grain mixture, results from the abundance of this component in respect of its dust grain number density (assuming $`n_\mathrm{D}`$ chemically different dust species) and the size distribution of the respective material:

$$\underset{j=1}{\overset{n_\mathrm{D}}{}}_{a_{\mathrm{min}}}^{a_{\mathrm{max}}}w_j(a)𝑑a=1,$$
(5)
where $`a_{\mathrm{min}}`$ and $`a_{\mathrm{max}}`$ are the minimum and maximum grain radius of the size distribution. The Stokes parameters as well as the extinction, absorption, and scattering cross section ($`C_{\mathrm{ext}}`$, $`C_{\mathrm{abs}}`$, $`C_{\mathrm{sca}}`$) are additive. Therefore, the representative values in case of a dust grain mixture can be derived on the basis of their weighted contributions (see, e.g., Martin 1978, S̆olc 1980):

$$C_{\mathrm{ext}}=\underset{j=1}{\overset{n_\mathrm{D}}{}}_{a_{\mathrm{min}}}^{a_{\mathrm{max}}}w_j(a)C_{\mathrm{ext}_j}(a)𝑑a,$$
(6)
$$C_{\mathrm{abs}}=\underset{j=1}{\overset{n_\mathrm{D}}{}}_{a_{\mathrm{min}}}^{a_{\mathrm{max}}}w_j(a)C_{\mathrm{abs}_j}(a)𝑑a,$$
(7)
$$C_{\mathrm{sca}}=\underset{j=1}{\overset{n_\mathrm{D}}{}}_{a_{\mathrm{min}}}^{a_{\mathrm{max}}}w_j(a)C_{\mathrm{sca}_j}(a)𝑑a,$$
(8)
and

$$\widehat{S}=\underset{j=1}{\overset{n_\mathrm{D}}{}}_{a_{\mathrm{min}}}^{a_{\mathrm{max}}}w_j(a)\widehat{S}_j(a)𝑑a,$$
(9)
where $`\widehat{S}`$ is the Mueller matrix which is used to describe the modification of the photon’s Stokes vector due to the interaction of a photon with the scattering/absorbing medium (dust grains; see Bickel & Bailey 1985, Bohren & Huffman 1983). For the albedo $`A`$ follows:

$$A=\frac{_{j=1}^{n_\mathrm{D}}_{a_{\mathrm{min}}}^{a_{\mathrm{max}}}w_j(a)C_{\mathrm{ext}_j}(a)A_j(a)𝑑a}{_{j=1}^{n_\mathrm{D}}_{a_{\mathrm{min}}}^{a_{\mathrm{max}}}w_j(a)C_{\mathrm{ext}_j}(a)𝑑a}=\frac{C_{\mathrm{sca}}}{C_{\mathrm{ext}}}$$
(10)
Since this formalism (Eq. 5-10) covers both grain size and chemical distributions, the conclusions to be derived for a grain size distribution but a single chemical component is also valid for grain ensembles consisting of grains with different chemical compositions. The mean dust grain parameters used in the following sections have been derived on the assumption of $`10^3`$ grain radii equidistantly distributed in the grain size range of the small/large grain ensemble.

In thermal equilibrium, the temperature of a dust grain ($`u`$th grain size or composition) can be estimated from the assumption of local energy conservation. If

$$dE_u^{\mathrm{abs}}=dt_0^{\mathrm{}}L_{\lambda ,u}^{\mathrm{abs}}𝑑\lambda $$
(11)
is the energy being absorbed by the $`u`$th dust grain component during the time interval $`\mathrm{d}t`$, then the local energy conservation can be written as

$$dE_u^{\mathrm{ree}}=dt_0^{\mathrm{}}L_{\lambda ,u}^{\mathrm{ree}}𝑑\lambda =dt_0^{\mathrm{}}4\pi a^2Q_{\lambda ,u}^{\mathrm{abs}}(a)\pi B_\lambda (T_u)𝑑\lambda =dE_u^{\mathrm{abs}},$$
(12)
where $`dE_u^{\mathrm{ree}}`$ is the total amount of energy being re-emitted during the time interval $`dt`$, $`L_{\lambda ,u}^{\mathrm{abs}}`$ and $`L_{\lambda ,u}^{\mathrm{ree}}`$ are the monochromatic absorbed/re-emitted luminosity at the wavelength $`\lambda `$, $`Q_{\lambda ,u}^{\mathrm{abs}}(a)`$ is the absorption efficiency and $`B_\lambda (T_u)`$ is the Planck function of the grain with the temperature $`T_u`$. For two particles of the same composition but different radii ($`a_1`$, $`a_2`$) the temperature is consequently only the same if

$$Q_\lambda ^{\mathrm{abs}}(a_1)=Q_\lambda ^{\mathrm{abs}}(a_2).$$
(13)
This is usually not the case. Thus, mean dust grain parameters may be used only in case of the simulation of dust scattering but not for the estimation of the dust grain temperature distribution and resulting observables, such as the spectral energy distribution and images in the mid-infrared to millimeter wavelength range.

In the following Sections (§4 and §5) it will be investigated under which conditions the approximation of the mean dust grain parameters provides nevertheless a good estimate for these observables. While the most crucial information can be derived on the basis of a one-dimensional model (§4), the study presented in §5 is aimed to reveal the influence of geometrical effects in case of two-dimensional disk-like structures. In each case, typical parameters for both the dust grain composition (see §2) and the circumstellar shells/disks are considered.

## 4 Temperature Distribution

### 4.1 Basic Model: Spherical Shell

In order to reveal the differences between the results of the RT based on real dust grain size distributions and the approximation of mean dust grain parameters representing the weighted mean optical parameters of the grain size distribution, a simple-structured but in respect of the astrophysical applications very useful model has been used for the following investigations. It consists of a spherical shell (outer radius: $`10^3`$ AU; inner boundary: $`0.2`$ AU) with a radial density profile described by a power law, $`\rho (r)r^\alpha `$, where $`\alpha `$ is a negative, constant quantity. The dust grain size distribution is assumed to be constant throughout the whole shell. The only heating source is an embedded, isotropically radiating star in the center of the shell with an effective temperature and radius identical to that of the Sun.

In order to study the influence of multiple scattering and (additional) heating due to the dust reemission, the optical depth of the shell $`\tau `$ (as seen from the star) and the exponent $`\alpha `$, describing the relative density gradient in the shell, have been considered as variable parameters of the model. Four values (-1.0, -1.5, -2.0, and -2.5) for the exponent $`\alpha `$ were chosen, covering a broad range of astrophysical objects ranging from circumstellar shells (see, e.g., Chini et al. 1986, Henning 1985, Yorke 1980) and disks (see, e.g., Men’shchikov & Henning 1997, Sonnhalter et al. 1995, Lopez et al. 1995) to the radial density profile being assumed for the dust density distribution in AGN models (see, e.g., Manske et al. 1998, Efstathiou & Rowan-Robinson 1995, Stenholm 1994). The optical depth $`\tau `$ was varied from the optical thin case ($`\tau `$=0.1) to the case of intermediate optical depth ($`\tau =1`$) and finally to $`\tau =10`$ as the optically thick case. For optical depths below this interval no remarkable differences to the case $`\tau =0.1`$ are expected since multiple scattering and reemission have a negligible influence on the dust grain temperature which is determined only by the distance from the star (attenuation of the stellar radiation field) and the absorption efficiency of the dust grains. For optical depths $`\tau 10`$ preparatory studies have shown that the temperature difference between the different dust grains becomes negligible, converging towards the temperature distribution obtained on the basis of the mean dust grain approximation. This finding is in agreement with results obtained by Krügel & Walmsley (1984 – investigation of dust and gas temperatures in dense molecular cores), Wolfire & Churchwell (1987 – study of circumstellar shells around low-mass, pre-main-sequence stars), and Efstathiou & Rowan-Robinson (1994 – multi-grain dust cloud models of compact H II regions). Furthermore, only the inner region of the circumstellar shell shows a broad distribution of dust grain temperatures at a given distance from the star.

In Fig. 1 and 2 the differences between the radial temperature distribution of the single dust grains $`T(r,a)`$ and that resulting from the approximation of mean dust grain parameters $`\overline{T}(r)`$ is shown from the inner boundary of the shell to a distance of 1 AU from the star. The main results are:

1. The temperature of the different grains spans a range of $`250`$ K ($`25`$% of the corresponding mean temperature) around the temperature of the mean grains at the inner boundary of the shell.
2. The temperature difference decreases

* towards larger distances from the star, and
    * with increasing optical depth as soon as the shell becomes optically thick.
3. The relativ density gradient, described by the exponent $`\alpha `$ is of minor importance. However, an increase of the absolute amount of this exponent results in an increased redistribution of energy between the different dust grain components at the inner boundary of the shell ($`\tau 1`$) and therefore in a decrease of the temperature dispersion.
4. The difference ($`T(r,a)\overline{T}(r)`$) strongly depends on both the absorption efficiency, and thus on the dust grain size distribution, and the temperature of the heating source since the combination of these parameters determines the amount of energy being absorbed. This is clearly illustrated by the different signs of the temperature difference of the large and the small grain size distribution to the temperature of the mean grains shown in Fig. 1 and 2.

The found complex, highly dispersed temperature structure in a dust grain mixture requires a large number (here: $`2^6`$) of considered grain radii to be correctly taken into account and cannot be adequately represented by a mean particle. This finding is of special interest for the simulation of chemical networks, since the existence of ice layers, the possibility of a certain reaction on the grain surface, surface reaction rates, and the temperature of the surrounding gas phase depend on the dust grain temperature (see, e.g., Cazaux & Tielens 2002, Markwick et al. 2002, Charnley 2001). The dust phase in the close stellar environment is of particular importance, since the chemical evolution takes place here on its smallest timescale. Furthermore it was found, that the approximation of mean particle properties is found to be justified for dust configurations of high ($`\tau 10`$) optical depth.

### 4.2 HH 30 circumstellar disk

Although all basic characteristics of the temperature structure in a dust grain mixture can be studied on the basis of a one-dimensional model (Sect. 4.1), a disk-like structure will be considered in the following. The motivation for this is to investigate the possible influence of geometrical effects which have been minimized in case of the spherical shell. Furthermore, two-dimensional models are of particular importance since they are widely applied in simulations of circumstellar disks, debris disks around evolved stars, active galactic nuclei, and galaxies (see, e.g., Chiang & Goldreich 1997, 1999; see also the references in Sect. 4.1).

In the following, a model of the nearly edge-on disk around the classical T-Tauri star HH 30 is considered (see, e.g., Burrows et al. 1996 for observational details). Following Cotera et al. (2001) and Wood et a. (2002), a flared geometry as described by Shakura & Sunyaev (1973) is adapted:

$$\rho =\rho _0\left(\frac{R_{}}{\stackrel{~}{\omega }}\right)^\alpha \mathrm{exp}\left(\frac{1}{2}\left[\frac{z}{h(\stackrel{~}{\omega })}\right]^2\right),$$
(14)
where $`\stackrel{~}{\omega }`$ is the radial coordinate in the disk midplane, $`\rho _0`$ is the density at the stellar surface, and the scale height $`h`$ increases with radius:

$$h=h_0\left(\frac{\stackrel{~}{\omega }}{R_{}}\right)^\beta .$$
(15)
In order to reproduce the high density (and therefore temperature) gradient in the inner region of the disk, a very high resolution of the model both in radial and vertical direction has been applied. The smallest resolved structure at the inner radius of the disk has a linear extent of about 29 % both in vertical and radial direction. As for the SED modelling performed by Wood et al. (2002), the following values have been used: stellar radius $`R_{}=1.2R_{\mathrm{}}`$, stellar effective temperature $`T_{}=3500\mathrm{K}`$, $`\alpha =2.25`$, $`\beta =1.25`$, $`h_0=0.017R_{}`$, inner radius of the disk $`R_0=6R_{}`$, and the outer radius of the disk amounts to 200 AU. In contrast to the simulations presented by Cotera et al. (2001) and Wood et al. (2002), only one chemical dust grain component (astronomical silicate as in Sect. 4.1) has been chosen in order to simplify the simulation analysis. The size distribution is specified using a power law with exponential decay:

$$n(a)a^p\mathrm{exp}\left(\frac{a}{a_\mathrm{c}}\right)$$
(16)
(see, e.g., Kim et al. 1994). Following Cotera et al. (2001), a maximum grain size $`a_{\mathrm{max}}=20\mu \mathrm{m}`$ has been chosen (other parameter values: $`a_\mathrm{c}=0.55\mu \mathrm{m}`$, $`p=3.0`$). A total number of 32 grain sizes with radii equidistantly distributed on a logarithmic scale between 0.05 $`\mu \mathrm{m}`$ and 20 $`\mu \mathrm{m}`$ are considered. The visual extinction of the disk amounts to $`A_\mathrm{V}=2800`$ at an inclination of $`84^\mathrm{o}`$ (see Wood et al. 2002, Tabl. 2). The RT simulations have been performed on the basis of the model geometry shown in Fig. 6\[B\] in Wolf, Henning, & Stecklum 1999. While the star could be assumed to be point-like in the one-dimensional model (§4.1), its real extend has to be taken into account now. The radiation characteristic at each point of the stellar surface is described by a $`I\mathrm{cos}(\theta )`$ law (star without limb-darkening; see Wolf 2002 for details).

The resulting temperature distribution for two selected grain sizes ($`a=0.05\mu \mathrm{m}`$ and $`a=20\mu \mathrm{m}`$) and their difference are shown in Fig. 3. On the largest, $``$10–100 AU scale the temperature of the small grains is up to 40 K higher than that of the large grains in the optically thin atmosphere above/below the disk. In agreement with the results found in optically thick spherical shells, the temperature difference decreases towards the optically extremly thick midplane where it amounts to only a few Kelvin.

However, the vertical temperature structure shows a more complex behaviour in the inner region of the disk, on a size scale of a few AU and smaller. Here, an inversion of the temperature difference in the upper layer of the optically thick disk occurs, i. e., the large grains are substantially warmer than the small grains. This remarkable effect can be explained by the more efficient heating of large grains by the mid-infrared to far-infrared reemission of the hot inner disk. This explanation is also supported by the comparison of the vertical temperature structure obtained at different distances from the star333In order to obtain the spatially high-resolved temperature distributions shown in Fig.4 and 5 on a PC with 1 GByte RAM, the RT has been simulated only in the inner disk region with a diameter of 8.8 AU, which is one to two orders of magnitude larger than the regions considered in the following discussions. The comparison with the results obtained on a temperature grid with a lower resolution (Fig. 3) showed that the resulting differences in the spatial temperature distribution is in the order of or smaller than the statistical noise of the results and therefore negligible.: Fig. 4 shows the difference between vertical temperature distribution of grains with increasing radii ($`a=0.01\mu \mathrm{m}`$, 0.1 $`\mu \mathrm{m}`$, 1.1 $`\mu \mathrm{m}`$, 9.0 $`\mu \mathrm{m}`$) and $`a=20\mu \mathrm{m}`$ grains. Furthermore, Fig. 5 shows the temperature difference $`\mathrm{\Delta }T=T(a)T(a=20\mu \mathrm{m})`$ as a function of the grain radius $`a`$ and radial distance from the star at different distances from the midplane. While at a radial distance of 0.05 AU the temperature of the 20 $`\mu `$m grains is about 130 K higher than that of $`0.01\mu \mathrm{m}`$ grains in the “inversion layer”, the amount of this temperature difference drops to less than $`20\mathrm{K}`$ at a radial distance of 0.5 AU from the star (see Fig. 4). Consequently, the vertical extent and location of the minimum of the temperature inversion region depends on the particular grain size (and the vertical density distribution in the disk). For instance it was found that the minimum of the $`\mathrm{\Delta }T=T(a=9.0\mu \mathrm{m})T(a=20\mu \mathrm{m})`$ distribution shown in Fig. 4 is shifted towards smaller distances from the midplane in all three considered cases ($`\stackrel{~}{\omega }`$=0.05 AU, 0.1 AU, and 0.5 AU). According to the explanation of the temperature inversion given above, a similar shift of the minima of the difference temperature distribution is expected in case of the other grain sizes (0.01 $`\mu \mathrm{m}`$ to 1.1 $`\mu \mathrm{m}`$) as well. However, the spatial distribution of the simulated temperature distribution is too low to allow the verification of this assumption.

## 5 Spectral Energy Distribution

Based on the dust grain temperature distribution shown in Fig. 1, the spectral energy distribution (SED) for the different (spherical) shells have been calculated assuming $`n_\mathrm{S}=2^{i_\mathrm{S}}`$ ($`i_\mathrm{S}=1,2,\mathrm{},6`$) different grain sizes. The relative difference

$$\mathrm{\Delta }S_\lambda =\frac{S_\lambda (n_\mathrm{S})\overline{S}_\lambda }{\overline{S}_\lambda }$$
(17)
is shown in Fig. 6 and 7. Here, $`S_\lambda (n_\mathrm{S})`$ represents the SED of the real grain size distribution, and $`\overline{S}_\lambda `$ is the SED resulting from the approximation of mean dust grain parameters.

As Fig. 6 illustrates, the results of the real grain size distribution converges towards the results based on the approximation of mean dust grain parameters. As it also shows, a minimum of $`2^5`$ to $`2^6`$ grain sizes have to be considered to achieve deviations of less then $`10\%`$. If too few single grain sizes are considered, the SED does not represent the observable SED of the real grain size distribution sufficiently but drastically overestimates

1. The absorption of the stellar radiation (and therefore underestimates the visual to near-infrared SED),
2. The depths of the absorption bands in general (as the silicate absorption feature at $`9.7\mu `$m demonstrates), and
3. The dust reemission spectrum in the near-infrared to millimeter wavelength range.

While the different absorption behaviour influences the temperature structure and overall SED in general, the influence of the different scattering behaviour depends on the optical depth. In the optically thin case, differences in the scattering behaviour result mainly in deviations of the short wavelength region of the stellar SED. At higher optical depth, the scattering behaviour has direct influence on the energy transfer in the dust envelope and thus on the temperature structure and the near- to far-infrared wavelength range as well. However, this effect is of importance only at intermediate optical depths ($`\tau 1\mathrm{}10`$), since for even higher optical depths the temperature structure, at least in case of a one-dimensional dust configuration, tends not to depend on the description of the dust grain ensemble (real dust grain size distribution / approximation of mean dust grain parameters).

## 6 Conclusions

In simulations of the radiative transfer in the circumstellar environment of young stellar objects it has been widely established to use mean values for those parameters which describe the interaction of the electromagnetic field with the dust grains (see, e.g., modelling efforts by Calvet et al. 2002: circumstellar disk around the young low-mass star TW Hya; Cotera et al. 2001 and Wood et al. 2002: circumstellar disk around the classical T Tauri star HH 30; Fischer, Henning, & Yorke 1994, 1996: polarization maps of young stellar objects). A main reason for this is given by the fact that the consideration of a large number of single grain sizes and chemically different grain species results in a nearly linearly increasing amount of required computer memory in order to store the separate temperature distributions. Furthermore, the calculation of the spatial temperature distribution for all grains of different sizes and chemical composition requires significantly more computing power since (a) the heating by the primary sources (e.g., the star embedded in a circumstellar envelope) has to be performed independently for each grain species, and (b) the number of computing steps required to model the subsequent mutual heating of the different dust grain species due to dust re-emission scales even as $`(n_\mathrm{D}\times n_\mathrm{S})^2`$, where $`n_\mathrm{D}`$ is the number of chemically different components and $`n_\mathrm{S}`$ is the number of separate dust grain radii considered. While these simulations are feasible in case of one-dimensional models (see, e.g., Chini, Krügel, & Kreysa 1986 and Efstathiou & Rowan-Robinson 1994: Dust emission from star forming region; Krügel & Siebenmorgen 1994 and Siebenmorgen, Krügel, & Zota 1999: Radiative transfer in galactic nulei) or simple-structured two-dimensional models (Men’shchikov & Henning 1997: Circumstellar disks; Efstathiou & Rowan-Robinson 1994: Disks in active galactic nuclei), it is hardly possible to handle two- and three-dimensional models with high density gradients and/or high optical depth which require high-resolution temperature grids.

In this study, the difference between the results of RT simulations based (a) on a mean dust grain parameter approximation and (b) real dust grain size distributions have been investigated. Based on a one-dimensional density distribution it was found that the temperature structure of a real grain size distribution shows a very complex behaviour in the inner, hot region of the shell depending on (1) the grain size distribution, (2) the effective temperature of the embedded star and the optical depth and therefore on the density distribution. However, the relative difference between the SED based on a real dust grain size distribution on the one hand and the approximation of mean dust grain parameters on the other hand was found to be smaller than about $``$10% if a minimum number of $`2^5`$ to $`2^6`$ grain sizes have been considered.

As the temperature structure in a circumstellar disk – based on the model for HH 30 – shows, the geometry of the density distribution is a significant parameter for the resulting temperature differences between grains of different size, too. In the inner region of the disk with a diameter of a few AU a temperature inversion layer was found where the sign of the temperature difference of the largest and smallest grains is reverted. As this and the results obtained on the basis of the spherical shell (§4.1) show, the dust grain temperature structure is not sufficiently represented by dust grains with mean optical parameters. On the one hand, this is of tremendous importance for the simulation of chemical networks since the largest deviations from results based on the approximation of mean dust grain parameters have been found in the inner hot, dense region of the shell/disk where the chemical evolution takes place on its smallest timescale. On the other hand, the complex temperature structure may significantly change the hydrostatic properties of the considered gas/dust density distribution itself. However, these questions have to be investigated in future studies in order to find out the influence on observable quantities such as SEDs, images, polarization maps, and visibilities. Furthermore, the influence on processes taking place in more evolved circumstellar disks, such as the dust settling and dust grain growth, have to be considered taking into account the temperature structure of real grain size distributions.

This research was supported through the HST Grant GO 9160, and through the NASA grant NAG5-11645. I wish to thank the referee E. Dwek who helped to improve the clarity of the presentation of the results.


text:
# Leptonic Observables in Ultra-Relativistic Heavy Ion Collisions

## 1 PHENIX Experiment at RHIC

The PHENIX experiment was designed with the specific capability of sampling high luminosity proton-proton, proton-ion and ion-ion collisions. This feature, combined with a suite of detectors for measuring both electrons in two central arm spectrometers and muons in two forward arm spectrometers, allows for a large range of interesting measurements that probe the earliest stages of the hot quark-gluon matter created at RHIC. Further details of the detector design and performance are given in .

The PHENIX experiment consists of two central spectrometer arms that were completely instrumented in Run II at RHIC. The central arms each cover pseudorapidity ($`|\eta |<0.35`$), transverse momentum ($`p_T>0.2`$ GeV/c), and 90 degrees in azimuthal angle $`\varphi `$. They are comprised from the inner radius outward of a Multiplicity and Vertex Detector (MVD), Drift Chambers (DC), Pixel Pad Chambers (PC), Ring Imaging Cherenkov Counter (RICH), Time Expansion Chamber (TEC), Time-of-Flight Scintillator Wall (TOF), and two types of Electromagnetic Calorimeters (EMC). This combination of detectors is necessary for the clean identification of electrons over a broad range in transverse momentum. In addition, we have excellent photon capabilities with the EMC and high precision hadron identification with our smaller coverage TOF system.

PHENIX has two forward muon spectrometers consisting of hadronic absorbers, a cathode strip chamber muon tracking system, and interleaved Iarocci tubes with steel plates for muon identification and triggering. Each spectrometer covers approximately $`1.2<|\eta |<2.2`$ and $`p_{tot}>2`$ GeV/c. The south muon spectrometer was commissioned and operational in Run II, while the north arm is currently being completed and will be ready for data taking in Run III.

Au-Au and proton-proton collisions are characterized with minimum bias triggers based on a set of zero degree calorimeters (ZDC), beam-beam counters (BBC), and a scintillator multiplicity counter (NTC) for larger coverage in proton-proton collisions.

## 2 Run II at RHIC

During the second running period at RHIC the accelerator achieved full design energy $`\sqrt{s_{NN}}=200`$ GeV for Au-Au collisions. RHIC delivered to the PHENIX intersection region approximately 42 $`\mu \mathrm{b}^1`$ within a vertex $`z`$ range ($`|z|<45`$ cm). PHENIX was successful in sampling this luminosity with a combination of “minimum bias” triggers and a full complement of software based Level-2 triggers. As shown in Figure 1, over 50% of the integrated luminosity was delivered during the last two weeks of the Au-Au running period. Thus, we are seeing just the beginning of high luminosity running at RHIC. In the analysis shown in these proceedings, we have scanned of order 26 million “minimum bias” events with $`|z|<30`$ cm.

The RHIC accelerator also commissioned polarized proton-proton collisions and had five weeks of experiment data taking. The collisions were at $`\sqrt{s}=200`$ GeV for optimal comparison with the Au-Au data set. The PHENIX experiment commissioned a first set of fully pipelined hardware Level-1 triggers including triggers on high energy photons, single electrons, single muons, and dielectrons and dimuons from $`J/\psi `$. We observed with “minimum bias” and Level-1 triggers approximately 150 nb-1.

## 3 $`\varphi `$ Meson

The $`\varphi (1020)`$ vector meson is particularly interesting because the restoration of approximate chiral symmetry at high temperature may modify the $`\varphi `$ mass and width . These modifications may result in a change in the branching fraction of $`\varphi K^+K^{}`$ and $`\varphi e^+e^{}`$ when the $`\varphi `$ decays in medium. Note that the $`\varphi `$ lifetime $`\tau 44`$ fm/c is longer than the expected lifetime of the coupled collision system, and thus only a fraction may decay in the hot fireball. It has also been hypothesized that final state interactions of kaons from $`\varphi `$ decay may lower the apparent measured branching fraction in the kaon channel . PHENIX has the unique ability to reconstruct the $`\varphi `$ meson in both the $`K^+K^{}`$ and $`e^+e^{}`$ channels.

PHENIX has excellent kaon particle identification using our high precision time-of-flight scintillator wall. Shown in Figure 3 is the $`K^+K^{}`$ invariant mass distribution after mixed event background subtraction. The distribution is for minimum bias (0-90% central) Au-Au collisions at $`\sqrt{s_{NN}}=200`$ GeV. The signal to background ratio is 1/12 and the mass peak and width values agree within errors with the Particle Data Group values. More details are shown in . After correcting for acceptance, efficiencies, and the branching fraction in vacuum $`\varphi K^+K^{}`$ (0.49), we find:

$$\varphi (\mathrm{via}\mathrm{K}^+\mathrm{K}^{}):\frac{dN}{dy}|_{y=0}=2.01\pm 0.22(\mathrm{stat})_{0.52}^{+1.01}(\mathrm{sys})[\mathrm{PRELIMINARY}]$$
(1)
This value is in agreement with the STAR experiment results reported at this conference  within errors. We have higher statistics using time information from our larger coverage electromagnetic calorimeter for particle identification. We expect to have transverse mass distributions for a full range of Au-Au centrality bins shortly.

PHENIX has excellent electron identification capabilities that are necessary to separate electrons from the much more abundant charged pions. The RICH yields a threshold selection for electrons, the TEC allows for pion-electron separation at low momentum due to different energy loss $`dE/dx`$, and the EMC is used to verify a tracking momentum-electromagnetic energy match. We show in Figure 3 the $`e^+e^{}`$ invariant mass distribution after mixed event background subtraction. There is an excess of counts at the $`\varphi `$ mass with a signal strength of $`101\pm 47(\mathrm{stat})_{20}^{+56}(\mathrm{sys})`$ and a signal to background ratio of 1/20. Within relatively large errors the mass peak and width values agree with the values from the Particle Data Group. Using the signal value quoted above, we have corrected for the acceptance and efficiencies, and assuming the branching fraction in vacuum $`\varphi e^+e^{}`$ = $`2.9\times 10^4`$, we find:

$$\varphi (\mathrm{via}\mathrm{e}^+\mathrm{e}^{}):\frac{dN}{dy}|_{y=0}=5.4\pm 2.5(\mathrm{stat})_{2.8}^{+3.4}(\mathrm{sys})[\mathrm{PRELIMINARY}]$$
(2)
Although the electron and kaon results differ, they are consistent within $`1\sigma `$ statistical errors. We need further statistics in the electron channel from future running.

## 4 Dielectron Continuum

In the dielectron invariant mass region below the $`\varphi `$, referred to as the Low Mass Region (LMR), there may be excess dielectrons from in-medium mass modification of the $`\rho `$ meson due to restoration of approximate chiral symmetry or collision broadening . In the region above the $`\varphi `$ and below the $`J/\psi `$, referred to as the Intermediate Mass Region (IMR), there may be a significant contribution from semi-leptonic decays of open charm D mesons at RHIC energies.

We show in Figure 5 dN/dm into the PHENIX acceptance as a function of dielectron invariant mass after mixed event background subtraction with statistical and systematic errors. Note that we have corrected for our efficiencies relative to our PHENIX specific acceptance, and will publish our simple acceptance definition shortly. The solid band represents the systematic error in the normalization of the mixed event background, while the individual bands are the point to point systematics. We show the same distribution in Figure 5 with a simple “cocktail” model of expected baseline contributions (black), including light meson decays (red) and a PYTHIA charm calculation scaled by the number of binary collisions (blue).

In the LMR ($`M_{inv}=0.31.0`$ GeV/c2), the dominant expected contributions are from meson Dalitz decays and low mass vector meson decays (including the $`\omega `$). Assuming no medium modifications we expect a yield of $`N=9.2\times 10^5`$ (counts/event). The PHENIX preliminary result is:

$$N(\mathrm{LMR},\mathrm{into}\mathrm{PHENIX})=13.4\pm 7.2(\mathrm{stat})_{8.4}^{+12.2}(\mathrm{sys})\times 10^5[\mathrm{counts}/\mathrm{event}]$$
(3)
The large statistical and systematic errors prevent an observation of the $`\omega `$ meson or conclusions about $`\rho `$ mass shifts at this time.

In the IMR ($`M_{inv}=1.12.5`$ GeV/c2), the dominant expected contribution is from the correlated semi-leptonic decays of D and associated $`\overline{D}`$ mesons. Using a PYTHIA charm calculation, assuming binary collision scaling for ion-ion reactions, we would expect $`N=1.5\times 10^5(\mathrm{counts}/\mathrm{event})`$. Our preliminary result in the IMR is:

$$N(\mathrm{IMR},\mathrm{into}\mathrm{PHENIX})=0.38\pm 2.60(\mathrm{stat})_{0.81}^{+1.40}(\mathrm{sys})\times 10^5[\mathrm{counts}/\mathrm{event}]$$
(4)
This result is consistent with the charm expectation within errors, but clearly requires further statistics to demonstrate a signal. PHENIX has many handles on the production of heavy quarks, and the first results are via the measurement of single electrons.

## 5 Charm via Single Electrons

The measurement of open charm and open beauty in heavy ion collisions is sensitive to the initial gluon densities in the incoming nuclear wavefunctions and possible thermal contributions from later gluon-gluon scattering. In addition, total charm production is a crucial comparison by which to gauge the possible suppression of $`J/\psi `$ in deconfined matter. Total charm production (integrated over all transverse momentum) may be expected to be suppressed relative to binary scaling due to shadowing of the gluon parton distribution in nuclei or enhanced due to thermal production in later stages. Also, it was first predicted that high transverse momentum charm quarks will lose energy in a dense gluonic medium, and then later postulated that this energy loss is tempered by the “dead-cone” effect that suppresses co-linear gluon emission for $`v<<c`$ partons . Also, if the charm mesons appear in medium they may undergo hadronic interactions.

Electrons resulting from semi-leptonic decays of charm D mesons ($`De+K+\nu `$) and beauty B mesons ($`Be+D+\nu `$) have significant contributions to inclusive measured electrons above 1.0 and 2.5 GeV, respectively. The dominant sources of “background” electrons are from neutral pion Dalitz decay ($`\pi ^0\gamma e^+e^{}`$) and photon conversions in material in the PHENIX aperture. We can account for these sources using our measured $`\pi ^0`$ distributions. In addition, although their contributions are relatively small, we must model additional electrons from $`\eta ,\eta ^{},\omega ,\varphi ,\rho `$.

We have published  our results from Run I Au-Au collisions at $`\sqrt{s_{NN}}=130`$ GeV and find a significant excess of electrons over “background” above $`p_T0.7`$ GeV/c that increases with $`p_T`$. The lower data points in Figure 6 are the minimum bias transverse momentum distribution of non-photonic electrons. The results are consistent within relatively large errors with a PYTHIA 6.152  charm calculation using CTEQ5L parton distribution functions and assuming that charm scales with the number of binary collisions.

The production mechanism within PYTHIA is dominantly the fusion of gluons to produce a $`c\overline{c}`$ pair. The charm quark then fragments into a D meson which can decay semi-leptonically yielding an electron. However, an alternative picture is that the charm quarks re-scatter in the medium and take part in hydrodynamic expansion. We have used PHENIX derived hydrodynamic model parameters from our measured pion, kaon and proton data, and calculated the shape of the D meson distribution. After simulating the D meson decays, we find equally good agreement with our Run I single electron transverse momentum shape.

In Run II, we have a more powerful method of extracting the single electron charm contribution. Most of the “background” sources are photonic in origin. Thus, in Run II we had a special set of runs with a brass photon converter wrapped around the beam pipe near the interaction point. This material increases by a fixed factor the number of electrons whose source is photonic. Since there is a fixed relation between the $`\gamma `$ from $`\pi ^0`$ decay (which is the dominant source of conversion photons), and the Dalitz contribution $`\pi ^0\gamma e^+e^{}`$, we can subtract out both contributions. This method has many advantages over the previously described “cocktail” subtraction method, and has completely independent systematic errors. Also, this method subtracts out any direct photon contribution to conversion electrons, while the “cocktail” method does not.

After subtraction of the photonic contributions, we show as the upper points in Figure 6 the minimum bias transverse momentum distribution of non-photonic electrons. Again, the data appear reasonably described by a PYTHIA calculation of the expected charm contribution assuming binary collision scaling. The electron data in four centrality bins appear to qualitatively follow binary collision scaling as described in . Significant work remains to reduce the systematic errors in order to exactly quantify this scaling relation. PHENIX is currently analyzing the proton-proton data set so that we can calculate the charm scaling without relying on the PYTHIA baseline.

At this point, we can make the following qualitative observations. (1) Our electron data are consistent with charm expectations and binary scaling within our current statistical and systematic errors. (2) The NA50 experiment at the CERN-SPS has observed an enhancement in intermediate mass dimuons . One possible interpretation is a factor of three enhancement of charm in central Pb-Pb collisions relative to proton-proton. We do not see such an effect at RHIC energies. (3) PHENIX reports a substantial suppression in high transverse momentum $`\pi ^0`$ relative to binary scaling . We do not see such a large suppression in the single electrons from charm. This could be due to a lower charm quark energy loss in the dense gluonic medium due to the “dead-cone” effect. Further quantification of these observations is forthcoming.

## 6 J/$`\psi `$ Physics

We expect a screening of the QCD attractive potential as we approach the deconfinement transition. The color screening may result in a decrease in the production of heavy quarkonia ($`J/\psi ,\psi ^{},\mathrm{{\rm Y}}`$, etc.) . Alternatively, there are models that predict an enhancement of heavy quarkonia due to $`c\overline{c}`$ coalescence as the collision volume cools . This enhancement has a large dependence on the total charm production.

### 6.1 Proton-Proton Results

We present here the first measurements of $`J/\psi `$ production in proton-proton collisions at $`\sqrt{s}=200`$ GeV. We show in Figure 9 the dielectron invariant mass distribution for proton-proton collisions with at least one 2x2 calorimeter tower tile above our Level-1 trigger threshold. This represents about half of our proton-proton statistics, some of which are from additional 4x4 calorimeter tower tile Level-1 triggers. We show in Figure 9 the dimuon invariant mass distribution for proton-proton collisions.

We have corrected for the central and muon arm acceptances, tracking and particle identification efficiencies, and Level-1 trigger efficiencies. We show in Figure 10 the $`J/\psi `$ branching fraction to leptons (B) times $`d\sigma /dy`$ as a function of rapidity. The central rapidity point is from our dielectron measurement and the muon data is divided into two rapidity bins. The data are fit to a simple Gaussian form and also to the predicted rapidity shape from PYTHIA using parton distribution function GRV94LO. Removing the branching fraction and integrating over rapidity, we find:

$$\sigma (ppJ/\psi +X)=3.8\pm 0.6(\mathrm{stat})\pm 1.3(\mathrm{sys})\mu \mathrm{b}[\mathrm{PRELIMINARY}]$$
(5)
### 6.2 Gold-Gold Results

Shown in Figure 9 is the dielectron invariant mass distribution for “minimum bias” (0-90% central) Au-Au collisions, along with a mixed event distribution. Further details are in . Although the statistics are limited, we have calculated $`B_{ee}dN/dy|_{y=0}`$ in three Au-Au centrality ranges, where $`B_{ee}`$ is the $`J/\psi `$ branching fraction to electrons. The ranges correspond to 0-20%, 20-40%, and 40-90% central with calculated average number of binary collisions of 791, 297, and 45, respectively. We show the $`B_{ee}dN/dy|_{y=0}`$ per binary collision as a function of the collision centrality (represented by the number of calculated participating nucleons) in Figure 11. Note that we have also included our midrapidity $`J/\psi `$ measurement in proton-proton collisions with one binary collision. The statistical errors are shown as lines and the systematic errors as brackets.

On the same figure, we show different possible models of the $`J/\psi `$ centrality dependence. All the models have been arbitrarily normalized to perfectly intersect the center of our proton-proton data point. This graphical presentation can be somewhat misleading since it effectively assumes no error on the proton-proton measurement. However, we can perform an unbiased measure of agreement using a statistical confidence level. The statistical confidence level is defined as the probability, based on a set of measurements, that the actual probability of an event is better than some specified level. We find that the statistical confidence level is 16%, 80%, and 75% for (1) a model assuming binary collision scaling of the $`J/\psi `$, (2) a model assuming binary collision scaling followed by “normal” nuclear absorption with a 7.1 mb cross section, and (3) the exact pattern of suppression relative to binary scaling observed by the NA50 experiment at lower energies , respectively. It is notable that although the binary collision scaling model is statistically less favorable, all three models have acceptable confidence levels without consideration of our current systematic errors. Also, we hope to measure the “normal” nuclear absorption in deuteron-nucleus reactions at RHIC in the next two years.

There are two extreme models of $`J/\psi `$ production at RHIC. One is a simple extrapolation from models used to describe the lower energy NA50 results that imply that almost all $`c\overline{c}`$ correlations are broken up and one should have over an order of magnitude more suppression . Another is that although all initial $`J/\psi `$ are destroyed, they are re-created via coalescence of uncorrelated $`c\overline{c}`$ in the late stages of the time evolution. This mechanism can lead to an enhancement of $`J/\psi `$ relative to binary collisions (which is statistically less likely given our data).

With the data set from Run II, we can expect a factor of two more statistics from our Level-2 dielectron triggered data in Au-Au, in addition to results at forward rapidity from our dimuon measurements. Also, we are now confident in our abilities to trigger at high RHIC luminosity, and eagerly await a larger data set.

## 7 Summary

In summary, the PHENIX program of leptonic observables has started with data from Run II. We have the first observation of $`J/\psi `$ at $`\sqrt{s_{NN}}=200`$ GeV in both proton-proton and Au-Au systems. We have shown results on $`\varphi `$ production, dielectron continuum, and single electrons from charm D meson decay. We have shown just the tip of the iceberg of the PHENIX leptonic observables, and are excited about prospects for future high luminosity running at RHIC.


text:
# Degenerated Third Order Linear Recurrences

## 1 Introduction

Let us consider a third order linear recurrence

$$U_n=a_1U_{n1}+a_2U_{n2}+a_3U_{n3},U_0=u_0,U_1=u_1,U_2=u_2.$$

Let us assume that the roots of the characteristic polynomial are real, simple, but two of them are equal in absolute value: we call this situation degenerated. If $`\{\lambda _i\}`$ are the roots let us assume that $`\lambda _3>0`$, $`\lambda _1=\lambda _3`$ and $`\lambda _1<\lambda _2<\lambda _3`$. Then the coefficients of the recurrence are $`a_1=\lambda _2,a_2=\lambda _3^2,a_3=\lambda _2\lambda _3^2`$.

The closed form (generalized Binet’s formula) is

$$U_n=c_1(\lambda _3)^n+c_2\lambda _2^n+c_3\lambda _3^n,$$

where $`\{c_i\}`$ are determined by the initial condition $`\{u_i\}`$, as well as by the roots $`\{\lambda _i\}`$.

## 2 Limits

Using Binet’s formula the ratio of two consecutive terms is given by

$$\frac{U_n}{U_{n1}}=\frac{c_1(\lambda _3)^n+c_2\lambda _2^n+c_3\lambda _3^n}{c_1(\lambda _3)^{n1}+c_2\lambda _2^{n1}+c_3\lambda _3^{n1}}.$$

After dividing both numerator and denominator by $`\lambda _3^{n1}`$ we get

$$\frac{U_n}{U_{n1}}=\frac{c_1\lambda _3(1)^n+c_2\lambda _2\left(\frac{\lambda _2}{\lambda _3}\right)^{n1}+c_3\lambda _3}{c_1(1)^{n1}+c_2\left(\frac{\lambda _2}{\lambda _3}\right)^{n1}+c_3}.$$

Now the limit for $`n\mathrm{}`$ depends on the parity of $`n`$. If $`n`$ is even

$$L_1=\underset{n\mathrm{}}{lim}\frac{U_n}{U_{n1}}=\lambda _3\frac{c_1+c_3}{c_1+c_3},$$
(1)
while if $`n`$ is odd

$$L_2=\underset{n\mathrm{}}{lim}\frac{U_n}{U_{n1}}=\lambda _3\frac{c_1+c_3}{c_1+c_3}.$$
(2)
Let us express $`c_1`$ and $`c_3`$ explicitly as functions of the initial conditions. Using Binet’s formula we have

$$\{\begin{array}{ccc}c_1+c_2+c_3\hfill & =\hfill & u_0,\hfill \\ c_1\lambda _3+c_2\lambda _2+c_3\lambda _3\hfill & =\hfill & u_1,\hfill \\ c_1\lambda _3^2+c_2\lambda _2^2+c_3\lambda _3^2\hfill & =\hfill & u_2.\hfill \end{array}$$

Let A be the following matrix

$$𝐀=\left[\begin{array}{ccc}1& 1& 1\\ \lambda _3& \lambda _2& \lambda _3\\ \lambda _3^2& \lambda _2^2& \lambda _3^2\end{array}\right].$$

Let $`𝐜`$ the column vector with elements $`\{c_i\}`$ and $`𝐮`$ the column vector with elements $`\{u_i\}`$. Then

$$𝐜=𝐀^1𝐮.$$

We obtain

$$𝐀^1=\left[\begin{array}{ccc}\frac{\lambda _2}{2(\lambda _2+\lambda _3)}& \frac{1}{2\lambda _3}& \frac{1}{2\lambda _3^2+2\lambda _2\lambda _3}\\ \frac{\lambda _3^2}{\lambda _3^2\lambda _2^2}& 0& \frac{1}{\lambda _3^2+\lambda _2^2}\\ \frac{\lambda _2}{2\lambda _3+2\lambda _2}& \frac{1}{2\lambda _3}& \frac{1}{2\lambda _3^22\lambda _2\lambda _3}\end{array}\right].$$

Then

$$c_1+c_3=\frac{\lambda _2^2u_0+u_2}{\lambda _3^2\lambda _2^2},$$

and

$$c_3c_1=\frac{\lambda _3^2(\lambda _2u_0+u_1)+\lambda _2(\lambda _2u_1+u_2)}{\lambda _3^3\lambda _3\lambda _2^2}.$$

It follows

$$L_1=\frac{\lambda _3^2(\lambda _2^2u_0u_2)}{\lambda _3^2(\lambda _2u_0u_1)+\lambda _2(\lambda _2u_1u_2)},$$

$$L_2=\frac{\lambda _3^2(\lambda _2u_0u_1)+\lambda _2(\lambda _2u_1u_2)}{\lambda _2^2u_0u_2}.$$

## 3 Consequences

If we write

$$\gamma =\frac{c_1+c_3}{c_1+c_3},$$
(3)
then

$$\frac{L_1}{\lambda _3}=\gamma ,\frac{L_2}{\lambda _3}=\frac{1}{\gamma }.$$

It follows

$`\lambda _3^2`$ $`=`$ $`L_1L_2`$ (4)
$`=`$ $`\underset{n\mathrm{}}{lim}{\displaystyle \frac{U_n}{U_{n1}}}\underset{n\mathrm{}}{lim}{\displaystyle \frac{U_{n1}}{U_{n2}}}`$
$`=`$ $`\underset{n\mathrm{}}{lim}{\displaystyle \frac{U_n}{U_{n1}}}{\displaystyle \frac{U_{n1}}{U_{n2}}}`$
$`=`$ $`\underset{n\mathrm{}}{lim}{\displaystyle \frac{U_n}{U_{n2}}}.`$
That is, this limit exists and it does not depend on the initial conditions.

Also

$`\gamma ^2`$ $`=`$ $`{\displaystyle \frac{L_1}{L_2}}`$ (5)
$`=`$ $`{\displaystyle \frac{\underset{n\mathrm{}}{lim}\frac{U_n}{U_{n1}}}{lim_n\mathrm{}\frac{U_{n1}}{U_{n2}}}}`$
$`=`$ $`\underset{n\mathrm{}}{lim}{\displaystyle \frac{U_n}{U_{n1}}}{\displaystyle \frac{U_{n2}}{U_{n1}}}`$
$`=`$ $`\underset{n\mathrm{}}{lim}{\displaystyle \frac{U_nU_{n2}}{U_{n1}^2}}.`$
That is, this limit exists and depends on the initial conditions.

## 4 Existence of a Limit

Under particular initial conditions our ratio converges. The equality $`L_1=L_2`$ entails

$$\lambda _3^2(\lambda _2^2u_0u_2)^2=[\lambda _3^2(\lambda _2u_0u_1)+\lambda _2(\lambda _2u_1u_2)]^2.$$

We can obtain $`u_2`$ as a function of $`u_0`$ and $`u_1`$. This leads to a quadratic equation in $`u_2`$ whose solutions are

$$u_2=\{\begin{array}{c}\lambda _2\lambda _3u_0+(\lambda _2\lambda _3)u_1\\ \lambda _2\lambda _3u_0+(\lambda _2+\lambda _3)u_1.\end{array}$$

If we take the first one then

$$L_1=L_2=\lambda _3.$$

On the other hand if we take the second one then

$$L_1=L_2=\lambda _3.$$


text:
# Saddle points in the chaotic analytic function and Ginibre characteristic polynomial

## References


text:
# Charge Particle Multiplicity and Transverse Energy Measurements in Au-Au collisions in PHENIX at RHIC

## 1 INTRODUCTION

$`N_{ch}`$ and $`E_T`$ are global variables which give excellent characterization of the high energy nucleus-nucleus collisions, thus providing information about the initial conditions . They help to constrain the wide range of theoretical predictions and discriminate among various mechanisms of particle production. $`N_{ch}`$ measurements also provide one with an opportunity to study high density QCD effects in relativistic nuclear collisions .

First results for $`N_{ch}`$ and $`E_T`$ at mid-rapidity in Au-Au collisions at $`\sqrt{s_{_{NN}}}=130`$ GeV measured with the PHENIX detector were published in . The same experimental techniques were used to analyze data obtained at $`\sqrt{s_{_{NN}}}=200`$ GeV. In order to make a more precise comparison of the results at the two beam energies, both data samples are processed through the same analysis procedure using a more restrictive event selection criteria. The trigger required the coincidence of two beam-beam counters and two zero degree calorimeters. The minimum bias trigger efficiency was found to be the same within 0.7% at $`\sqrt{s_{_{NN}}}=130`$ GeV and 200 GeV, equal to $`91.4_{3.0}^{+2.5}`$%. The corrections for the measured $`E_T`$ and $`N_{ch}`$ at both $`\sqrt{s_{_{NN}}}=130`$ GeV and 200 GeV related to particle composition and mean transverse momentum are performed based on PHENIX results obtained at $`\sqrt{s_{_{NN}}}=130`$ GeV , instead of HIJING  used in .

The part of the PbSc electromagnetic calorimeter used for the $`E_T`$ measurements covers the pseudorapidity range $`|\eta |0.38`$ with an azimuthal aperture of $`\mathrm{\Delta }\varphi =44.4^{}`$ in $`\sqrt{s_{_{NN}}}=130`$ GeV data and $`\mathrm{\Delta }\varphi =112^{}`$ in $`\sqrt{s_{_{NN}}}=200`$ GeV data. The pad chambers used for $`N_{ch}`$ measurements have a fiducial aperture of $`|\eta |0.35`$ and $`\mathrm{\Delta }\varphi =90^{}`$ and $`180^{}`$ in the two data sets respectively.

The results at $`\sqrt{s_{_{NN}}}=130`$ GeV presented in this paper and published earlier in  are consistent within systematic errors (excluding the 4% scale shift in $`E_T`$ due to its new definition222 In our previous publications , $`E_T`$ was defined using the kinetic energy for nucleons and the total energy for all other particles. This approach does not take into account the mass of nucleons resulting from pair production. In this paper, $`E_T`$ was defined with $`EM_N`$ for baryons, $`E+M_N`$ for antibaryons, and $`E`$ for all other particles, where $`E`$ is the total energy of the particle and $`M_N`$ is the nucleon mass. The revised definition of $`E_T`$ increases its value by about 4%, independent of centrality.).

## 2 RESULTS

Fig. 1 shows the centrality dependence of $`dE_T/d\eta `$ and $`dN_{ch}/d\eta `$ per participant pair measured at $`\sqrt{s_{_{NN}}}=130`$ GeV and 200 GeV. Both values show a steady rise with $`N_p`$.

The ratio of the pseudorapidity densities measured at $`\sqrt{s_{_{NN}}}=130`$ GeV and 200 GeV ($`R_{200/130}`$) for each centrality bin, corresponding to 5% of the nuclear interaction cross section, is shown in Fig. 2. For the most central bin the transverse energy increases by $`14\pm 2`$% and the charged particle multiplicity increases by $`17\pm 3`$%. The centrality dependence of the ratios is consistent with a constant.

Fig. 3 shows the comparison of our results to different model predictions. A powerful test for theoretical models is the comparison to the measured ratios $`R_{200/130}`$, since many systematic errors in experimental measurements cancel out. The increase of $`dN_{ch}/d\eta `$ with centrality is in contrast to the predictions of the EKRT model . HIJING  is in qualitative agreement with this scenario, however the strong centrality dependence of the ratio $`R_{200/130}`$ predicted by HIJING is excluded by the data. Our experimental results are well described by high energy QCD gluon saturation  and two-component mini-jet model  calculations.

Rapidity densities per pair of participants obtained in different experiments are shown in Fig. 4. The data are measured in the center-of-mass system for the RHIC experiments (BRAHMS, PHENIX, PHOBOS and STAR) or the laboratory system for all other experiments. It is assumed that $`dX/dydX/d\eta `$ in the laboratory system ($`X`$ stands for $`E_T`$ or $`N_{ch}`$), and a factor of $`1.20`$ derived from the HIJING generator was applied to account for the $`\eta y`$ transformation in the center-of-mass system. All RHIC results are in good agreement. The data points for $`E_T`$ and $`N_{ch}`$ are consistent with a logarithmic rise with $`\sqrt{s_{_{NN}}}`$ over a broad range of collision energies.

$`E_T`$ and $`N_{ch}`$ behave in a very similar manner such that the mean $`E_T`$ per charged particle remains unchanged over a broad range of centralities (see Fig. 5, left). The same behaviour of $`<E_T>/<N_{ch}>`$ with centrality was observed by the WA98 Collaboration at $`\sqrt{s_{_{NN}}}=17.2`$ GeV . The ratio stays almost unchanged also as a function of $`\sqrt{s_{_{NN}}}`$ (see Fig. 5, right).