text:
# Numerical simulation of the surface flow of a companion star in a close binary system

## 1 Introduction

A semi-detached binary is a system consisting of a primary star and a Roche-lobe filling secondary star. In cataclysmic variables (CVs), for instance, the primary star is a white dwarf, while the Roche-lobe filling star is a low-mass main-sequence star, from which gas is supplied to the primary star via an accretion disc (see Warner 1995 for a review).

Accretion discs have been the main target of astrophysical research, while gas supplying companion stars has not been investigated fully from both theoretical and observational viewpoints. Knowledge of the properties of the companion star is required in order to understand the evolutionary mechanism of the binary systems. One reason that the companions have not attracted much attention from astrophysicists is that they had little observable constraints. Recently, however, Dhillon and Watson (dhi00 (2000)) have demonstrated the technique of imaging the companion stars in CVs using Roche tomography.

The companion stars in CVs are similar to low-mass main-sequence stars in their gross properties. However, companion stars exist in a different environment from those of isolated stars. They might be influenced by effects like irradiation from the accretion disc around the primary star, its rapid rotation, Roche-lobe shape, mass loss, and so on. These effects could in turn affect the evolution of the binary systems.

### 1.1 Astrostrophic wind of Lubow and Shu

Lubow and Shu (lub75 (1975), lub76 (1976)) conducted pioneering work in the study of gas dynamics in semi-detached binaries. In their paper, they analyzed the surface layer of the companion star, as well as the L1 region (Lubow and Shu lub75 (1975)) and the L1 stream (Lubow and Shu lub76 (1976)). In their analysis they solved the hydrodynamic equations by a matched expansion technique and reduced partial differential equations to ordinary differential equations. Their method was superior to numerical methods at that time considering the poor capability of computers. Since then, the progress in computer power has been considerable, so that it is worthwhile to investigate the same problem numerically. Our main target in the present paper is to investigate the surface layer, the L1 region and the L1 stream numerically. The region around the accreting compact star has been investigated separately by our group (Makita, Miyawaki and Matsuda mak00 (2000), Matsuda et al. mat00 (2000), Fujiwara et al. fuj01 (2001)).

Lubow and Shu (lub75 (1975)) stated that “the horizontal component of the fluid velocity is parallel to isobars in the astrostrophic approximation. Like the large-scale circulation patterns in the atmosphere and oceans of the Earth, the flow does not simply proceed from high pressures to low, but is directed around the contours of equal pressure. In the northern hemi-lobe, the flow proceeds by keeping the region of high (low) pressure to the right (left); in the southern hemi-lobe, by keeping high (low) pressure to the left (right).” On a given equi-potential surface they expected “the pressure to be largest at the poles (for the same reason that meridional circulation occurs in rotating stars) and smallest at L1 (because of the mass-loss flow)”; therefore “the excess circulation should be counter to the sense of the orbit rotation.”

They predicted that “on each equipotential surface, the flow may be astrostrophic, i.e., parallel to isobars, but a net transfer of matter within the entire surface layer to the equator is possible … if there is an outward flux of matter at the bottom of the surface layer.” They also considered that “another effect which reinforces the concept that there must be a slow drift of matter from the poles to the equator, and thence to L1, occurs if the envelope of the contact component is convective. In this case, the effective frictional drag exerted on the upper layers by the lower layers must produce the “Ekman effect” (Ekman ekm22 (1922)) where the induced secondary flow corresponds to a drift across isobars from high pressures to low.”

One of the aims of the present paper is to investigate the flow pattern of the gas in the surface layer of the companion numerically and to compare the results with the qualitative prediction of Lubow and Shu (lub75 (1975)). This kind of problems may be termed stellar meteorology. The structure of the L1 region and the L1 stream are also considered.

### 1.2 Our previous work

Two-dimensional finite-volume hydrodynamic simulations of flows in a semi-detached binary were done by Sawada, Matsuda and Hachisu (saw86 (1986), saw87 (1987)). Three-dimensional calculation was done by Sawada and Matsuda (saw92 (1992)). They were concerned mainly by the accretion disc and did not pay much attention to the surface flow on the secondary companion.

Recently, Fujiwara et al. (fuj01 (2001)) performed three-dimensional hydrodynamic simulations of a semi-detached binary system including the companion star. The equation of state is that of an ideal gas characterized by the specific heat ratio $`\gamma `$ and the case with $`\gamma =1.01`$, meaning that the gas is almost isothermal, was mainly studied. The reason for the choice of $`\gamma =1.01`$ is to mimic a radiative cooling effect occurring in an accretion disc. Although their main concern was the hydrodynamics of the accretion disc, they also obtained the flow pattern on the surface of the companion. Their results showed that gas flows from the pole region to the equatorial region rather directly, although small eddy patterns appeared on the surface of the companion star. One, however, has to be careful to select the values of the specific heat ratio. In the surface region of the companion, where the gas is expanding outward, the temperature of the gas is reduced due to the adiabatic expansion. Thus the case with $`\gamma =1.01`$, which results in the input of energy into the expanding gas, may not be appropriate for the present purpose.

This paper is organized as follows. In $`\mathrm{\S }`$ 2 we describe the numerical method and physical assumptions. In $`\mathrm{\S }`$ 3 we show the results of our numerical simulations. In $`\mathrm{\S }`$ 4 we discuss the dynamics of the obtained eddy system. A summary is given in $`\mathrm{\S }`$ 5.

## 2 Assumptions and Numerical Method

### 2.1 Assumptions

We consider the present study as a first step in the investigation of the surface flow of the companion, and therefore do not take into account such complex effects as viscosity (except numerical), magnetic fields, and irradiation from the accretion disc. The equation of state considered here is that of an ideal gas characterized by a specific heat ratio $`\gamma `$, and the cases with $`\gamma =5/3`$ and 7/5 are mainly studied, although the cases of $`\gamma =1.2`$ and 1.01 are also considered to compare the present results with those of Fujiwara et al. (fuj01 (2001)).

In the present study, the binary system consists of a primary star with mass $`M_1`$ and a companion star with mass $`M_2`$ rotating counter-clockwise when viewed from the north pole. The mass ratio of the two stars, $`q=M_2/M_1`$, is assumed to be unity in the present work. We normalize all physical quantities as follows: the length is scaled by the separation $`A`$ between the centres of the two stars. The time is scaled by $`1/\mathrm{\Omega }`$, where $`\mathrm{\Omega }`$ denotes the orbital angular velocity, and therefore the orbital period becomes $`2\pi `$. The density at the inner boundary is taken to be unity. The gravitational constant $`G`$ is eliminated using the above normalization.

### 2.2 Numerical method

We use a Cartesian coordinate system. The origin of the coordinate is located at the centre of the companion star. The $`x`$axis coincides with the line joining the centre of the two stars. We set the $`xy`$ plane as the orbital plane, thus the $`z`$axis is perpendicular to the orbital plane and is oriented in the same direction as the angular momentum vector of the orbital rotation. The L1 point is located at ($`0.5,0,0`$) and the primary star is located at ($`1,0,0`$), which is out of the present computational region.

The computational region is a rectangular box of size $`0.5<x<0.7`$, $`0.5<y<0.5`$ and $`0<z<0.5`$. Note however that part of this box contains the mass-losing star, for which we assume that the inner part contains uniform density gas at rest (see Sect. 2.3). With the assumption of symmetry of physical quantities around the orbital plane, the calculations are performed only in the upper half region above this plane. This region is divided in $`121\times 101\times 51`$ grid points. In order to test the assumption of symmetry about the orbital plane, we performed a test in which the whole space was taken into account in the computation. We confirmed that the assumption of the symmetry was good enough.

In the present study, we use the same scheme as in Makita, Miyawaki, and Matsuda (mak00 (2000)), Matsuda et al. (mat00 (2000)) and Fujiwara et al. (fuj01 (2001)): the simplified flux splitting (SFS) finite-volume method proposed by Jyounouchi et al. (jou93 (1993)) and Shima and Jyounouchi (shi94 (1994)). We refer to Makita, Miyawaki, and Matsuda (mak00 (2000)) for the numerical and the SFS schemes. With this MUSCL-type technique, we can keep the spatial and temporal accuracy at second-order levels.

### 2.3 Boundary and initial conditions

We apply the same boundary conditions as in Fujiwara et al. (fuj01 (2001)) except for the inner boundary location of the companion star. The inner boundary is assumed to be an equi-potential surface slightly smaller than the critical Roche lobe: the mean radius of the inner boundary is chosen to be 0.45 in our dimensionless unit. Note that this inner boundary does not necessarily mean the real surface of the companion, but is chosen only for numerical purposes. We have also made simulations where the inner boundary was chosen to be 0.40 and 0.50 and the results are very similar. In the rest of this paper, we will thus concentrate on the case with an inner boundary of 0.45.

What we call here the “surface” of the companion star might not correspond to the photosphere, where optical depth reaches unity. Lubow and Shu (lub75 (1975)) pointed out that “the depth over which the optical depth at visual wavelengths reaches unity is not a priori related to the depth over which appreciable mass flow occurs.” In the present study, we are interested in the surface flow of the companion, and therefore we define the “surface” as the region where appreciable flow occurs.

The inside of the companion star is filled with gas having zero velocity, density $`\rho _0=1`$, and sound speed $`c_o`$, where $`c_o`$ is a free parameter. We assume $`c_0=0.05`$ in the present work, although a case with $`c_0=0.1`$ is also calculated and does not show much difference in the results. Such a value for the sound speed might, at first glance, look very large, if one would only take into account the effective temperature of the star and scale the corresponding speed with the orbital velocity. This is however not correct because as we noted above, the inner boundary is not the ”surface” of the star but may well lie more inside the star, where the gas is much hotter. Temperatures of several 105 K corresponding to sound speeds of the order of 50 km/s or more are therefore not unrealistic.

The assumption that the gas inside the companion star has zero velocity does not mean that there is no gas outflow from the companion star. The gas does flow out when the pressure above the inner boundary is lower than that inside the companion. The outflow flux is calculated by solving Riemann problems at each step.

The outside of the outer boundary is assumed to be always filled with a gas with velocity $`0`$, density $`\rho _1`$, and sound speed $`c_1`$, where $`\rho _1`$ and $`c_1`$ are parameters. Use of these boundary conditions ensures a stable calculation. Note that inflow or outflow of the gas through the outer boundary is also possible, and the flux is calculated by solving Riemann problems between the inside and outside of the computational region at each step.

At the initial time $`t=0`$, the entire region, except the inside of the companion, is occupied by a gas with velocity $`0`$, density $`\rho _1=10^5`$, and sound speed $`c_1=\sqrt{10}`$.

To know how gas elements from the companion build the L1 stream is interesting, because this stream does bring the material from the companion to the primary via the L1 point and an accretion disc. Thus we investigate the structure of the L1 stream as well as the surface flow in the present work.

Simulations are run up to $`t=62.8`$ ($`10`$ orbital periods) except for the finer grid case, in which the simulation is run up to $`t=12.56`$. We find that the main flow patterns reach a steady state by these times.

## 3 Numerical results

### 3.1 Astrostrophic wind

Figures 1(a)-(d) show the equi-density surface of the companion star and the streamlines starting from the equi-density surface; Fig. 1(a) depicts the equi-density surface of $`\mathrm{log}_{10}\rho =1`$, (b) $`1.5`$, (c) $`2`$, and (d) $`2.5`$. In order to draw the streamlines, we use the same integration time step. Therefore, the higher the density of gas, the shorter the streamlines as the velocity is lower. We observe that there are ascending motions of gas in the higher density region (deeper inside the star) and a remarkable circulating flow in the lower density region (surface region). The streamlines of lighter color have higher velocity than those of darker ones. As will be shown later, there is a high pressure around the north pole, and so we may call this circulatory flow the north pole circulatory flow or simply an H-eddy.

In order to see the surface flow more closely, Fig. 2 shows the equi-density surface ($`\mathrm{log}_{10}\rho =2`$) of the companion star and the streamlines starting from the equi-density surface of $`\mathrm{log}_{10}\rho =2.5`$ seen from various view angles. We confirm the presence of the H-eddy again. Gas in the H-eddy seems to drift gradually downward and its velocity increases with decreasing latitude. This phenomenon will be discussed later.

Figure 3 shows the iso-velocity lines plotted on the equi-density surface of $`\mathrm{log}_{10}\rho =2.5`$ seen from the north pole. The numbers in the figure show the magnitude of the velocity (note that the velocity is normalized by $`A\mathrm{\Omega }`$, where $`A`$ is the separation and $`\mathrm{\Omega }`$ is the rotational frequency of the binary.)

In the L1 region, which is close to the L1 point, there is a strong low pressure due to the mass outflow through the L1 point. Because of the presence of the L1 low pressure, there is a circulatory flow, around the L1 point, rotating counter-clockwise, although this eddy does not form a complete circle for geometrical reasons. We call this flow the L1-eddy, although, strictly speaking, this is not a complete eddy.

The other remarkable feature is the presence of two eddies rotating counter-clockwise at the opposite side of the L1 point. These two eddies merge to form one larger eddy at different density levels. Since this region is close to the L2 point, we may call this eddy the L2 eddy. Since these H, L1 and L2 eddies are explained in terms of the astrostrophic wind (Lubow and Shu lub75 (1975)), the pressure distribution on the surface of the companion is important.

### 3.2 High / low pressures on the companion surface

If there is no flow on the companion, the isobaric surface should coincide with the equi-potential surface as a result of the hydrostatic balance. In our case, the isobaric surface does not coincide with the equi-potential surface in an exact sense because of the presence of the surface flow. Fig. 4 shows the isobaric lines plotted on the equi-potential surface of the companion. This equi-potential surface nearly corresponds to the equi-density surface with $`\mathrm{log}_{10}\rho =2.5`$ except near the L1 region, where high speed gas motion is present. Note that the isobaric lines on the companion surface have concentric patterns around the north pole.

As was described earlier, the flow pattern on the surface of the rotating (in the inertial frame) companion can be explained in terms of the astrostrophic wind, which is a counterpart of the geostrophic wind in the atmosphere of the Earth. If viscosity can be neglected, the direction of the flow is parallel to the isobaric lines of low/high pressure. In the present case of a counter-clockwise rotating binary system (as viewed from the north pole), the flow proceeds by keeping the high (low) pressure to the right (left) in the northern hemisphere. Therefore, the H-eddy circulates clockwise, while the L1 eddy and L2 eddy circulate in the counter-clockwise direction.

If the pressure is only a function of the density of gas, the gas is called barotropic. Otherwise, the gas is called baroclinic. In the case of $`\gamma =1.01`$, the gas is nearly iso-thermal, hence almost barotropic. In the other cases, the gas is baroclinic. In a barotropic gas, an isobaric surface must coincide with an equi-density surface. Fig. 5 shows the baroclinicity of the gas with $`\gamma =5/3`$, whose property is responsible for the generation of vorticity. This will be discussed later.

### 3.3 Ekman boundary layer

In the Earth’s atmosphere, an Ekman boundary layer exists close to the surface of the Earth, where (turbulent) viscosity cannot be neglected and where the direction of the flow deviates from that of a pure geostrophic wind. In low pressure regions, gas drifts inwards and eventually ascends near the centre of the low pressure regions (Pedlosky ped87 (1987)).

In our simulations, there is no (physical) viscosity assumed. Nevertheless, numerical viscosity arising from a truncation error acts as a kind of a turbulent viscosity, leading to gas drifting. This is the reason why the circulating flow about the north pole gradually migrates towards the equatorial region. In order to test the above consideration, we perform a calculation with the mesh size halved. In this case the magnitude of the (numerical) viscosity is expected to be $`1/8`$ of the coarser grid case. We find a similar flow pattern as before, although the circulating motion is tighter than before (see Figs. 6). The main characteristics of the flow pattern are not changed, though.

In the surface region of low-mass main-sequence stars, there is a convection zone, which provides the necessary turbulent viscosity. We may therefore expect that the present flow pattern is qualitatively correct.

It is however not trivial from our simulations to derive the exact value of our (numerical) viscosity. We can only make a very tentative guess that the typical Reynolds number in our simulations is probably of the order of 10 to 100. In a future work, we plan to solve this problem using the Navier-Stokes equation instead of the Euler ones, and using the molecular viscosity as a rough estimate of the turbulent viscosity which should be at play here.

In our case the inner boundary does not play the role of the surface of the Earth. The Ekman layer is present at the top of the atmosphere, a situation similar to a water vortex. Gas circulating around the high pressure region at the north pole migrates gradually towards the equator, switches to the counter-clockwise motion close to the L1 low pressure, and finally escapes from the companion towards the primary star. This flow pattern is important to understand the structure of the L1 stream to be discussed later.

### 3.4 Effect of $`\gamma `$

Until now we have shown the results for $`\gamma =5/3`$, corresponding to ionized hydrogen gas. In order to study the case of molecular hydrogen, we calculate the case $`\gamma =7/5`$. The results are essentially the same as for $`\gamma =5/3`$, so that we do not discuss them any further.

However, the results for $`\gamma =1.01`$ are very different from the cases of $`\gamma =5/3`$ and 7/5, and are essentially the same as that of Fujiwara et al. (fuj01 (2001)): there is no H eddy but the flow originating at the pole region migrates rather directly towards the equatorial region. More precisely, such a flow is rather weak, because the pressure on an equi-potential surface is rather uniform except in the region close to the L1 point. The L1 flow mainly originates from the region close to the L1 point. The case of $`\gamma =1.2`$ is intermediate to the above two extreme cases. We found the L2 eddy in this case.

We may explain the difference of the flow patterns between the cases of $`\gamma =5/3,7/5`$ and that of $`\gamma =1.01`$ in terms of the baroclinicity of the gas. Vorticity is generated by the term containing $`\rho \times p`$, which is identically zero for a barotropic gas (Pedlosky ped87 (1987)). Therefore, local vorticity of the gas in the case of $`\gamma =1.01`$ is almost zero everywhere. Moreover, the Taylor-Proudman theorem, which states that the horizontal velocity components $`u`$ and $`v`$ are independent of $`z`$, holds not only for an incompressible flow but also for a barotropic gas (Busse bus70 (1970); Pedlosky ped87 (1987)). Since $`u`$ and $`v`$ are almost zero on the inner boundary, they must also be small in the upper atmosphere. This problem, however, needs further investigation. We note that the cases of $`\gamma =5/3`$ and 7/5 are more appropriate to compare with (future) observations.

### 3.5 L1 region and L1 stream

Lubow and Shu (lub75 (1975)) paid lots of effort to solve the flow in the L1 region by their matched expansion technique. In this subsection we summarize our numerical results.

In order to make the influence of the outer boundary to the L1 stream small, we slightly expand the computational region of the $`x`$-component from $`0.7`$ to $`0.9`$ in this subsection. The primary star is still out of the computational region. Thus the L1 stream reaches the boundary at $`x=0.9`$ and escapes from the calculation region. In the present calculation no accretion disc is formed, and so there is no collision between the L1 stream and the accretion disc around the primary star.

Lubow and Shu (lub75 (1975)) predicted that the sonic point, where subsonic gas is accelerated to supersonic velocity, must occur close to the L1 point. Figure 7 shows the iso-Mach number lines on an equi-density surface and the numbers on the figure show the Mach number. It can be seen that the sonic line on the surface is certainly close to the L1 point.

Figs. 8(a)-(c) show the velocity vectors, the streamlines, and the equi-density lines of the gas in the orbital plane, respectively. The flow pattern in the orbital plane near the L1 point shows some difference with the one depicted in Fig. 3 in Lubow and Shu (lub75 (1975)). They predicted that the streamlines of the lower portion of the L1 point turn back to the companion, while in our case all gas flow towards the primary compact component. There are two differences between thier assumptions and ours: 1) they assumed iso-thermal gas; 2) they start their calculation at the L1 point, where the gas is assumed to be at rest. The assumption 2) may be the main reason. We solve the whole region including both the L1 region and the surface layer of the companion. The gas at the L1 point is already supersonic as was described above.

From Fig. 8(b) we can measure the deflection angle $`\theta _s`$ of the L1 stream, which is the angle between the central flow line in the L1 stream and the $`x`$axis. It is about -10 degrees. This angle was calculated to be about -20 degrees for the case of $`q=1`$ by Lubow and Shu (lub75 (1975)). In Fujiwara et al.’s result (fuj01 (2001)) this is about -15 degrees. This difference is also due to the difference in boundary conditions. Lubow and Shu assumed that the fluid element was at rest at the L1 point. The effect of the Coriolis force is measured by the Rossby number defined by $`U/A\mathrm{\Omega }`$, where $`U`$ is a typical velocity of gas. The lower the initial velocity of the gas, the lower the Rossby number, the larger the Coriolis’ effect and, hence, the deflection angle.

In Fig. 8(c) we found that the atmosphere of the secondary is thicker in the leading hemisphere (negative $`y`$) than the trailing hemisphere. This can be easily explained by the direction of the Coriolis force, which acts to the right side of the flow direction in our system.

Because of this asymmetry of the atmosphere, the structure of the L1 stream is also asymmetric about the central flow line. In Figs. 9 we show the cross section of the L1 stream, which is bow shaped convex to the positive $`y`$ direction. We obtained a similar pattern in Fujiwara et al. (fuj01 (2001)), and this was explained by the interaction of the L1 stream and the circulating circum-disc flow on the orbital plane about the primary. In the present work we realize that this explanation is not correct. The bow shape of the L1 stream should be explained by the structure of the surface flow on the companion.

Lubow and Shu (lub75 (1975)) predicted a symmetric Gaussian density distribution at the centre of the stream as

$$\rho (r,\theta )\frac{C}{r}\mathrm{exp}[Br^2(\theta \theta _s)^2],$$
(1)
where $`r`$ is the distance from the L1 point in the scaled coordinate, $`\theta _s`$ is the angle of the stream centre from the line joining the centre of the two stars, and $`C`$, $`B`$ are the coefficients depending on the mass ratio. Fig. 10 shows that our density profile nearly has a Gaussian density distribution. The density at the centre of the stream obtained by them decreases as $`r^1`$, while our result decreases as $`r^{0.8}`$. Fig. 11 shows the equi-density surface of $`\mathrm{log}\rho =2.5`$ of the companion star and the L1 stream.

## 4 Discussion

In our simulations, we obtained a system of eddies: the H eddy rotating clockwise, and counterclockwise rotating L1 and L2 eddies. We now discuss what we believe is the reason that this eddy system is formed (see Fig. 12).

If we consider a time sequence, the L1 low pressure is formed first, because of the mass loss through the L1 point. Since the gas near the equator does not feel much resistance in flowing towards the L1 point, the equatorial region becomes a low pressure zone, and a high pressure zone builds up at the north pole: the H eddy is thus formed.

The above explanation does not explain the appearance of the third eddy: the L2 eddy. This can be explained as follows. Consider a force exerted by an eddy on another eddy. The L1 eddy tries to convect the H eddy to negative $`y`$-direction, while the H eddy tries to convect the L1 eddy to negative $`y`$-direction as well. However, the L1 eddy must be fixed to the L1 region, so that the H eddy must drift to negative $`y`$-direction. The system is therefore not steady. In order to stabilize the eddy system, we need the third eddy, which tries to convect the H eddy to positive $`y`$-direction. The L2 eddy thus appears and evolves until the three eddy system is stabilized. The exact location and the strength of these eddies can be obtained only through numerical simulations.

It is interesting to note that Davey and Smith (dav92 (1992)) have found that the trailing hemisphere has stronger line absorption than the leading hemisphere in a dwarf nova (see also Dhillon and Watson dhi00 (2000)). This asymmetry can be naturally explained in terms of the H-eddy current, which convects heat due to irradiation to the leading hemisphere.

If this kind of astrostrophic wind conveys heat effectively from an irradiated hemisphere to the other, the unirradiated side of the companion becomes hotter than previously thought. Such an example may be observed in the 2000 outburst of the recurrent nova CI Aql because, in nova outbursts, a companion star should be strongly irradiated by a hot white dwarf. The inclination angle of the CI Aql orbit is about 80 degrees and the eclipse is almost total in the sense that the companion completely occults the bright accretion disc around the white dwarf (e.g., Hachisu & Kato 2001). In eclipse minima, therefore, we see only the unirradiated side of the companion. Hachisu, Kato and Schaefer (2002) showed that the eclipse minima of CI Aql in outburst is about 0.6 mag brighter than in quiescence and could reasonably reproduce the orbital light curve (folded by the orbital period) when the temperature of the unirradiated hemisphere is about 1000 K hotter than that in quiescence.

In the present study, we do not take into account such complex effects as (real) viscosity, magnetic fields, and irradiation from the accretion disc. Activities of magnetic fields on the companion surface would lead to star-spots, and the irradiation from the accretion disc would lead to an increase of the gas temperature of the companion. These effects might more or less change the surface flow patterns. The present study is a first step to examine the surface flow structure. Thus to study these phenomena is beyond the scope of the present paper, and is left for future work.

## 5 Summary

We performed numerical simulations of the surface flow on the companion star in a semi-detached binary system. The case with $`\gamma =5/3`$ was mainly studied, although the cases of $`\gamma =7/5,1.2`$ and 1.01 were also investigated. The results are summarized as follows:

1. We obtain an eddy configuration composed of an H-eddy, an L1-eddy, and an L2-eddy associated with high/low pressures. If there is no viscosity, the flow lines must be closed. The above eddies are explained in terms of the astrostrophic wind.
2. Some kind of viscosity (e.g., turbulent viscosity originating from a surface convection zone) would lead to a downward drifting motion. In our simulation the numerical viscosity acts as turbulent viscosity.
3. The sonic line on an equi-density surface occurs close to the L1 point.
4. The deflection angle of the L1 stream is about -10 to -15 degrees. This small angle is due to the larger fluid velocity at the L1 point.
5. The density distribution in the L1 stream is nearly symmetric about the central flow line. The cross section of the L1 stream is bow shaped convex to the positive $`y`$ direction. This is an intrinsic nature of the stream and is not caused by the interaction with the accretion disc gas.

###### Acknowledgements.

The authors would like to thank Prof. Takahiro Iwayama for his valuable discussions on meteorological aspects of the present study. The referee, Dr. S. Lubow, is gratefully acknowledged for valuable comments. T.M. was supported by the grant in aid for scientific research of the Japan Society of Promotion of Science (13640241). H.B. would like to thank the Kanbara-Tosao foundation for its financial support. Calculations were carried out on SGI Origin 3800 at the Information Processing Center of Kobe University.


format: text

text:
# Phenomenology of Jet Quenching in Heavy Ion Collisions

## Abstract

We derive an analytical expression for the quenching factor in the strong quenching limit where the $`p_T`$ spectrum of hard partons is dominated by surface emission. We explore the phenomenological consequences of different scaling laws for the energy loss and calculate the additional suppression of the away-side jet.

It is commonly believed that the yield for “hard” observables in high-energy nuclear reactions scales as the number of binary nucleon-nucleon (NN) collisions occurring during the encounter of the two nuclei. This expectation applies to processes that are characterized by a high virtuality $`q^2`$, for which final state interactions are negligible, such as lepton pair production or the total yield of heavy flavor quarks. There is no a priory reason to expect this rule to hold for high-$`q^2`$ processes, in which the final state can be strongly modified by interactions with comovers. A well documented example is the production of heavy vector mesons, such as the $`J/\mathrm{\Psi }`$, which is found to be “anomalously” suppressed in collisions of heavy nuclei at the CERN-SPS GH99 .

The yield of hadrons produced with high transverse momentum $`p_T`$ in Au+Au collisions at the Relativistic Heavy Ion Collider (RHIC) has recently been shown to be significantly suppressed in comparison with the cumulative yield of NN collisions PHEN02 ; STAR02 . This effect, called “jet quenching”, was predicted to occur as a result of energy loss by the hard scattered partons due to interactions with the surrounding dense medium Bjo82 ; Tho91 ; Gyu94 . The theory of this energy loss has been a topic of intense research over the past few years Zak96 ; BDMPS ; Wie00 ; BDMS ; BSZ00 . The present consensus is that the dominant mechanism for the energy loss in QCD is collisionally induced radiation of gluons by the fast parton.

It is difficult to measure the energy loss of a scattered parton directly in heavy ion reactions, because the large multiplicity of emitted hadrons makes it almost impossible to isolate the resulting jet by kinematic cuts. However, the energy loss of the parton is imprinted as an equivalent loss of energy of the leading hadron produced in its fragmentation Wan92 . This is what has been observed in the RHIC experiments. Generally, it is assumed that the fragmentation occurs after the parton has left the comoving medium, and thus is described by the measured vacuum fragmentation functions. We will, therefore, not be concerned with the conversion from partons to hadrons, but focus directly on the $`p_T`$-spectrum of scattered partons.

Preliminary data from Run 2 at RHIC (Au+Au at a center-of-mass energy of 200 GeV per nucleon pair) confirm the effect observed in Run 1 and its interpretation as jet quenching QM02 . For pions with $`p_T5`$ GeV/c the measured suppression factor $`Q(p_T)`$ is about 1/5. There is not yet complete agreement between different experiments about the scaling of the hadron yields with the number of nucleons participating in the reaction, $`N_{\mathrm{part}}`$, or with the binary collision number $`N_{\mathrm{coll}}`$. The PHOBOS collaboration has presented evidence that the yield of charged hadrons with $`p_T4`$ GeV/c scales like $`N_{\mathrm{part}}`$ PHOB-QM . This finding is contradicted by data from the PHENIX collaboration PHEN-QM1 . It is not clear whether the discrepancy is due to different normalization methods, different experimental acceptance, or other reasons.

In addition, the experiments have found that the supression factors for mesons and baryons are quite different up to transverse momenta of 5 Gev/c. It was recently proposed that this phenomenon can be attributed to a competition between different hadron formation processes Fries03 ; Greco , with parton recombination dominating at lower $`p_T`$ and fragmentation at higher $`p_T`$. The participant or binary collision scaling may be complicated in the transition region. The situation is predicted to simplify for $`p_T>6`$ GeV/c, where hadrons are overwhelmingly produced by fragmentation of an energetic parton.

We will first show analytically that the spectrum of high-$`p_T`$ partons, in the limit of large energy loss, is dominated by partons emitted from the surface of the collision zone and thus scales like the surface rather than the volume of the interaction region Shu01 . We will then explore the scaling of the quenching factor $`Q(p_T)`$ with participant number and $`p_T`$. We finally calculate the additional suppression of the away-side jet and the azimuthal anisotropy of the parton yield in noncentral collisions.

We begin by considering the loss of energy by an energetic parton traversing a homogeneous, static medium of thickness $`L`$. We assume that the geometry is given by a cylinder with radius $`R`$, as in the boost-invariant Bjorken model Bjo83  for a nuclear collision with impact parameter $`b=0`$, and that the parton moves in the transverse plane in the local rest frame of the medium. We further assume that the effective energy loss, defined as the shift of the momentum spectrum of fast partons, depends on $`p_T`$ and $`L`$ in the following general way:

$$\mathrm{\Delta }p_T=\eta p_T^\mu L,$$
(1)
where $`\mu `$ is a scaling exponent. The linear dependence on $`L`$ holds when the loss occurs in subsequent, independent interactions with the medium. It has also been shown to be valid, when multiple interactions are suppressed by the LPM effect, due to the steep fall-off of the parton spectrum with $`p_T`$ BDMS . The traversed pathlength inside the medium for a parton created at transverse position $`𝐫`$ with an angle $`\varphi `$ relative to the radial direction is

$$L(\varphi )=(R^2r^2\mathrm{sin}^2\varphi )^{1/2}r\mathrm{cos}\varphi \frac{z^2}{2R\mathrm{cos}\varphi },$$
(2)
where $`r=|𝐫|`$, $`z^2=R^2r^2`$, and the approximation is valid near the surface ($`zR`$). Perturbative QCD predicts that the parton spectrum at moderatedly large values of $`p_T`$ has the form FMS02

$$\frac{dN}{d^2p_T}=N_0\left(1+\frac{p_T}{p_0}\right)^\nu $$
(3)
with a power $`\nu 8`$ and $`p_01.75`$ GeV/c. The quenched spectrum is given by

$`{\displaystyle \frac{d\stackrel{~}{N}}{d^2p_T}}`$ $`=`$ $`Q(p_T){\displaystyle \frac{dN}{d^2p_T}}`$ (4)
$`=`$ $`{\displaystyle \frac{1}{2\pi ^2R^2}}{\displaystyle _0^{2\pi }}𝑑\varphi {\displaystyle _0^R}d^2r{\displaystyle \frac{dN(p_T+\mathrm{\Delta }p_T)}{d^2p_T}}.`$
Replacing the integration over $`r`$ with an integration over $`z`$, using the approximation (2), and formally extending the range of the integration to infinity, one finds:

$$Q(p_T)\frac{2(p_0+p_T)}{\pi R\eta (\nu 1)p_T^\mu }.$$
(5)
Two things are remarkable about this result. First, the factor $`R`$ in the denominator reduces the scaling of the parton yield with the size of the reaction zone by one power of $`R`$, from a volume to a surface dependence. Second, the dependence of $`Q(p_T)`$ on $`p_T`$ is determined by the power $`\mu `$ governing the $`p_T`$-dependence of the energy loss. For $`\mu =1`$ the quenching factor $`Q`$ falls slowly with increasing $`p_T`$; for $`\mu =1/2`$ it grows with $`p_T`$, implying less quenching at higher $`p_T`$.

We have confirmed the range of validity of the approximate analytical expression (5) by comparing it with the exact integral (4). As seen in Fig. 1, the analytical approximation $`Q_{\mathrm{anal}}`$ deviates from the exact result by less than 5% when the quenching factor $`Q0.2`$. An even better agreement is found when the analytical result is divided by the correction factor $`(1+Q_{\mathrm{anal}}^2)`$. The excellent agreement suggests that we may extend the calculation to noncentral collisions. Generalizing the surface-to-volume ratio of the cylindrical geometry ($`S/V=2/R`$) to the geometry of a collision with impact parameter $`b`$, one finds

$$Q(p_T,b)=Q(p_T,0)\alpha _b/(\alpha _b\mathrm{sin}\alpha _b)$$
(6)
with $`\alpha _b=2\mathrm{arccos}(b/2R)`$. This approximation is valid with about the same accuracy as (5), except for very peripheral collisions.

In order to be able to address the experimental data, we need to relax some of the geometrical oversimplifications used above.

1. The transverse profile of the primary jet yield is not homogeneous, but proportional to the binary NN collision profile

$$T(𝐫,𝐛)=\rho _1(𝐫)\rho _2(𝐫𝐛)$$
(7)
where $`\rho _i(𝐫)`$ is the longitudinally integrated density of nucleus $`i`$.
2. The density of the comoving medium is also not homogeneous in the transverse plane. We assume that it is proportional to the local density of participant nucleons which, in the Glauber model, is given by

$`\rho _{\mathrm{part}}(𝐫)`$ $`=`$ $`\rho _1(𝐫)\left(1e^{\sigma \rho _2(𝐫𝐛)}\right)`$ (8)
$`+\rho _2(𝐫𝐛)\left(1e^{\sigma \rho _1(𝐫)}\right),`$
where $`\sigma `$ denotes the inelastic NN cross section.
3. The comover medium expands and its density decreases with time. Here we assume

$$\rho (𝐫,\tau )=C\rho _{\mathrm{part}}(𝐫)/(\tau +\tau _0),$$
(9)
which is modeled to represent a longitudinal, boost invariant expansion. As an approximation to the phenomenology at full RHIC energy we use the values $`C3`$ and $`\tau _0=1`$ fm/c.

Finally, we need to return to equation (1) for the energy loss. Perturbative QCD predicts that the radiative energy loss depends quadratically on the medium thickness $`L`$ BDMPS . As pointed out by Baier et al. BDMS , this holds for the average energy loss $`\overline{\mathrm{\Delta }E}`$ of a given parton, but the average energy loss of observed partons with fixed transverse momentum $`p_T`$ has a different scaling. This is so, because the energy loss distribution $`D(ϵ)`$ is strongly skewed toward small values of $`ϵ`$ by the steeply falling $`p_T`$-spectrum of fast partons. In fact, the average shift of the spectrum due to the energy loss, here called the effective energy loss, is given by BDMS

$$\mathrm{\Delta }p_T\alpha _sL\sqrt{\pi \widehat{q}p_T/\nu },$$
(10)
where $`\widehat{q}`$ encodes the “scattering power” of the medium, which is proportional to the density. For an expanding medium, the expression $`\widehat{q}L^2`$ must be replaced with

$$\widehat{q}_0L_{\mathrm{eff}}^2=2\widehat{q}_0_0^L\tau 𝑑\tau \frac{\rho (𝐫(\tau ),\tau )}{\rho (𝐫,0)},$$
(11)
where $`𝐫(\tau )=𝐫+𝐯\tau `$ denotes the position of the fast parton in the medium at time $`\tau `$, and $`\widehat{q}_0`$ is a function of the transverse position $`𝐫`$ at which the jet is produced. We thus can make contact with (1) by writing

$$\mathrm{\Delta }p_T=\eta ^{}L_{\mathrm{eff}}\sqrt{\rho (𝐫,0)p_T/\nu }$$
(12)
with the constant $`\eta ^{}=\alpha _s\sqrt{\pi \widehat{q}_0/\rho (0)}`$, which does not depend on $`𝐫`$.

We will denote the scaling law (12) for the energy loss as BDMS. In our following numerical study we have explored two other scaling laws. The first one is the Bethe-Heitler (BH) scaling law Jeon02 ; Mus98

$$\mathrm{\Delta }p_T=\eta p_T_0^L𝑑\tau \rho (𝐫(\tau ),\tau )\eta p_T(L\rho )_{\mathrm{eff}}$$
(13)
corresponding to $`\mu =1`$. The second scaling law is

$$\mathrm{\Delta }p_T=\eta p_T\sqrt{(L\rho )_{\mathrm{eff}}},$$
(14)
which we will call the (RW) scaling law. It could be interpreted as describing a random walk in $`p_T`$ as the fast parton traverses the medium, with some interactions resulting in an energy gain and others in a loss of energy.

We begin the discussion of our numerical results for the quenching factor $`Q`$ with its dependence on the transverse momentum of the fast parton, shown in Fig. 2. The QCD-motivated BDMS law (solid line) and the other two scaling laws exhibit clearly different behaviors. This reflects the different $`p_T`$ scaling of the energy loss in these models (linear for BH and RW; square-root for BDMS). The data from the RHIC experiments PHEN02 ; STAR02 ; PHEN-QM1  suggest that the quenching first becomes stronger with increasing momentum, reaches a minimum, and finally begins to diminish. This would indicate that the BDMS law only applies at high $`p_T`$, and that other laws govern the energy loss at lower $`p_T`$ Jeon02  or hadron production is not dominated by parton fragmentation in this kinematic region. We note that the dependence of $`Q`$ on $`p_T`$ is well described by the analytical formula (5).

The impact parameter dependence of the quenching factor is shown in Fig. 3, plotted as the yield per half the number of particpant nucleons against the particpant number $`N_{\mathrm{part}}`$. The unquenched jet yield, scaling with the number of binary NN collisions, would increase relative to $`N_{\mathrm{part}}`$. As the figure shows, the quenching counteracts this increase, and the yield per participant actually falls for the BDMS and the BH laws as the collision centrality increases. An approximately flat behavior, as observed in the PHOBOS experiment PHOB-QM , is only found for the RW scaling law.

For noncentral collisions, the quenching factor $`Q`$ is a function of the azimuthal emission angle, because the geometry is not axially symmetric. This is known to lead to an angular asymmetry of a quadrupole shape in the spectra of high-$`p_T`$ particles Wan01 ; Gyu01 . The elliptic flow parameter $`v_2`$ is defined as the Fourier component proportional to $`\mathrm{cos}(2\varphi )`$ of the angular distribution of particles with respect to the scattering plane Pos98 . We find (see Fig. 4) that the values of $`v_20.1`$ obtained for all three scaling laws are significantly smaller than the measured values ($`v_20.2`$) for semicentral and peripheral collisions STAR02a . However, the calculated $`v_2`$ for partons would be large enough to explain the measured elliptic flow of hadrons, if the hadrons were produced by recombination in the $`p_T`$-range of the RHIC data Molnar .

As observed hadrons from hard partons preferentially originate from the surface region facing the detector, the parton emitted in the opposite direction has to traverse more material and thus endures an even larger energy loss. This leads to an additional suppression of the away-side jet and its leading hadron spectrum. The dependence of the incremental away-side hadron suppression factor $`Q_{\mathrm{asj}}`$ on the primary suppression factor $`Q`$ is shown in Fig. 5 for the BDMS scaling law for two impact parameters ($`b=0`$ and $`b=8`$ fm). There appears to be a universal relationship, which is linear for $`Q0.2`$.

Several conclusions can be drawn from our results. First, the momentum dependence of the BDMS energy loss formula does not seem to agree with some of the RHIC data. A linear dependence of the average energy loss on $`p_T`$ is in better agreement with the data in the $`p_T`$ region explored so far Jeon02 . Also, the dependence of the effective momentum loss $`\mathrm{\Delta }p_T`$ on the medium thickness $`L_{\mathrm{eff}}`$ predicted by BDMS does not yield a scaling of the charged hadron yield with participant number as seen in the PHOBOS data up to $`p_T=4.25`$ GeV/c; only the RW scaling law yields such a dependence. This may indicate that energy gain and loss mechanisms are competing in this momentum range WW01 . Another explanation may be that hadrons at intermediate momenta are not produced by fragmentation of fast partons, but by other processes, such as parton recombination. The observed magnitude of the elliptic flow lends support to this interpretation. Finally, we have found a universal relation between the same-side and away-side suppression factors for the BDMS law, which can be tested experimentally.

Acknowledgments: This work was supported in part by a grant from the U. S. Department of Energy (DE-FG02-96ER90945). I thank S. A. Bass, X.-N. Wang, and K. Rajagopal for valuable comments on the manuscript.


text:
# Migration and Accretion of Protoplanets in 2D and 3D Global Hydrodynamical Simulations

## 1. Introduction

Migration has entered the puzzling scenario of planetary formation as the favorite mechanism advocated to explain the extremely short orbital period of many extrasolar planets.

It is known that any planet-like body is forced to adjust its distance from the central star because of gravitational interactions with the circumstellar material. However, the dispute about how fast migration proceeds is far from being over. Numerical methods have been employed to evaluate gravitational torques exerted on embedded planets. We have performed a series of simulations modeling both two and three dimensional disks, varying the mass $`M_p`$ of the protoplanet in the range from 1 Earth-mass to 1 Jupiter-mass.

The physics of the problem demands that the flow in the protoplanet’s neighborhood should be accurately resolved. In order to achieve sufficient resolution, even for very small planetary masses, we use a Nested-Grid technique (D’Angelo, Henning, & Kley 2002a). This paper addresses the issues of flow circulation around protoplanets, orbital migration, and mass accretion.

## 2. Numerical Model

We assume the protostellar disk to be a viscous fluid (viscosity $`\nu =10^{15}\mathrm{cm}^2\mathrm{s}^1`$) and describe it through the Navier-Stokes equations (Kley, D’Angelo, & Henning 2001; D’Angelo et al. 2002a). The set of equations is integrated over a grid hierarchy, as shown in Figure 1. The planet is supposed to move on a circular orbit at $`r_p=5.2\mathrm{AU}`$ around a solar-mass star. The disk has an aspect ratio $`h=0.05`$ and the mass within the simulated region is $`M_D=3.5\times 10^3M_{}`$.

## 3. 2D Simulations of Uranus and the Earth

A circumplanetary disk forms inside the Roche lobe of massive as well as low-mass protoplanets. Such structures are characterized by a two-arm spiral wave perturbation. They are detached from the circumstellar disk spirals, which arise outside of the planet’s Roche lobe. Indicating with $`l`$ the distance from the planet normalized to $`r_p`$, for a wide range of planetary masses the spiral pattern can be approximated to

$$\mathrm{\Theta }\mathrm{\Theta }_0=2k\left(1/\sqrt{l_0}1/\sqrt{l}\right),$$
(1)
where $`k=\zeta \sqrt{M_p/M_{}}/h`$ and $`\zeta 1`$. The ratio $`k/\sqrt{l}`$ represents the Mach number of the circumplanetary flow. Figure 2 (left panel) demonstrates how equation (1) fits to the spiral perturbation around a Uranus-mass planet.

Even an Earth-mass planet induces two weak spirals, which wrap around the star for $`2\pi `$, but no circumplanetary disk is observed (Figure 2, right panel).

## 4. 3D Simulations of Saturn and Neptune

A more complete description of the flow near protoplanets is provided by 3D computations (D’Angelo, Kley, & Henning 2002b). The major differences between 2D and 3D modeling arise in the vicinity of the planet because the latter can account for the vertical circulation in the circumplanetary disk. Instead, the two geometries decently agree on length scales larger than the disk scale height.

Two examples of our simulations are illustrated in Figure 3. The images represent the logarithm of the density close to a Saturn-mass planet (left) and Neptune-mass planet (right), in two orthogonal planes (see Figure caption for details). The velocity field is overplotted to display the flow features. From the top panels it is clear that the spiral perturbations are weaker and more open than in 2D simulations. The bottom panels indicate the presence of vertical shock fronts, which are generally located outside the Hill sphere of the protoplanet.

## 5. Migration and Accretion

In general, torques exerted by disk material cause a protoplanet to migrate toward the star (Ward 1997). Yet, nearby matter can be very efficient at slowing down its inward motion (Tanaka, Takeuchi, & Ward 2002; D’Angelo et al. 2002a). The migration time scale can be defined as $`\tau _\mathrm{M}r_p/|\dot{r}_p|`$, where the migration drift $`\dot{r}_p`$ is directly proportional to the total torque acting on the planet. Since we also measure the rate $`\dot{M}_p`$ at which the planet accretes matter from its surroundings, an accretion time scale can be introduced as well: $`\tau _\mathrm{G}M_p/\dot{M}_p`$.

Some of our 2D and 3D outcomes for both time scales are shown in Figure 4.

## 6. Summary

Circumplanetary disk forms around protoplanets. The spiral wave pattern which marks such disks is less accentuated when the full 3D structure is simulated. Vertical shock fronts develop outside the Hill sphere of the planet.

The estimated values of $`\tau _\mathrm{M}`$ in 3D are longer than those predicted by analytical linear theories because of non-linearity effects. When $`M_p<30M_{}`$, both $`\tau _\mathrm{M}`$ and $`\tau _\mathrm{G}`$ are longer in 3D computations then they are in 2D ones.

## References

D’Angelo, G., Henning, Th., & Kley, W. 2002a, A&A, 385, 647

D’Angelo, G., Kley, W., & Henning, Th. 2002b, ApJ, submitted

Kley, W., D’Angelo, G., & Henning, Th. 2001, ApJ, 547, 457

Tanaka, H., Takeuchi, T., & Ward, W. 2002, ApJ, 565, 1257

Ward, W. 1997, Icarus, 126, 261


text:
# Radiation driven atmospheres of O-type stars: constraints on the mass–luminosity relation of central stars of planetary nebulae

## 1. Introduction

A lot of work on model atmospheres of PN central stars (CSPNs in what follows) has been motivated by the desire to obtain information about the basic properties of CSPNs (surface temperature, mass, luminosity, abundances), so as to be able to test predictions from post-AGB evolutionary calculations. The earlier efforts, based on plane-parallel non-LTE models, could not achieve a completely independent test, in the following sense: since the plane-parallel model fits to H and He photospheric absorption lines can only produce information about surface temperature, He abundance and $`\mathrm{log}g`$, we cannot derive stellar masses or luminosities, but only $`L/M`$ ratios. This is exactly the same problem we face when dealing with low-gravity early-type “supergiant” stars at high Galactic latitudes: are they luminous and massive, or are they evolving away from the AGB? We need some independent evidence to settle the issue – for example, the distance to the star. Unfortunately, we lack reliable distances to most CSPNs.

So what could be done was to plot the positions of CSPNs in the $`\mathrm{log}g`$$`\mathrm{log}T_{\mathrm{eff}}`$ diagram, and compare them with plots of post-AGB tracks, translated from the $`\mathrm{log}L`$$`\mathrm{log}T_{\mathrm{eff}}`$ diagram. After doing this translation it is possible to read the stellar mass in the $`\mathrm{log}g`$$`\mathrm{log}T_{\mathrm{eff}}`$ diagram. From this, we can derive $`L`$ and, if we know the visual dereddened apparent magnitude, a so-called “spectroscopic distance”. All this work, however, is based on assuming that the evolutionary models give us the correct relation between stellar mass and luminosity. It is not a real test of the evolutionary models, but only a consistency check.

In the last 10 years there has been a lot of progress in the modelling of stellar winds driven by radiation pressure. Many CSPNs show spectroscopic evidence of winds, in the form of emission lines in the visible spectrum and especially P-Cygni-type profiles of resonance lines in the ultraviolet between 1000 and 2000 Å.

The existence of these wind features provides both a challenge and an opportunity. The challenge is to model them. The opportunity is to use the information about the geometrical extension of the atmosphere and the forces of gravity and radiative pressure, implicit in the wind profiles (from which the terminal velocity and mass loss rate can be derived), to obtain the physical size of the star, which is the key to derive the stellar luminosity and mass. (The idea is described together with a first application by Pauldrach et al. 1988.) Thus the successful modelling of the wind features opens the way for a real test of the mass–luminosity relation of CSPNs.

In this review we would like to present the current situation of the project and the rather surprising results obtained up to now. Section 2 describes a model analysis of a massive Population I star, using state-of-the-art hydrodynamically consistent, spherically symmetric model atmospheres to demonstrate the power of the technique and to show how successfully we can reproduce the ionizing fluxes and observed spectra of such stars. In Section 3 we introduce the wind-momentum–luminosity relation and describe a previous attempt to determine whether or not CSPNs follow this relation. In Section 4 we show two examples of the application of our hydrodynamically consistent wind models to CSPNs (fits to IUE and HST spectra and derived parameters). In Section 5 we add the results from 6 other CSPNs similarly analyzed and discuss the consequences.

## 2. Modelling a massive O supergiant: $`\alpha `$ Cam

The analysis method is described in detail by Pauldrach et al. (2001).111 This paper is also available on the Web, at
http://www.usm.uni-muenchen.de/people/adi/adi.html. Here we can only give a brief summary. The analysis method is based on modelling a homogeneous, stationary, extended, outflowing, spherically symmetric radiation driven atmosphere. A complete model atmosphere calculation involves solving the hydrodynamics and the NLTE problem (rate equations, radiative transfer). The solution of the total interdependent system of equations is obtained iteratively. This permits the calculation of the predicted or synthetic spectrum, which is then compared to the observed UV spectrum. The process is repeated with different stellar parameters until a satisfactory fit is obtained.

The UV spectrum between 1000 and 2000 Å carries a lot of information: P-Cygni-type profiles of resonance lines of several ions of C, N, O, Si, S, P, as well as hundreds of strongly wind-contaminated lines of Fe iv, Fe v, Fe vi, Cr v, Ni iv, Ar v, Ar vi. But the information about the stellar parameters can be extracted only after careful analysis. A very important recent improvement of our method concerns the development of a substantially consistent treatment of the blocking and blanketing influence of all metal lines in the entire sub- and supersonically expanding atmosphere. All the results we will present are based on this new generation of models.

The operational procedure is as follows. A preliminary inspection of the visual or UV spectrum of the star to be analyzed gives an initial estimate of $`T_{\mathrm{eff}}`$. From the UV spectrum, the terminal wind velocity $`v_{\mathrm{}}`$ can be measured directly. Now an initial value for the stellar radius $`R`$, defined at a Rosseland optical depth of $`2/3`$, is assumed. Using this $`R`$ and $`v_{\mathrm{}}`$ we obtain an estimate of the mass ($`v_{\mathrm{}}`$ scales with $`\sqrt{M/R}`$ as explained by the theory of radiation driven winds). With the current values of $`R`$, $`T_{\mathrm{eff}}`$, $`M`$, and assuming a set of abundances, we can solve the model atmosphere and calculate the velocity field, the mass loss rate $`\dot{M}`$, and the synthetic spectrum. Now this predicted spectrum is compared to the observed one. If the fit is not satisfactory, we need to modify $`\dot{M}`$ via a change of $`R`$ (since $`\mathrm{log}\dot{M}\mathrm{log}L`$, according to the radiation driven wind theory). The change in $`R`$ forces us to change the mass, too, in order to keep $`v_{\mathrm{}}`$ consistent with the observed value. The new model is calculated and the process is repeated until we obtain a good fit to all features in the observed spectrum.

Figure 1 shows the result of applying this procedure to the O supergiant HD 30614 ($`\alpha `$ Cam). The final parameters are $`T_{\mathrm{eff}}=29000\mathrm{K}`$, $`\mathrm{log}g=3.0`$, $`R=30R_{}`$, $`\dot{M}=4\times 10^6M_{}\mathrm{yr}^1`$, $`v_{\mathrm{}}=1500\mathrm{km}/\mathrm{s}`$. (Note that the consistent treatment of the hydrodynamics in our procedure will be described in a forthcoming paper (Pauldrach and Hoffmann 2002); the method is also illustrated in these proceedings by Hoffmann and Pauldrach.) This implies a stellar mass of $`33M_{}`$ and a spectroscopic distance of 1.2 kpc, in good agreement with the distance estimate of 1 kpc (Scuderi et al. 1998). Thus we show that our current models produce satisfactory results for massive Population I stars. Can we apply the same procedure to CSPNs?

## 3. The relation between wind-momentum loss rate and luminosity

The radiatively driven wind theory predicts, for fixed abundances, a simple relation between the quantity $`\dot{M}v_{\mathrm{}}`$, which has the dimensions of a momentum loss rate, and the stellar luminosity:

$$\dot{M}v_{\mathrm{}}R^{0.5}L^{(1/\alpha )}$$

where $`\alpha `$, the power law exponent of the line strength distribution function, is $`2/3`$ (slightly dependent on temperature and metallicity; see, for example, Puls et al. 1996). It is practical to plot the log of $`\dot{M}v_{\mathrm{}}R^{0.5}`$ as a function of $`\mathrm{log}L`$. In this kind of plot the theory predicts, in first approximation, a linear relation, which is indeed followed by all kinds of massive hot stars, as shown in Figure 2.

An initial attempt to verify if CSPNs follow the wind-momentum–luminosity relation was partly successful (see Figure 3 in Kudritzki et al. 1997 and also our Figure 2). In that paper $`Q`$-values (a quantity relating mass loss rate and stellar radius, $`Q\dot{M}(Rv_{\mathrm{}})^{3/2}`$) were derived from observed H$`\alpha `$ profiles, and the stellar masses were derived from $`T_{\mathrm{eff}}`$ and $`\mathrm{log}g`$, using post-AGB tracks plotted in the $`\mathrm{log}g`$$`\mathrm{log}T_{\mathrm{eff}}`$ diagram. The stellar radii (and thus, mass loss rates) and luminosities were then obtained from the masses and the post-AGB mass–luminosity relation. The CSPNs were found to be at the expected position along the wind-momentum–luminosity relation, indicating a qualitatively successful prediction by the theory of radiatively driven winds. However, the situation was not satisfactory because there appeared to be a large dispersion in wind strengths at a given luminosity (strong-winded and weak-winded CSPNs) and some of the CSPN masses and luminosities were very high ($`M>0.8M_{}`$), in contradiction with theoretical post-AGB evolutionary speeds

Thus, at that point we had a qualitative positive result, namely that in principle the CSPN winds obey the same physics as the massive O star winds; but we also had some unsolved problems which we now want to rediscuss using the improved model atmospheres.

As a first step we have used our models to calculate the terminal velocities and mass loss rates for a grid of stars following the current theoretical post-AGB evolutionary tracks (see, for instance, Blöcker 1995); the resulting wind momenta are also plotted in Figure 2. The numerical models do nicely follow the expected theoretical wind-momentum–luminosity relation, although the spread is much smaller than that in the values derived by Kudritzki et al. Furthermore, the location of the observed sample in the diagram indicates masses between $`0.6`$ and $`0.95M_{}`$, with a clear absence of CSPNs with masses below $`0.6M_{}`$.

This result is entirely unexpected from the standpoint of current evolutionary theory. To understand this discrepancy, we must compare the relations of the individual dynamic quantities, $`v_{\mathrm{}}`$ and $`\dot{M}`$. This is done in Figures 4 and 4: Figure 4 shows our calculated terminal velocities and the observed values, Figure 4 shows our computed mass loss rates and those derived by Kudritzki et al. for their sample.

Here, too, a fundamental discrepancy immediately becomes obvious: whereas the positions of the observations in the diagram showing the terminal velocities cluster at rather small CSPN masses (between $`0.5`$ and $`0.6M_{}`$), their mass loss rates point to a majority of masses above $`0.7M_{}`$.

A detailed look at the positions of individual CSPNs in the plots reveals even more alarming discrepancies. Take, for example, He 2-131. Its terminal velocity would indicate a mass of about $`0.6M_{}`$ (circle 1A in Figure 4). But this mass is completely irreconcilable with its mass loss rate: it is found not at the position labelled 1A in Figure 4, but at 1B, with $`\dot{M}`$ a factor of hundred higher, suggesting a mass of above $`0.94M_{}`$! The reverse is true for NGC 2392. Its terminal velocity points to a mass of about $`0.9M_{}`$ (circle 2A in Figure 4), but its observed mass loss rate is much too small for this mass (circle 2B in Figure 4), indicating a mass of approximately $`0.6M_{}`$.

We still face a problem if we take the mass loss rate determinations of Kudritzki et al. to be correct: then our calculations would place these two stars at the positions labelled 1B and 2B in Figure 4 – with terminal velocities differing by a factor of 2 to 3. But this is clearly ruled out by the observations. ($`v_{\mathrm{}}`$ is a directly measurable quantity.)

We are therefore left with the conclusion that the analysis of Kudritzki et al. is at odds with the mass loss rates computed by our models based on the post-AGB evolutionary tracks. To determine whether the reason for this discrepancy lies with the evolutionary tracks on the one hand or our hydrodynamical models and the analysis by Kudritzki et al. on the other requires further observational evidence. This is given to us by the UV spectra of the stars.

## 4. Examples: UV analysis of the CSPNs He 2-131 and NGC 2392

Figure 5 (top left) shows the synthetic UV spectrum of the model corresponding to 1A in Figures 4 and 4. It is clearly incompatible with the observed spectrum of He 2-131 (middle), since its mass loss rate is obviously too small, as evidenced by the presence of almost only purely photospheric lines hardly influenced by the thin wind – indicating that this CSPN must have a much larger luminosity, because $`L`$ is the major factor determining the mass loss rate. Thus, we have calculated a series of models with increasing luminosity – and therefore increasing mass loss rate – (at the same time adjusting the mass to keep the terminal velocity at its observed value) to see whether one of these models could reproduce the numerous strongly wind-contaminated iron lines observed especially in between $`1500`$ and $`1700\mathrm{\AA }`$. Indeed, a model which due to its luminosity yields approximately the high mass loss rate of 1B gives a much improved fit – see Figure 5, bottom left. (The parameters of this star are given in Table 1.)

The situation is reversed with NGC 2392. The synthetic spectrum of model 2A is incompatible with the observed UV spectrum (Figure 5, right top and middle), since it produces many strongly wind-contaminated lines, which are, however, not observed. Instead, almost only photospheric lines are produced by the star. Again the problem is the luminosity, which in this case is much too high. Decreasing the luminosity and thus the mass loss rate yields a model with a much better agreement with the observed spectrum (Figure 5, bottom right).

What does this mean for the derived stellar parameters, which by virtue of this analysis now have the status of observed quantities? Let us consider first the weak-winded CSPN, NGC 2392. We determine a $`T_{\mathrm{eff}}`$ of 40000 K from the ionization equilibrium of Fe ions in the stellar UV spectrum, not too different from the value obtained by the ionization equilibrium of He i and He ii (absorption lines in the optical stellar spectrum). The very low terminal velocity in the wind of $`400\mathrm{km}\mathrm{s}^1`$ leads, together with a decreased luminosity (in order to reduce the predicted mass loss rate until the predicted and observed spectra are in good agreement), to a small radius. From this radius ($`1.5R_{}`$) and $`v_{\mathrm{}}`$ we get a stellar mass of only $`0.41M_{}`$, a value much smaller than if we believe in the classical mass–luminosity relation – a high mass of $`0.9M_{}`$ was the result found by Kudritzki et al. (But due to the smaller $`R`$ and $`T_{\mathrm{eff}}`$ our luminosity is also smaller.) We remark at this point that our error in the mass is extremely small ($`0.1M_{}`$) due to the sensitive dependence on $`v_{\mathrm{}}`$ and the small error in this value ($`10\%`$). Furthermore, our predicted values of $`v_{\mathrm{}}`$ are in agreement within 10% with the observed values for the case of massive O stars (cf. Hoffmann and Pauldrach, these proceedings). Thus the systematic error is almost negligible.

Next we consider the central star of He 2-131. In this case the terminal velocity of $`500\mathrm{km}\mathrm{s}^1`$ ($`T_{\mathrm{eff}}=33000\mathrm{K}`$) would appear to suggest, according to the classical post-AGB mass–luminosity relation, a stellar mass of about $`0.6M_{}`$ (cf. Figure 4). However, the wind features observed in the UV spectrum forced us to increase the stellar $`R`$ and $`L`$, which in turn increased $`\dot{M}`$ until a good fit was obtained. From the corresponding large radius – $`5.5R_{}`$ – and $`v_{\mathrm{}}`$ we derive a stellar mass of $`1.39M_{}`$, a value very close to the Chandrasekhar mass of relevance for type Ia supernovae. Thus, in this case the resulting mass is even more extreme than the value of $`0.9M_{}`$ obtained by Kudritzki et al.

## 5. Interpretation of CSPN winds: results and discussion

Table 1 shows the result of applying the method and UV analysis described in the previous section to eight CSPNs, including the two exemplary objects above. The problem which Table 1 points to is obvious: according to the post-AGB evolutionary timescales we should not find so many extremely luminous CSPNs, because according to this theory they are expected to fade very quickly. We shall note, however, that the sample of objects chosen here is most likely not a representative one.

Figure 6 shows the relation between stellar mass and luminosity obtained from our model atmosphere analyses, in comparison with the mass–luminosity relation of the evolutionary tracks, represented by the values from Kudritzki et al. 1997. This plot shows that the problem already indicated is indeed very disturbing: the derived masses and luminosities do not agree with the classical post-AGB mass–luminosity relation. There is a very large spread in masses, between 0.4 and 1.4 $`M_{}`$, and there is no well-defined relation between CSPN mass and luminosity.

In Figure 7 we now again show the wind-momentum–luminosity relation for both massive hot stars and CSPNs, but this time based on the parameters derived in our analysis. Our models give wind momenta of the right order of magnitude and within the expected luminosity range (there may be too many CSPNs at $`\mathrm{log}L/L_{}>4`$, but not so many as in Kudritzki et al. 1997). The CSPNs are found along the extrapolation of the wind-momentum–luminosity relation defined by the massive hot stars, and the CSPNs show a smaller dispersion, i. e., a tighter correlation of wind-momentum with luminosity, than was the case in Kudritzki et al. (1997). And, most important, this was achieved by fitting the multitude of diagnostic features in the CSPN UV spectra by means of up-to-date hydrodynamically consistent models.

How then does our wind-momentum–luminosity relation compare to that found by Kudritzki et al.? The answer is, quite favorably. If we drop the assumption made by Kudritzki et al. that the stars obey the theoretical post-AGB mass–luminosity relation, and instead scale their mass loss rates to our radii222 Additionally allowing for their different effective temperatures by requiring that the observed visual flux ($`R^2T_{\mathrm{eff}}`$) stay constant. – keeping $`Q`$, the real observational quantity, fixed – then their wind momenta match ours to within about a factor of two. Furthermore, their sample with the radii thus scaled now also shows a much tighter correlation of the wind momentum to luminosity than before (see Figure 7).

All this is strong evidence that all these winds are radiatively driven. It would be extremely difficult to explain the wind-momentum–luminosity relation if there was another mechanism driving the winds.

What makes the CSPN mass discrepancy problem found most intriguing is the fact that the same model atmospheres, with exactly the same physics, obviously work also perfectly well for massive O stars, as we have shown. We are reluctant to conclude that there is something basic we do not understand about either the winds and photospheres of O-type stars, or how to produce CSPNs and what their internal structure is. Both alternatives are difficult to believe. But we need to explain why we obtain such a spread in the masses.

We cannot offer a fair solution to this paradox; all we can do right now is to present our method, the results obtained from our analysis, and the corresponding problem in the clearest possible way, which is usually the first step along the road that leads to the solution.

## References

Blöcker T., A&A 299, 755 (1995)

Hoffmann T. L., Pauldrach A. W. A., these proceedings

Kudritzki R.-P., Méndez R. H., Puls J., McCarthy, J. K., in IAU Symp. 180, Planetary Nebulae, eds. H. J. Habing & H. J. G. L. M. Lamers, p. 64 (1997)

Méndez R. H., Kudritzki R.-P., Herrero A., A&A 260, 329 (1992)

Pauldrach A. W. A., Puls J., Kudritzki R.-P., et al., A&A 207, 123 (1988)

Pauldrach A. W. A., Hoffmann T. L., Lennon, M., A&A 375, 161 (2001)

Pauldrach A. W. A., Hoffmann T. L., A&A in press (2002)

Puls J., Kudritzki R.-P., Herrero A., et al., A&A 305, 171 (1996)

Scuderi S., Panagia N., Stanghellini C., et al., A&A 332, 251 (1998)


text:
# Partonic Energy Loss and the Drell-Yan Process

## Abstract

We examine the current status of the extraction of the rate of partonic energy loss in nuclei from $`A`$ dependent data. The advantages and difficulties of using the Drell-Yan process to measure the energy loss of a parton traversing a cold nuclear medium are discussed. The prospects of using relatively low energy proton beams for a definitive measurement of partonic energy loss are presented.

Characterizing the energy loss of partons in nuclear matter has been a matter of considerable theoretical  and phenomenological  interest over the past few years. In spite of this activity there is no agreed to value for partonic $`dE/dx`$ (GeV/fm) in either cold or hot nuclear matter, though there is general agreement that the rate of energy loss in the latter exceeds that of the former.

Much of the problem stems from the inability to directly measure this energy loss, and the lack of a common agreement of the processes and mechanisms to be included in specifying the energy loss. The purpose of this communication is to briefly present the current status of the problem, present a clear idea of quark energy loss in the nuclear medium, and delineate an experimental program capable of determining the parameters to adequately characterize quark energy loss in cold nuclei. We do not deal with the issue of energy loss in relativistic heavy ion collisions (hot nuclear matter), because the complexity of these reactions makes quantitative extraction of partonic energy loss highly uncertain with present knowledge. We will further adhere to the conventional assumption that the rate of energy loss is independent of parton momentum at sufficiently high energies , though this assumption has been questioned .

Table I demonstrates the lack of agreement among the extracted values for quark energy loss that exist at the present time. Admittedly, much of the difficulty may be matters of definition. There is variability in the implicit definitions of what is included in the energy loss and the correct path length over which the energy loss takes place. To reduce some of the confusion in Table I, an average path length is employed using the conventional value, $`L=3/4(1.2A^{1/3})`$ fm. There are more elaborate treatments  of the nuclear path length that yield appreciably shorter path lengths and hence higher rates of energy loss.

Doubtless some of the discrepancies in Table I would be reduced if consistent definitions, analysis and proper attention to systematic errors were applied. The nuclear dependence of fragmentation appears to us to be too uncertain at the present time to permit a reliable extraction of partonic energy loss. The $`A`$ dependence of the $`p_T`$ broadening of the $`pA`$ Drell-Yan (D-Y) yield may be an excellent measure of the incident quark energy loss due to gluon radiation but the small value observed for the broadening indicates that non-perturbative physics is involved . The above procedures will not take account of the significant energy loss ($`1`$ GeV/fm) due to the color confining interaction which is operative even in vacuum. The same criticism may be levied against any direct comparison of jet energy loss in lepton-nucleus to lepton-proton DIS.

We believe that the best measure of the quark energy loss in cold nuclear matter is the nuclear dependence of the $`pA`$ D-Y cross section. Viewed in the infinite momentum frame the energy loss of the incident quark is terminated within the nuclear medium by the $`q\overline{q}`$ annihilation process. This approach has been taken by refs.  shown in Table I. The extracted energy losses are seen to be discrepant and are discussed below. First, for illustrative purposes a simple description of how nuclear energy loss affects the measured D-Y yield is presented. The $`pA`$, leading-order D-Y cross section per nucleon is given by

$`{\displaystyle \frac{d\sigma _{pA}^{DY}}{dx_1dx_2}}={\displaystyle \frac{4\pi \alpha _{em}^2}{9s}}{\displaystyle \frac{x_1x_2}{x_1+x_2}}{\displaystyle \underset{i}{}}e_i^2[q_i^p(x_1)\overline{q}_i^A(x_2)+\overline{q}_i^p(x_1)q_i^A(x_2)].`$ (1)
The fractional momenta $`x_1`$ and $`x_2`$ of the partons of flavor $`i`$ in the beam and target respectively are determined from the mass of the lepton pair, $`M^2=sx_1x_2`$, combined with the c.m. longitudinal momentum ($`p_l`$) of the pair $`x_1x_2=2p_l/\sqrt{s}`$. For large $`x_1`$ ($`x_1>0.5`$) and for an interval in $`x_2`$ where $`\overline{q}^A(x_2)/\overline{q}^D(x_2)1`$, which occurs for $`0.1<x_2<0.3`$ as shown in , the ratio of the D-Y yields from nucleus $`A`$ to that from deuterium can be expressed as

$`({\displaystyle \frac{d\sigma _{pA}^{DY}}{dx_1dx_2}})/({\displaystyle \frac{d\sigma _{pD}^{DY}}{dx_1dx_2}}){\displaystyle \frac{q_u^p(x_1^A)}{q_u^p(x_1^D)}}{\displaystyle \frac{(1x_1^A)^3}{(1x_1^D)^3}}.`$ (2)
The cubic dependence on $`(1x_1)`$ comes from quark counting rules which are approximately correct as $`x_11`$. The energy of the incident quark is equal to $`E_px_1`$. Assuming that the incident quark experiences an energy loss per unit length, $`\alpha `$, in nuclear matter, and that there is no such energy loss in deuterium, then

$`x_1^D=x_1^p,x_1^A=x_1^p+\alpha {\displaystyle \frac{L_A}{E_p}},`$ (3)
where $`L_A`$ is the average path length of the incident quark in the nucleus $`A`$. Because of the energy loss, $`x_1^A`$ must originate from larger values of $`x_1^p`$ than is the case for deuterium. Eq. 2 can then be expressed to lowest order in $`L_A/E_p`$ as

$`({\displaystyle \frac{d\sigma _{pA}^{DY}}{dx_1dx_2}})/({\displaystyle \frac{d\sigma _{pD}^{DY}}{dx_1dx_2}}){\displaystyle \frac{(1x_1^p\alpha \frac{L_A}{E_p})^3}{(1x_1^p)^3}}1{\displaystyle \frac{3\alpha L_A}{E_p(1x_1^p)}}.`$ (4)
Eq. 4 shows that the effects of partonic energy loss are the largest at lower proton energies and at larger $`x_1`$. Under readily achievable conditions with $`A=200`$, $`E_p=120`$ GeV and $`x_1=0.7`$, Eq. 4 yields 0.78 for $`\alpha `$= 0.5 GeV/fm, a readily measurable suppression due to partonic energy loss. Of course, for a more quantitative analysis the approximations used in this paragraph should not be employed. The features presented in Eq. 4 remain when measured parton distributions are employed.

Why then, is there such difficulty in extracting a reliable value for nuclear partonic energy loss? First, there has not been experimental activity focused on elaborating partonic energy loss. The available data come from a variety of sources  that focused on other aspects of nuclear dependence in the D-Y process. There is the additional serious problem of separating the effects of shadowing from those of energy loss. We recall that nuclear shadowing effect was observed in Deep Inelastic Scattering for $`x`$ less than $``$0.07. This effect is interpreted in the infinite momentum frame as a reduction of parton density at small $`x`$ due to recombination of quarks and gluons from different nucleons. In the target frame, nuclear shadowing occurs when the coherence length, $`L_c1/(2m_Nx_2)`$, grows larger than the average distance between nucleons. Both the energy loss and the shadowing effects lead to an $`A`$-dependent reduction of the D-Y cross section. In principle, they are readily separable as shadowing reduces the yield at small values of $`x_2`$ while the effects of energy loss are expected to be most pronounced for large values of $`x_1`$. However, the acceptance of typical fixed-target detectors are biased to favor small values of $`x_2`$ in conjunction with large values of $`x_1`$. References  in Table I represent different analyses of the nuclear dependence of $`E_p=800`$ GeV, $`pA`$ D-Y yield. Reference  used $`A`$-dependent data from FNAL E866 which had had a large fraction of its data at $`x_2<0.05`$ and thus required careful attention to nuclear shadowing. The shadowing corrections employed a phenomenological shadowing scheme due to Eskola et al.  which had determined the shadowing of sea quarks using the FNAL E772 data with no account of the effects of parton energy loss. Such an analysis naturally produces the small energy loss shown in Table I.

References  analyzed data from FNAL E772  and FNAL E866. They utilized a very different reference frame and prescription for calculating shadowing than used in . Specifically, the nuclear rest frame was employed so the lepton pair originates from a heavy photon bremsstrahlunged off the incident fast charged quark. This description has appreciably less shadowing and consequently requires greater parton energy loss to account for the reduced D-Y yield with increasing $`A`$. Reference  has a complete description of this analysis procedure and demonstrates that it correctly accounts for shadowing in DIS. The energy loss per unit length reported in  is more than twice as large as that shown in Table I because in  the parton’s path length, $`L_{A=184}`$ in Tungsten is only 2.1 fm. This rather short parton path length is due to the requirement in  that the incident proton needs to travel a certain distance (mean-free-path $``$ 2.5 fm) before an interaction can liberate the high-$`x`$ parton. The results shown in Table I use the commonly defined parton path length to make the results more comparable. In , the incident proton energy is high ($`E_p`$ = 800 GeV) so that fractional energy loss is small, making the extracted energy loss strongly dependent on the applied shadowing corrections.

The strongest evidence for a small rate of parton energy loss would appear to come from the nuclear dependence of $`\pi ^{}N`$ D-Y production. Two experiments, one at 225 GeV/c  and one at 150 GeV/c  reveal little in the way of nuclear dependence. An analysis  of the data from  is shown in Table I. The value obtained, $`dE/dx=0.20\pm 0.15`$ GeV/fm, is in serious disagreement with . Apart from explanations involving the low statistics of  and incomplete knowledge of the pion’s parton distribution, the difference might lie with a smaller $`L_A`$ for partons from the pion. The $`\pi `$-N cross section is only $`60\%`$ that of N-N resulting in a longer nuclear mean-free-path for the pion and a corresponding shorter $`L_A`$ for the partons.

Thus, due to the multiplicity of significant physical effects and the sparse database available, there is not yet a reliable value for the partonic energy loss in cold nuclei. To remedy this situation we propose a focused set of $`pA`$, D-Y measurements at $`E_p100`$ GeV. The lower energy enhances the effects of energy loss, accesses a region of $`x_2`$ where shadowing is negligible, and using a proton beam allows the accumulation of adequate statistics in a reasonable time. The use of several targets of differing $`A`$ will be necessary in order to extract the effects of the appropriate parton path length, as well as the relative contributions of linear and quadratic dependence of the energy loss on path length.

To assess this issue more quantitatively we have examined the projected sensitivity for D-Y experiments using parton distributions fit to the world data and have incorporated nuclear effects in addition to parton energy loss. In particular, we consider using the proton beams from both the 120 GeV Fermilab Main Ring Injector and the future 50 GeV Japan Hadron Facility. Fig. 1 shows how parton energy loss would affect the $`(p+W)/(p+d)`$ per nucleon D-Y yields at these beam energies. In this calculation, we use the kinematic acceptance of a dimuon spectrometer proposed for Fermilab Experiment 906  and a similar spectrometer considered  at the Japan Hadron Facility. Only D-Y muon pairs with $`M>4.2`$ GeV are used to avoid contamination from charmonium decays. The coverages in $`x_2`$ are $`0.1<x_2<0.35`$ and $`0.2<x_2<0.5`$, respectively, for beam energies of 120 GeV and 50 GeV, well outside the small-$`x_2`$ region of nuclear shadowing. The parton distribution functions are taken from  and their modifications in the nuclear medium are from . Leading-order D-Y formula given in Eq. 1 was used in this calculation. The effect of the next-to-leading order diagrams is to introduce an overall normalization factor which is expected to cancel in the D-Y cross section ratios. The average path length in nucleus was taken as $`3/4(1.2A^{1/3})`$. The various curves in Fig. 1 correspond to a partonic energy loss, $`dE/dx`$, of 0.0, 0.1, 0.25 and 0.5 GeV/fm, respectively. As expected, these curves clearly demonstrate the great sensitivity of the D-Y cross section ratios to partonic energy loss when the incident proton energy is as low as shown in these examples. The projected statistical errors for a 60-day run on each target are also shown in Fig. 1. A partonic energy loss as small as 0.1 GeV/fm is seen to be readily measurable.

$`A`$-dependent D-Y data can further determine whether the energy loss is linear or quadratic with the path length. This is illustrated in Fig. 2 where D-Y cross section ratios are shown as a function of $`A`$. The solid circles correspond to a linear energy loss of 0.25 GeV/fm. The open squares correspond to the case where the energy loss is proportional to the path length squared and matched to give the same energy loss when $`A=184`$ (Tungsten target). Fig. 2 shows that one can easily distinguish between $`L`$ and $`L^2`$ dependence of the energy loss even when it is as small as 0.25 GeV/fm.

Thus it appears that many of the presently confusing issues of partonic energy loss can be resolved via focused measurement of the p-A dependence of the Drell-Yan cross section. A comparison of the yields from several nuclear targets relative to that from deuterium should determine the amount of energy loss and the appropriate path length, thereby firmly establishing the rate of energy loss. In addition to solving an interesting problem in its own right, these measurements will have a direct impact on nuclear shadowing and on the interpretation of many hard scattering processes in nuclei.