text:
# Observational Tests and Predictive Stellar Evolution

## 1 INTRODUCTION

Prior to any rigorous investigation of thermonuclear yields, mixing, rotation, mass loss, and other complex phenomena in the evolution of stars, it is necessary to insure that the methodology used can reproduce observations at the current state of the art. Detached, double-lined eclipsing binaries provide the most accurate source of information on stellar masses and radii, and as such provide a crucial test for models of stellar evolution (Andersen, 1991).

The predictive power of such simulations depends upon the extent to which the validity of their oversimplifications can be tested. Phenomenology is particularly pernicious in that good tests, which could give rise to falsification of the models and thus to progress, can be nullified by parameter adjustment (“calibration”). Similarly, the inclusion of new processes in the simulations is often contentious in that there may be several candidates put forward as the cause of the puzzling data. The treatment of the boundaries of convection zones (“overshooting”) is ripe for re-examination in terms of the underlying hydrodynamics (Asida & Arnett, 2000). To aid this development, we wish to understand just how well a standard convective treatment can do (i.e., we turn overshooting off). A more physically sound procedure will be presented subsequently. Parameterized overshooting has been widely discussed (see Maeder (1975) for an early discussion, and Schröder et al. (1997) and Schaller et al. (1992); Bressan et al. (1993); Dominguez et al. (1999) for a recent ones with extensive references). Rotation also may cause mixing (Meynet & Maeder, 2000). Within the context of a plasma, rotation and convection may generate magnetic fields, which by their buoyancy and angular momentum transport may provide additional causes for mixing.

This paper represents a first step toward testing our extensively revised stellar evolution and hydrodynamics code, TYCHO. Our goals are (1) the understanding the predictive capability of stellar evolution theory by the critical re-evaluation of its assumptions, and of its underlying basis in observations and in laboratory data, (2) the examination of the tricky problem of mixing in stars to help design experiments (Remington et al., 1999) and numerical tests (Asida & Arnett, 2000), and (3) the development of an open source, publicly available stellar evolution code with modern capabilities for community use.

### 1.1 Choice of Binaries

The most comprehensive list of binary systems with accurately measured masses and radii is given in the review by Andersen (1991). A subset of the original binaries was chosen for this exploratory effort. The systems with the smallest uncertainties were picked such that the range of masses from 23 to 1.1 $`M_{}`$ was well sampled. The upper mass limit is the largest mass present in the data, while the lower limit is safely above the point at which the equation of state used in the modeling becomes inaccurate. This is primarily due to approximate treatment of Coulomb contributions to the pressure; we use only the weak screening limit for the plasma. The Coulomb interaction leads to a negative pressure correction of $``$ 8% in the outer part of the convection region and $``$ 1% in the core for a star of 1 $`M_{}`$ (Däppen & Nayfonov, 2000), and is less important for more massive stars. We obtain a comparable correction to the central pressure for a solar model (-1.7%).

Aside from the exclusion of stars of still lower mass, no bias was applied in the selection process. For stars of $`1\mathrm{M}_{}`$ or less, the issue of possible overshooting in the convective core is moot because their cores are radiative. Also, because they rotate rapidly only for a brief part of their lives, rotational mixing is expected to be less than for our selected stars. The binaries used in the study along with their fundamental parameters are presented in Table 1. Ribas et al. (2000) have revised the Andersen (1991) temperature estimates, and revised estimates for the masses of EM Car (Stickland, Lloyd, & Corcoran, 1994) and CW Cep (Stickland, Koch & Pfeiffer, 1992) have been presented. Latham et al. (1996) have revised the parameters for DM Vir. All these changes are incorporated into Table 1.

### 1.2 The Mixing Length

If the purpose of a stellar evolution code is to make testable predictions of the behavior of stars, then the adjustment of the mixing length by fitting stellar data is repugnant. Alternatives are to constrain it by experiment or by simulation. At present we know of no definitive experimental results which determine the mixing length appropriate to stars, although a variety of experiments do test other aspects of stellar hydrodynamics and the codes used to simulate them (Remington et al., 1999). However, it is becoming possible to simulate turbulent, compressible convection with sufficient realism to constrain the range allowed for the mixing length (Rosenthal et al., 1999; Porter & Woodward, 1994, 2000; Elliott, Miesch, & Toomre, 2000; Asida, 1998). Porter & Woodward (2000) find a mixing length $`\alpha _{ML}=2.68`$ in units of pressure scale height; this is based upon simulations having mesh resolutions as high as $`512\times 512\times 256`$ and corresponding Rayleigh numbers as high as $`3.3\times 10^{12}`$. Rosenthal et al. (1999) based their work on resolutions up to $`253\times 253\times 163`$, but with a more realistic treatment of radiative transfer and ionization. They were able to synthesize the line profile of FeII $`\lambda `$5414 which compared well with that observed. Their results agreed with standard 1-D models, although they suggest that this might be “the right result for the wrong reason.” Standard models use $`\alpha _{ML}=1.5`$ to $`2`$, which is significantly smaller than the value of Porter & Woodward (2000). For a red giant, Asida (1998) inferred $`\alpha _{ML}=1.6`$, based upon 2-D simulations but with fairly realistic microphysics.

Canuto & Mazzitelli (1991, 1992) have proposed a serious model to replace mixing length theory; this has had the salutary effect of shifting the debate to the physics of convection and away from the best choice of mixing length. Finally, Asida & Arnett (2000) have shown from 2-D simulations that the underlying physical picture for stellar convection is incomplete, even in the deep interior where the complication of radiative transfer is minor.

We simply choose $`\alpha _{ML}=1.6`$, and look forward to the convergence of these efforts to provide a convection algorithm which is independent of stellar evolutionary calibrations.

### 1.3 The TYCHO Code

The evolutionary sequences were produced with the TYCHO stellar evolution code. The code was originally developed for one dimensional (1D) hydrodynamics of the late stages of stellar evolution and core collapse (Arnett, 1996). It is being completely rewritten as a general purpose, open source code for stellar evolution and hydrodynamics. The present version is written in structured FORTRAN77 and is targeted for Linux machines. It has been successfully ported to SunOS and SGI IRIX operating systems. It has extensive online graphics using PGPLOT, an open source package written by T. J. Pearson (tjp@astro.caltech.edu). A library of analysis programs is being built (modules for apsidal motion, pulsational instability, reaction network links, and history of mass loss are now available). The code is being put under source code control to allow versioning (this will allow particular versions of the code—for example the one used in this paper, to be resurrected accurately at later times), and to improve the reproducibility of results.

Knowledge of the radiative opacity of plasma at stellar conditions has changed qualitatively in the past decade. Historically, solar and stellar atmospheres provided much of the empirical data on hot plasmas. For example, Kurucz (1991) tabulates a range of effective temperatures from 2,000 to 200,000 K. These temperatures and the corresponding (low) densities characterized what was directly observable. Higher temperatures could be found at lower densities (in non-LTE) or indirectly inferred by use of theoretical models. Terrestrial tests involved explosions which were difficult to quantify with adequate precision to determine opacity. Measurements of opacity in a well characterized, hot, dense, laser-produced plasma have become possible (Perry et al., 1991; Springer et al., 1992; Mostovych et al., 1995; Davidson et al., 2000). The first experiments to simultaneously quantify temperature and density with good precision (Perry et al., 1991; Springer et al., 1992) involved temperatures $`T7\times 10^5\mathrm{K}`$ and densities $`\rho 2\times 10^2\mathrm{g}\mathrm{cm}^3`$, which are directly relevant to stellar evolution and to apsidal motion. The range of conditions which are experimentally accessible is expanding with the development of new instruments and techniques. Not only can direct measurements be made, but complex and sophisticated theories of the physical state of the plasma can be tested, giving more reliable extrapolations into conditions not yet experimentally accessible (Perry et al., 1996; Davidson et al., 2000). The conditions just quoted are encountered in stars of about $`1M_{}`$, and are important for stellar evolution (see Ch. 7, Arnett (1996)).

The opacities used here are from Iglesias & Rogers (1996) and Kurucz (1991), for a solar abundance pattern (Anders & Grevesse, 1989). The Iglesias & Rogers (1996) opacities were computed with 21 elemental species; Iglesias et al. (1995) have shown that the remaining elements are so rare as to have only a marginal effect on the Rosseland mean opacities (for solar relative abundances of the heavier elements). While the OPAL opacities were constructed for astrophysical use, the underlying experiments and theoretical models are determined by the inertial confinement fusion (ICF) community, reducing the danger of unconscious bias from astronomical puzzles leaking back into the construction of opacities. Extension of the opacity library to lower temperatures and lower entropies is planned.

TYCHO is designed to use an adaptable set of reaction networks; for these calculations, two networks were used. At higher temperatures ($`T10^7\mathrm{K}`$), an 80 element reaction network was solved. The reaction rates were from F. K. Thielemann (private communication); see also (Thielemann, Arnould, & Truran, 1988). For lower temperatures this was replaced by a 15 element network which was designed for deuterium, lithium, beryllium, and boron depletion. The reaction rates were from Caughlan & Fowler (1988). The match at the temperature boundary was sufficiently good as to require no smoothing. More recent compilations of nuclear rates are available (Rauscher & Thielemann, 2000; Angulo et al., 1999), but were not used here to simplify the comparison with previous work.

The outer boundary condition was determined by use of the Eddington approximation to a grey, plane parallel atmosphere, integrated in hydrostatic equilibrium inward to a fitting point for the interior. For the most extended model considered here, the ratio of mean free path to radius was $`\lambda /R5\times 10^3`$, so that spherical effects are negligible. The most vigorous mass loss considered was so mild that the ram pressure $`\rho v_{loss}^2`$ at the photosphere was $`10^6`$ of the total thermal pressure, which is consistent with the hydrostatic assumption. Such integrations were used to define the pressure and temperature at the fitting points, $`T_f(L,R)`$ and $`P_f(L,R)`$, for the stellar luminosity $`L`$ and radius $`R`$. Their derivatives with respect to stellar luminosity and radius were approximated by finite differences constructed between three such integrations (at $`L,R`$, $`L+\delta L,R`$, and $`L,R+\delta R`$). Typically, of order 200 to 400 steps were used in the envelope integration.

We used Schwarzchild convection as formulated by Kippenhahn & Wiegert (1990), and our treatment of convective overshooting was turned off.

Mass loss was included and based on the theory of Kudritzki et al. (1989) for $`T_{eff}7.5\times 10^3\mathrm{K}`$ and the empirical approach of de Jager, Nieuwenhuijzen, & van der Hucht (1988) for lower effective temperatures. R. Kudritzki kindly provided appropriate subroutines for the hotter regime. Even for EM Car, the most massive system in the list, the effects of mass loss were modest ($`0.6\mathrm{to}0.7\mathrm{M}_{}`$).

The equation of state was that discussed in Timmes & Arnett (1999), augmented by the solution of the ionization equilibrium equations for H, He, and a set of heavier elements scaled from the solar abundance pattern. Both the equation of state and the thermonuclear reaction rates are affected by coulomb properties of the plasma. Only weak screening was necessary here. Extension to include both weak and strong screening consistently in the equation of state and thermonuclear reaction rates is planned; previous versions of the code included strong screening as well.

Models were run for each mass, starting with a fully convective initial model on the Hayashi track and ending well beyond hydrogen depletion in the core. A more realistic approach would have been to form the stars by accretion ( A. G. W. Cameron, private communication; Norberg & Maeder (2000)). We justify our choice by its simplicity, and by noting that only the last stages of the pre-main sequence are relevant here, for which the two cases give similar results.

Zoning in the interior typically ranged from 200 to 500 zones. All runs had solar heavy element abundance (Anders & Grevesse, 1989) and used a ratio of mixing length to scale height of $`\alpha =1.6`$ for convection. This choice gave a reasonably good solar model when compared to Bahcall & Pinsonneault (1998) and Christensen-Dalsgaard (2000); inclusion of element settling by diffusion and adjustment of the helium abundance would give improved consistency for the present-day sun, but diffusion would have less time to operate in the more massive stars considered here. Rotational mixing was turned off.

### 1.4 Related Investigations

The stars are all at or near the main sequence, so that the possible list of citations is enormous; the efforts of the Padua group (Bressan et al., 1993), the Geneva group (Schaller et al., 1992), and the FRANEC group (Dominguez et al., 1999) have comparable input physics and form a useful context. We focus discussion on Pols et al. (1997b) and Ribas et al. (2000), who consider many of the same binaries, and Claret & Giménez (1993) who examined the apsidal motion test.

The largest differences in microphysics between Pols et al. (1997b) and this paper are our use of Iglesias & Rogers (1996) rather than Rogers & Iglesias (1992) opacities, a more realistic nuclear network, and a different approach to the equation of state, but none of this seems to be particularly significant for the issues here. We do include mass loss, but these effects are not large. Pols et al. (1997b) do not calculate the pre-main sequence evolution (which is relevant to several of our binaries). They define an overshooting parameter which is fixed by previous work on $`\zeta `$ Aurigae binaries (Schröder et al., 1997). They construct a grid of models in mass (0.5 to $`40\mathrm{M}_{}`$) and heavy element abundance ($`Z=0.01`$, 0.02 and 0.03), assuming $`X=0.763.0Z`$ and $`Y=0.24+2.0Z`$ for the abundances of hydrogen and helium. They minimize a $`\chi ^2`$ error estimator in four parameters: the masses $`M_A`$ and $`M_B`$, the age $`t`$ of the binary, and $`Z,`$ the heavy element abundance.

As a test of consistency for later evolution, we have reproduced the $`4`$ and $`8\mathrm{M}_{}`$ sequences of Pols et al. (1997a), which did not use overshooting. The notoriously sensitive blue loops were reproduced to graphical accuracy (their Figure 4) for the same input physics. The codes seem highly consistent.

Ribas et al. (2000) used the models of Claret (1995, 1997b) and Claret & Giménez (1995, 1998), which use Rogers & Iglesias (1992) opacities and a 14 isotope network and include overshooting and mass loss. They too interpolated in a grid, minimizing a $`\chi ^2`$ error estimator. This procedure was more complex than that used by Pols et al. (1997b), and need not be described here. Both heavy element abundance and helium abundance were freely varied.

Our strategy differs from both Pols et al. (1997b) and Ribas et al. (2000), which may provide a useful contrast. Here we are interested in isolating the possible inadequacy of the standard formulation of stellar evolution, so we avoid optimization of parameters as much as possible. By using (1) solar abundances and (2) the measured masses, we reduce the degrees of freedom, and hopefully make the possible flaws in our stellar evolution prescriptions easier to see. By the same token, our models should fit the data less well because we have not optimized abundances or masses. Mathematically, optimization will almost always improve the fit, but not necessarily for the correct reasons. However, the actual abundances may be different from our assumptions, and the masses do have error bars.

Claret & Giménez (1993) used Rogers & Iglesias (1992) opacities, solar abundances, a mixing length ratio $`\alpha _{ML}=1.5`$, and overshooting of $`\alpha _{OV}=0.2`$ pressure scale heights (that is, essentially the same physics as the models used by Ribas et al. (2000)), and computed structure constants for apsidal motion for seven of the binaries we consider (EM Car, CW Cep, QX Car, U Oph, $`\zeta `$ Phe, IQ Per, and PV Cas).

Detailed comparisons will appear in the discussion below.

## 2 FITTING MODELS TO BINARIES

The first step in comparing the binary data with the computations is the choice of the best models. This was done by examining a $`\chi ^2`$ quantity for each binary pair, defined by

$`\chi ^2`$ $`=`$ $`((\mathrm{log}L(m_A,t)\mathrm{log}L_A)/\sigma _{LA})^2`$ (1)
$`+`$ $`((\mathrm{log}L(m_B,t)\mathrm{log}L_B)/\sigma _{LB})^2`$
$`+`$ $`((\mathrm{log}R(m_A,t)\mathrm{log}R_A)/\sigma _{RA})^2`$
$`+`$ $`((\mathrm{log}R(m_B,t)\mathrm{log}R_B)/\sigma _{RB})^2,`$
where $`A`$ and $`B`$ denote the primary and the secondary star, respectively. Here $`L_A`$ and $`R_A`$ are the observationally determined luminosity and radius of the primary, with $`\sigma _{LA}`$ and $`\sigma _{RA}`$ being the observational errors in $`\mathrm{log}L_A`$ and in $`\mathrm{log}R_A`$. We convert the observational data for the radii to logarithmic form for consistency. Correspondingly, $`L(m_A,t)`$ and $`R(m_A,t)`$ are the luminosity and radius of the model. This $`\chi ^2`$ was evaluated by computing two evolutionary sequences, one for a star of mass $`m_A`$ and one for $`m_B`$, and storing selected results from each time step. Then these files were marched through, calculating $`\chi ^2`$ at consistent times ($`t_A=t_B`$ to a fraction of a time step, which was a relative error of a few percent at worst). The smallest $`\chi ^2`$ value determined which pair of models was optimum for that binary. Note that if the trajectories of both $`A`$ and $`B`$ graze the error boxes at the same time, $`\chi ^24`$. These error parameters along with the corresponding uncertainties from the observations are presented in Table 2.

Figure 1 displays the resulting $`\chi ^2`$ for each binary pair, in order of descending mean mass. The binaries fall into three separate groups: ten have excellent fits ( $`\chi ^2<4`$; EM Car, V478 Cyg, CW Cep, CV Vel, U Oph, PV Cas, RS Cha, PV Pup, DM Vir and UX Men), six are marginal ($`16\chi ^24`$; QX Car,$`\zeta `$ Phe, IQ Per, MY Cyg, EK Cep, and V1143 Cyg), and two are poor fits ( $`\chi ^2>16`$, denoted offscale in Figure 1; VV Pyx and AI Hya). The boundaries between these groups are indicated by vertical lines.

### 2.1 Global Aspects of the Errors

The weakness of a $`\chi ^2`$ measure is that it is most meaningful if the errors have a gaussian distribution around the mean (Press et al., (1992), chapters 14 and 15), which does not seem to be the case here. In particular, systematic shifts in the empirical data, due to new analyses, can give significant shifts in the error estimation. Ribas et al. (2000) have re-estimated the effective temperatures of 13 of the 18 binaries we have examined. Five (QX Car, U Oph, PV Cas, AI Hya, and RS Cha) were changed by more than twice the error estimates of either Ribas et al. (2000) or Andersen (1991). Further, Stickland, Koch & Pfeiffer (1992) and Stickland, Lloyd, & Corcoran (1994) have analyzed additional data (from IUE) and find masses of CW Cep and EM Car which lie beyond twice the error estimates. This is to be expected if the errors are dominated by systematic effects, and warns us to distrust all but our most robust inferences.

Because the fractional errors in mass and in radius are much more restrictive, it strongly supports the need for renewed efforts to pin down the effective temperatures of these stars. The choice of $`L`$ and $`R`$ rather than $`L`$ and $`T_{eff}`$ in our definition of $`\chi ^2`$ is significant: the smaller errors for $`R`$ make the $`\chi ^2`$ more discriminating. Pols et al. (1997b) use $`R`$ and $`T_{eff}`$ which has the slight advantage here of involving less propagation of observational errors, but because $`R`$ is much more precise than $`T_{eff}`$, the effect is small for the present data.

We have chosen to update the original data of Andersen (1991), incorporating the changes made by the Ribas et al. (2000) effective temperatures and the Stickland, Koch & Pfeiffer (1992) and Stickland, Lloyd, & Corcoran (1994) masses. We have used the new data for DM Vir (Latham et al., 1996). Our general conclusions are unaffected by which of these sets of data we use.

The comparison of observed and computed stars may be presented as an goodness of fit vector, which has the advantage of being directly representable in the HR diagram for the stars. The observed points with error bars are plotted along with an arrow indicating the distance and direction to the best model point. The way in which the models differ from the observations can then be taken in at a glance. Figure 2 shows the goodness of fit vectors, from the observed points (shown with error bars) to the best model star (chosen as described above). The largest discrepancy is the secondary of VV Pyx. Of the eight binaries (QX Car, $`\zeta `$ Phe, IQ Per, VV Pyx, AI Hya, EK Cep, MY Cyg, V1143) which have mediocre or poor fits, seven (QX Car is the exception) have at least one component lying in the range $`1.7<\mathrm{M}/\mathrm{M}_{}<2.6`$. Andersen, Nordström, & Clausen (1990) noticed similar behavior.

Figure 3 shows the luminosity differences between the models and the stars. The vertical axis is mass in solar units; binary components are connected by a line. The two binaries with $`\chi ^2>16`$ (VV Pyx and AI Hya) are denoted by crosses; they are poor fits and should be given little weight. Considering the best fits, $`\chi ^2<4`$ (solid squares), there is a dramatic trend: the highest mass models (for example, EM Car) are underluminous relative to the actual binaries, while the lower mass models are not.

Given the indications that the errors may be dominated by systematic effects, we approach a statistical discussion with caution. The two binaries which have $`\chi ^2>16`$ are eliminated from this statistical discussion on the basis that these fits are too poor to be meaningful. In principle, the mean errors could show a systematic shift in the models relative to the data, but because we choose an optimum pair of models, the choice masks any absolute shift. The error should reappear as a larger RMS difference. For luminosity, the first moment of the difference between model and stellar logarithmic luminosity is just the mean of this difference, which is $`0.017`$ in the base ten logarithm (the models are too dim by this amount). The shift is smaller than the RMS error of the observational data, which is 0.056. If there were a bad global mismatch, the RMS difference in “model minus star” would be much larger than the average error in the observations. However, the RMS difference between the models and stars is 0.054, which is almost the same as the observational error. The luminosity is basically a measure of the leakage time for radiation, which is dominated by the value of the opacity in the radiative regions. It samples the whole star, including the deep interior. In a global sense, our mean leakage rate seems correct to the level of the statistical and observational error.

The shifts in radius between the models and the stars are shown in Figure 4. The vertical axis and the symbols are the same as in the previous figure. The mean shift is 0.0053 in the logarithm (the models are too large by this small amount); the corresponding standard error of the Andersen data for which the fits are acceptable is 0.016, to be compared with an RMS differnce between models and stars of 0.014. Except for a few outlying cases, the distribution is fairly uniformly distributed around zero. If only the best fits (squares) are considered, a subtle trend might be inferred: 9 of 12 of the models above $`4\mathrm{M}_{}`$ have radii which are too large.

The corresponding mean shift in $`\mathrm{log}T_{eff}`$ is $`0.007`$ (the models are too cool by this amount). Again, this is small in comparison to the standard error of the observations ($`0.014`$). The corresponding RMS difference between models and stars is $`0.017`$. The effective temperature is a surface quantity, and is more sensitive to the outer layers which contain little of the stellar mass.

These numbers suggest that standard stellar evolutionary sequences of these stages should be able to produce luminosities, radii, and effective temperatures within 11, 3 and 4 percent, respectively, of good observational data. Otherwise, new physics is indicated. Because the standard stellar evolutionary models do this well, small “improvements” may contain no information. We will emphasize systematic trends, and those implications which emerge from several independent tests.

### 2.2 Massive Binaries

Figure 5 shows the evolution of the model stars in log luminosity and log radius, for EM Car, V478 Cyg, CW Cep, QX Car, CV Vel, and U Oph, corresponding to a mass range from 23 to 4.6 $`\mathrm{M}_{}`$. Except for QX Car ($`\chi ^2=11.3`$), models of these binaries have $`\chi ^2<4`$, and so represent good fits. The error bars are centered on the observed stars; the arrows point from them to the optimum models. Notice the the fits can be multivalued because the trajectories may pass through the error boxes multiple times. This is shown occurring first as the model descends from the pre-main sequence (pre-MS), and again during main sequence hydrogen burning. If the stellar masses are significantly different, this ambiguity is removed by the condition that both components have the same age. All the model stars are too dim (all the arrows point downward), a signal that the standard stellar evolution prescription is systematically wrong.

EM Car, being the most massive system, also has the most significant mass loss. The model evolutionary sequences are set by the choice of initial mass, but the observational constraint on mass is applied after the best fitting model is determined. This is an implicit function of the choice of initial mass, and iteration is required. Initial masses of 22.91 and 20.91 produce masses at fit of 22.25 and 20.12 $`\mathrm{M}_{}`$, respectively. This mass ratio of 0.904 is consistent with the observational value of $`0.910\pm 0.011`$ (Stickland, Lloyd, & Corcoran, 1994), and the masses agree with observation to within the estimated errors ($`\pm 0.3\mathrm{M}_{}`$).

However, even this loss is still small. A loss of $`0.7\mathrm{M}_{}`$ is about twice the uncertainty in mass determination, $`\pm 0.32\mathrm{M}_{}`$. Such a change in mass, since $`LM^4`$, corresponds to a shift in luminosity of $`\mathrm{\Delta }\mathrm{log}L0.05`$, to be compared to the observational error in luminosity of $`\mathrm{\Delta }\mathrm{log}L=0.1`$, which is still larger. This is due to the fact that effective temperature is less well determined than the radius. A concentrated effort to refine the effective temperature determinations for EM Car, V478 Cyg, CW Cep, QX Car, CV Vel, and U Oph would translate directly into much sharper constraints on massive star evolution. Ribas et al. (2000) have revised the effective temperature for QX Car (upward) by twice the quoted error, so that the inferred luminosity increases. Prior to this revision, the fit to QX Car was good ($`\chi ^2<4`$). This larger discrepancy for QX Car is in the same sense as noted for the other massive systems; the models are dimmer than the stars. Taking the larger masses from Andersen (1991), with or without mass loss, still gives good fits, and the models are still dimmer than the stars. The result seems robust.

Pachoulakis et al. (1996) have used high resolution spectral images obtained with the International Ultraviolet Explorer (IUE) to study the winds from CW Cephei (HD 218066). They place upper limits on the mass-loss rates of $`1.0\times 10^8\mathrm{M}_{}\mathrm{yr}^1`$ for the primary and $`0.32\times 10^8\mathrm{M}_{}\mathrm{yr}^1`$ for the secondary. The model masses start at 12.9 and 11.9 and decrease only to 12.8 and 11.88 respectively, to be compared to $`12.9\pm 0.1`$ and $`11.9\pm 0.1\mathrm{M}_{}`$ (Stickland, Koch & Pfeiffer, 1992). The mass loss predicted by the Kudritzki et al. (1989) theory for our models at the point of minimum error is $`0.66\times 10^8\mathrm{M}_{}\mathrm{yr}^1`$ for the primary, and $`0.43\times 10^8\mathrm{M}_{}\mathrm{yr}^1`$ for the secondary. If the upper limits were detections, this could be considered good agreement, considering the complexity of the problem of interpreting the system (Pachoulakis et al., 1996). The net loss of mass up to this point is no larger than the error in mass determination, $`\pm 0.1\mathrm{M}_{}`$. Because the mass loss rate is restricted by these observations to be at or below the value we use, the mass loss process should have no larger effect than we compute. Hence, the remaining discrepancy must come from some other effect.

Table 3 gives the instantaneous mass loss rates from the models, at the point of optimum fit, for the most massive binary systems. At lower masses, the mass loss rates are smaller still. Additional observational data on mass loss for these systems could prove crucial in clarifying the role of mass loss in stellar evolution.

Ribas et al. (2000) estimate ages for EM Car and CW Cep. Their procedure not only gives ages, but also error estimates for those ages. Our ages agree with theirs to within these errors, even though we use no overshooting and they do. It may be that the convective region in high mass stars is sufficiently large that the gross evolutionary properties of stars on the main sequence are not greatly affected by the overshooting correction. The understanding of the physics of overshooting is still too preliminary to do more than speculate on this issue.

Pols et al. (1997b) find acceptable fits for EM Car, V478 Cyg, CW Cep, QX Car, and U Oph (it is probable that QX Car would not have been a good fit with the revised effective temperatures), but with increasingly lower heavy element abundances with increasing mass. EM Car and V478 Cyg have fits at the limit of the heavy element abundance range. Ribas et al. (2000) find a similar effect: the heavy element abundances of their massive binaries are marginally smaller than those of the less massive ones. Their effect is not quite as obvious as in Pols et al. (1997b), perhaps because Ribas et al. (2000) do not force the helium abundance to correlate with heavy element abundance, and it fluctuates for these systems. The added degree of freedom may allow the fitting procedure to obscure the trend.

This behavior could be interpreted as a galactic evolutionary effect, which would be extremely interesting, but there is another possibility. The problem with the massive models (see Figure 2 and Figure 3) is that they are too dim. Lower heavy element abundance gives higher luminosity because of reduced opacity. The fitting algorithms, having little freedom for mass variation (thanks to the high quality of the data), must find lower heavy element abundance preferable, whether or not the heavy element abundances are actually smaller. More effective mixing, giving larger cores, also results in higher luminosities even if the abundances are unchanged. It is crucial to obtain spectroscopic information to decide the issue. Guinan et al. (2000) have recently examined V380 Cyg, which is a binary of disparate masses ($`11.1\pm 0.5\mathrm{M}_{}`$ and $`6.95\pm 0.25\mathrm{M}_{}`$) and evolutionary state. They conclude that more mixing is needed ($`\alpha _{OV}=0.6\pm 0.1`$). However, this system is complicated. Guinan et al. estimate that the system is approximately ten thousand years from Roche lobe overflow by the primary. Thus, conclusions based on this system should be made with caution.

The overshoot parameter $`\alpha _{OV}`$ may be a function of mass (at least). Similar behavior can also be seen in Claret & Giménez (1991), which finds different values of the best fit overshoot parameter for five different masses. Alternatively, rotational mixing might be increasingly effective for larger masses. Phenomenological prescriptions are valuable if they capture the essential physics of the phenomena; if fitted parameters turn out to be variable, a new formulation is needed. We have at least three causes for one effect; sorting this out is an interesting theoretical and observational challenge.

### 2.3 Intermediate Mass Binaries

Figure 6 shows $`\zeta `$ Phe, IQ Per, PV Cas, AI Hya, VV Pyx, and RS Cha, the group which has some of the most challenging binaries. The masses range from 3.9 to 1.1 $`\mathrm{M}_{}`$.

Both $`\zeta `$ Phe and IQ Per have a mass ratio significantly different from one: 0.65 and 0.49 respectively. Because the more massive components will evolve more rapidly, common age is a stringent constraint. In both cases, the error is dominated by the less massive component. For $`\zeta `$ Phe, the $`2.55\mathrm{M}_{}`$ secondary is brighter than the model; Pols et al. (1997a) have the same problem. Ribas et al. (2000) avoid it by using a lower heavy element abundance (0.013) and a higher helium abundance (0.29). The heavy element abundance might be tested by high resolution spectroscopy. For IQ Per, the $`1.74\mathrm{M}_{}`$ secondary is too blue; its evolutionary track never gets so hot. Pols et al. (1997b) attribute the difficulty in fitting $`\zeta `$ Phe and IQ Per to problems in determining $`T_{eff}`$ of the secondary (Andersen, 1991), which is much dimmer because of the large mass ratios of the components.

For $`\zeta `$ Phe and IQ Per, it is clear that much of our “difficulty” is due to the relatively small error bars; see Figure 2. Consequently, small changes may improve the $`\chi ^2`$ significantly, even if they do not correspond to the physics of the system. In the case of IQ Per, use of the Ribas et al. (2000) value of effective temperature improves the fit, compared to Pols et al. (1997a), as does adjustment of the abundances.

AI Hydrae is particularly interesting because the primary is fitted by a model which is swiftly evolving, so that catching it in such a stage is unlikely. Overshoot from the convective core would broaden the main sequence band and increase the age of the fast evolving primary, allowing the possibility of a fit with a more probable, slower evolutionary stage. This is consistent with the conclusions of Pols et al. (1997b), who find that AI Hya is the only binary for which the overshooting models give a greatly improved fit.

Some of the AI Hya behavior can be attributed to heavy element abundance effects. Both members are classified as peculiar metal line stars. The heavy element abundance of this system (from multi-color photometry) is 0.07 (Ribas et al., 2000), which is 3.5 times the value used in the models. The true interior composition cannot be this metal rich. We have examined a sequence which had a heavy element abundance of 0.03 rather than 0.02. This modest change gave a dramatic shift toward lower luminosity ($`\mathrm{\Delta }\mathrm{log}L0.06`$, which is three times the observational error) and cooler effective temperatures ($`\mathrm{\Delta }\mathrm{log}T_e0.027`$, which is also three times the observational error). However, if the heavy element abundance were high only near the surface of the star, the opacity effects would produce a shift to the red in the evolutionary tracks, which would bring the models much more in line with the observations.

VV Pyx has almost identical components, so that their coeval origin has almost no effect on the fit. They track the same path at essentially the same time. The fit is simply the point that the observational error box is most closely approached, and should be viewed with caution, especially as the models give a poor fit.

Only PV Cas and RS Cha have good fits ($`\chi ^2<4`$), and they are pre-main sequence (pre-MS). RS Cha has previously been suggested to be in a pre-main sequence stage (Mamajek, Lawson & Feigelson, 1999). The pre-MS identification would have important theoretical implications. If true, it implies that the error in the models occurs after the core convection has been established in these stars. In any case, convection is an interesting possible cause for the problem; these binaries have at least one component with convective core burning. PV Cas has sufficiently different masses to require us to examine the pre-MS fit seriously.

### 2.4 Is PV Cas Pre-Main Sequence?

Questions have been raised about the evolutionary status of PV Cas since Popper (1987). Previous attempts to fit the system to main sequence models (Pols et al., 1997a) have been unsatisfactory, mainly due to a large and irreconcilable age discrepancy between the members. Fitting both components to pre-MS models, however, produces excellent agreement.

To test the case for PV Cas being pre-MS, we looked for other observation clues. The double-lined eclipsing binary system RS Cha was recently found by Mamajek, Lawson & Feigelson (1999, 2000) to be pre-MS. Not only were pre-MS tracks for RS Cha a better fit than post-MS tracks, but two other observations strengthened the argument: (1) RS Cha had several nearby ROSAT All-Sky Survey X-ray sources nearby which were found to be very young, low-mass, weak-lined T Tauri stars, and (2) RS Cha’s proper motion matched that of the T Tauri stars, suggesting a genetic tie. PV Cas is at a distance of 660 pc (Popper, 1987), and a young stellar aggregate or membership within an OB association could have been previously overlooked.

Searching the Hipparcos and Tycho-2 catalogs, as well as examining PV Cas on the Digitized Sky Survey, we found no evidence for PV Cas being a member of a known OB Association. More massive members of a putative association would be included in the Hipparcos catalog with proper motions similar to that given in the Tycho-2 entry for PV Cas, but none were found. We searched for known groups of young stars with Vizier at CDS: the compilations of OB Associations by Ruprecht, Balazs, & White (1982), Melnick & Efremov (1995), and de Zeeuw et al. (1999), and open clusters by Ruprecht, Balazs, & White (1983) and Lyngå (1987). The only possible known associations that PV Cas could belong to are Cep OB3 (d=840 pc, $`\mathrm{\Delta }\theta =4.0\mathrm{°},v_R=23\mathrm{km}/\mathrm{s}`$) and Cep OB2 (d=615 pc, $`\mathrm{\Delta }\theta =9.5\mathrm{°},v_R=21\mathrm{km}/\mathrm{s}`$), however their projected separations from PV Cas are large (600 pc and 100 pc, respectively), and their average radial velocities are far from Popper’s value for PV Cas ($`v_R=3\mathrm{km}/\mathrm{s}`$). Hence, PV Cas does not appear to be connected to any known OB Associations or clusters which help us to infer its nature.

The ROSAT All-Sky Survey (RASS) Bright Source Catalog (BSC; Voges et al. (1999)) and Faint Source Catalog (FSC; Voges et al. (2000)) were searched to see whether there was any evidence for a clustering of X-ray-emitting T Tauri stars in the vicinity of PV Cas. No concentration of sources near PV Cas was detected, although the sensitivity of RASS at 660 pc is about L<sub>x</sub> $``$ 10<sup>31</sup> erg s<sup>-1</sup>, corresponding to the very high end of the X-ray luminosity function for T Tauri stars (Feigelson & Montmerle, 1999). Only one RASS-BSC source was within 30 ′ ($``$6 pc projected) of PV Cas, but its f<sub>X</sub>/f<sub>V</sub> ratio was 2 magnitudes too high to be a plausible T Tauri star candidate. The only RASS-FSC source within 30 ′ of PV Cas appeared to be related to a galaxy cluster on the Digitized Sky Survey.

We conclude that we currently have no evidence for a pre-MS aggregate around PV Cas which could strengthen the argument for its pre-MS status. However, the Taurus clouds also are forming low mass pre-MS stars without high mass cluster counterparts.

### 2.5 Lower Mass Binaries

Figure 7 shows EK Cep, MY Cyg, PV Pup, DM Vir, V1143 Cyg, and UX Men, whose masses range from 2.03 to 1.12 $`\mathrm{M}_{}`$. The Ribas et al. (2000) effective temperatures improve the fit for UX Men.

EK Cep has a large mass ratio. The optimum fit occurs as EK Cep B is still on the pre-MS track, in agreement with Martin & Rebolo (1993) and Claret, Giménez, & Martin (1995). We find that the surface abundance of $`\mathrm{Li}^6`$ is depleted to about $`10^4`$ of its initial value, while $`\mathrm{Li}^7`$ is depleted from $`1.47\times 10^9`$ to $`0.393\times 10^9`$. This corresponds to a depletion of elemental lithium of about 0.57 dex (base 10 logarithm). This is somewhat larger than found by Martin & Rebolo (1993) (0.1 dex), but may be due to differences in the nuclear reaction rates used. In this range, the depletions are almost linear in the net cross section for $`\mathrm{Li}^7`$ destruction. A careful analysis with a variety of rates is warranted: Martin & Rebolo (1993) suggest that the observations are in conflict with pre-MS models giving a Li depletion greater than 0.3 dex.

Although EK Cep has a large $`\chi ^2`$ ($`\chi ^2=10.3`$) if the radii are used in determining the fitting function, the situation is different for $`\mathrm{log}L`$-$`\mathrm{log}T_{eff}`$, the conventional HR plane. The observational errors are now larger, and the corresponding $`\chi ^2`$ approaches 4. This confirms the importance of using the radii directly as a discriminant(Andersen, 1991).

MY Cyg and UX Men are found to be well into main sequence hydrogen burning. MY Cyg is underluminous relative to the models. A higher heavy element abundance would remove the discrepancy; observational tests of this are needed. Pols et al. (1997b) found $`Z=0.024`$ and Ribas et al. (2000) found $`Z=0.039`$, which are consistent with this suggestion.

PV Pup and V1143 Cyg are on the pre-MS/MS boundary. The fitting procedure chooses the cusp at which the star settles down to main sequence burning. This cusp shifts with small changes in abundance, so these fits would benefit from independent measurement of the abundances in these binaries.

DM Vir has been updated for Latham et al. (1996). Although the changes were small, the new fit is in the middle of main sequence hydrogen burning instead of pre-MS contraction. The track lies well within the error bars; the previous data also had $`\chi ^2<4`$, although a much younger age estimate. Because the masses are almost equal, the coeval birth requirement has little effect, and the ages have a corresponding uncertainty.

### 2.6 Roche Lobes

Observational selection favors binaries with a small separation. In order to determine the true usefulness of these systems as tests of models of single star evolution it is necessary to know to what extent these systems are detached (noninteracting), and how far into the past and the future this condition is satisfied.

In order to answer this question to first order for the systems in our sample, the average Roche lobe radius for each star was calculated using

$$R_{Roche}=\frac{0.49a}{0.6+q^{2/3}\mathrm{ln}(1+q^{1/3})}$$

where $`a`$ is the binary separation and $`q`$ is the mass ratio with the star in question in the numerator (Lewin, van Paradijs, & van den Heuvel, 1995). This average radius was then compared to the model radii to estimate when each star overflows its Roche lobe. Dynamical evolution of the orbits was not taken into account. None of the models indicated significant mass transfer prior to the ages of the models closest to the observed points since early in the pre-Main Sequence evolution.

Four of the binaries in the sample have at least one member which, according to the model radii, will overflow their Roche lobes when they are between 1.3 and 2 times older than their current age. The results of the Roche lobe comparisons are given in Table 4. All stars labeled “Post” do not exceed their Roche Lobe radius until well into their post-main sequence evolution. Two stars are in contact early in the pre-main sequence evolution (EK Cep B and MY Cyg A). The times given for these stars correspond to when they contract below the critical radius and mass transfer ends. These numbers should be taken as a rough guideline at best, since the dynamical evolution of protostars is undoubtedly more complex than the simple algorithm used here. The models corresponding to the primaries EM Car and V478 Cyg exceed their Roche lobe radii in less than $`3\times 10^6`$ years.

These values are approximate in that dynamical evolution is not taken into account, the model radii do not match exactly the observed radii, and an approximate Roche lobe geometry was used to facilitate comparison to the spherically symmetric models.

## 3 APSIDAL MOTION

Apsidal motion in binaries allows us to infer constraints on the internal mass distributions (Schwarzschild, 1957). Apsidal motion, that is, rotation of the orientation of the orbital ellipse relative to an inertial frame, does not occur for binary orbits of point particles interacting by Newtonian gravity. Levi-Civita (1937) showed that the general relativistic expression for the periastron shift of a double star is the same as for the perihelion shift of Mercury. Following Weinberg (1972) (see pages 194-7), the shift is

$$\left(P/U\right)_{GR}=3G(M_A+M_B)P/a(1e^2)c^2,$$
(2)
where $`c`$ is the speed of light and $`G`$ the gravitational constant. Using units of solar masses and radii, and with the period $`P`$ in days, this dimensionless number becomes

$$\left(P/U\right)_{GR}=6.36\times 10^6(M_A+M_B)P/a(1e^2),$$
(3)
apsidal orbits per orbit. Tests of general relativity have reached high precision (Will, 1998); the perihelion shift has now been tested to about $`3\times 10^3`$. There has been some controversy as to a possible breakdown of general relativity because of a discrepancy between observations and predictions of the apsidal motion of some systems. This has been clarified by Claret (see Claret (1997, 1998) for a recent discussion), who pointed out errors in theoretical models and difficulties in observations, especially for systems whose apsidal periods are too long for much to be measured with modern equipment. We adopt the point of view that general relativity is better tested than subtleties in the evolution of binary stars, and ascribe errors to other causes (tidal effects not included, rotational effects, and systematic errors in observational interpretation, for example).

Tides induced by each companion give an additional interaction which is not purely inverse square in the separation and cause apsidal motion. Quataert, Kumar & On (1996) have discussed the validity of the classical formula, which we use,

$$\left(P/U\right)_{CL}=(15/a^5)[k_1R_1^5M_2/M_1+k_2R_2^5M_1/M_2]f(e),$$
(4)
where $`P`$ is the period of the orbit, $`U`$ the period of apsidal motion, $`M_i`$ the mass and $`R_i`$ the radius of the star $`i`$, and

$$f(e)=(1+\frac{3}{2}e^2+\frac{1}{8}e^4)/(1e^2)^5,$$
(5)
where $`e`$ is the eccentricity of the orbit. The separation of the pair in solar radii is

$$a=4.207P^{\frac{2}{3}}(M_1+M_2)^{\frac{1}{3}},$$
(6)
if the period $`P`$ is measured in days and the masses in solar units. The classical apsidal motion formula gives accurate results when the periods of the low-order quadrupole $`g`$, $`f`$ and $`p`$modes are smaller than the periastron passage time by a factor of about 10 or more (Quataert, Kumar & On, 1996). For EM Car, the lowest order pulsational mode of the primary has a period of $`0.324`$ days compared with the orbital period of 3.414 days and an eccentricity of $`0.0120\pm 5`$, so that this condition is just satisfied.

If we assume that the observed apsidal motion is due only to these two effects, classical simple tides and general relativity, we have

$$\left(P/U\right)_{OBS}\left(P/U\right)_{GR}=\left(P/U\right)_{CL}.$$
(7)
We use the products $`k_iR_i^5`$ directly for greater precision, but quote the apsidal constants $`k_i`$ for comparison. Petrova (1995) has pointed out that accuracy problems may exist because the relevant parameter is $`k_iR^5`$, where $`k_i`$ is the apsidal constant and $`R`$ is the stellar radius, not just $`k_i`$ alone.

Figure 8 shows the integrand of the apsidal constant, which approaches an asymptotic value as the integration exceeds about 0.7 of the radius. Inner regions contribute little because of their small radii; outer regions have little mass. The change from the interior (Henyey) integration to envelope integration occurs around $`r/R=0.5`$, and is visible in the change in the density of points. At the join, the temperature is about $`T6\times 10^6\mathrm{K}`$, and the density $`\rho 2.0\times 10^2\mathrm{gcm}^3.`$ This temperature is about ten times the value attained in the early opacity experiments (Perry et al., 1991, 1996) on the NOVA laser, and is about half the goal for the National Ignition Facility (NIF). For such main sequence (and pre-main sequence) stars, the apsidal constants are most sensitive to the range of density and temperature which is directly accessed by high energy density laser experiments (see Remington et al. (1999) and discussion above). In this range, the new opacities show significant deviation from those previously used in astrophysics (Rogers & Iglesias, 1992; Iglesias & Rogers, 1996).

The Petrova & Orlov (1999) catalog contains orbital elements for 128 binaries, including most (11 of 18) of the binaries in our list. Table 5 gives apsidal constants $`k_i`$, as well as the observed and predicted ratios of $`P/U`$. Given the significant improvement in the opacities, a critical re-examination of these data seems warranted.

Figure 9 shows the dimensionless rate of apsidal motion, $`(P/U)_{CL}=(P/U)_{OBS}(P/U)_{GR}`$, which would be due to classical apsidal motion, plotted versus log of half the total binary mass. $`P`$ is the orbital period and $`U`$ the apsidal period. The observational data (corrected for general relativity) are shown as diamonds, with vertical error bars. The model predictions are shown as solid squares ($`\chi ^2<4`$) for the best fits, open triangles for $`4<\chi ^2<16`$, and crosses for $`\chi ^2>16`$. The massive binaries with good fits (EM Car, V478 Cyg, and CW Cep; $`(M_A+M_B)/2>10M_{}`$) have predicted apsidal motion in excess of that observed, and QX Car also follows that trend. These models are not as centrally condensed as the stars. This may be related to the underluminosity of these models found above. Additional mixing would give more massive, convective cores, which would result in both greater luminosity and more centrally condensed structure.

Of the lower mass binaries with measured apsidal motion, only PV Cas has a good fit model. Its predicted apsidal motion is also larger than that observed (the stars are more centrally condensed). The other binaries need better fitting models before the tests can be convincing. Note that at the lowest tick mark in Figure 9, the apsidal period is measured in centuries.

Claret & Giménez (1993) have shown that inclusion of (a) overshooting, (b) variation in heavy element abundance, and (c) rotation can produce models consistent with the apsidal data; see also (Claret, 1999). This represents good progress toward establishing the apsidal motion data as a useful test of stellar evolution. Our results, while not based upon identical assumptions, are consistent. The challenge is that of correctly determining the relative importance of the several different small effects which can give consistency with the observations.

## 4 CONCLUSIONS

Standard stellar evolution, without embellishments such as overshooting and rotation, does fairly well on these quantitative tests. While clear discrepancies exist, they are relatively subtle. This makes it difficult to uniquely identify exactly which additional physics is needed. We find a detailed consistency with similar calculations by Pols et al. (1997b) and Ribas et al. (2000). It is important to test observationally the abundance variations implied by $`\chi ^2`$ optimization, as such procedures may hide missing physics in parameter variation. Laser experiments now explore the regions of temperature and density which are relevant not only to conventional stellar evolution, but also to apsidal motion tests.

Massive stars require more mixing than given by standard stellar evolution, and probably more than the prescriptions for overshooting used by Pols et al. (1997b) and Ribas et al. (2000). Rotational mixing, overshooting which is mass dependent, or something else is needed. Our mass loss prescription is near the observational upper limit, so that additional mass loss is an unlikely solution.

Lower mass stars with convective cores are not well fit by standard stellar evolution. Again, additional mixing is a promising answer. Several of these binaries seem to be pre-main sequence; this will allow some interesting tests of depletion of light nuclei and the mixing processes.

We find two serious challenges: (1) disentangling conflicting solutions of the relatively subtle discrepancies, and (2) controlling shifts in the observational “target areas” due to systematic errors, which seem to be larger than the statistical errors. Improved determinations of effective temperature, and of heavy element abundances (e.g., $`[\mathrm{Fe}/\mathrm{H}]`$), would greatly improve these tests.

Stimulating discussions with Brigitta Nordstrom and Johannes Andersen are gratefully acknowledged. Special thanks are due to Marcel Arnould, Shimon Asida, John Bahcall, Audra Baleisis, Al Cameron, Joergen Christensen-Dalgaard, Paul Drake, Rolf Kudritzki, Robert Kurucz, T. J. Pearson, Ted Perry, Thomas Rauscher, Marc Rayet, Bruce Remington, Friedel Thielemann, Frank Timmes, and Cliff Will, whose work and/or advice contributed to aspects of this paper. This work was supported in part by DOE, grant number DE-FG03-98DP00214/A001, and a subcontract from ASCI Flash Center at U. of Chicago.

EM Car (catalog ) V478 Cyg (catalog ) CW Cep (catalog ) QX Car (catalog ) CV Vel (catalog ) U Oph (catalog )

zeta Phe (catalog ) IQ Per (catalog ) PV Cas (catalog ) AI Hya (catalog ) VV Pyx (catalog ) RS Cha (catalog )

EK Cep (catalog ) MY Cyg (catalog ) PV Pup (catalog ) DM Vir (catalog ) V1143 Cyg (catalog ) UX Men (catalog )


format: text

text:
# Strong ferromagnetism and weak antiferromagnetism in double perovskites: Sr2FeMO6 (M=Mo, W, and Re)

## Abstract

Double perovskites Sr2Fe$`M`$O6 ($`M`$=Mo and Re) exhibit significant colossal magnetoresistance even at room temperature due to the high Curie Temperature (419K and 401K). However, such a high Curie Temperature is puzzling, given the large separation between magnetic elements (Fe). Moreover, with $`M`$=W, the electronic and magnetic properties suddenly change to insulating and antiferromagnetic with the Néel temperature of only 16$``$37 K. Based on detailed electronic structure calculations, a new mechanism is proposed which stabilizes the strong ferromagnetic state for $`M`$=Mo and Re and is passivated for $`M`$=W.
PACS number: 75.30.Et, 75.30.Vn, 71.20.Be

Intensive studies on the perovskite transition-metal oxides (TMO), particularly manganites, have revealed a variety of novel phenomena only for half a decade . Among those phenomena, colossal magnetoresistance (CMR) has been attracting strong attention not only as a challenging subject of fundamental science but also as an important phenomenon for potential technological application. With regard to the latter aspect, materials with not only the half-metallic nature but also Curie temperature ($`T_c`$) much higher than room temperature is strongly desired in order to realize strong CMR effects at room temperature. It was demonstrated that some of the double perovskite TMO such as Sr2FeMoO6 (SFMO) and Sr2FeReO6 (SFRO) are suitable candidates . They are half metallic according to the band structure calculations and their $`T_c`$’s are 419K and 401K.

The present work deals with two fundamental problems in these double perovskite TMO by performing detailed electronic structure calculations. In both of SFMO and SFRO, the magnetic moments of Fe are aligned ferromagnetically and the induced moments on Mo and Re are coupled antiferromagnetically to Fe moments. Therefore these materials can be regarded as ferrimagnetic. However, we regard them ferromagnetic (FM) because Mo and Re are intrinsically non-magnetic in the sense that their magnetic polarization cannot be sustained spontaneously by the exchange potential on these atoms. Actually their negative moments (i.e., antiparallel to Fe moments) are induced by Fe moments through the 4d(5d)-3d hybridization. Now the first fundamental question is why $`T_c`$ is so high despite the fact that Fe atoms are very much separated with non-magnetic elements (Mo, Re) sitting in between. We will point out that a FM stabilization mechanism proposed by Kanamori and Terakura  operates in SFMO and SFRO. The same problem was treated recently also by Sarma et al. for SFMO . The second question concerns the striking difference of Sr2FeWO6 (SFWO) from SFMO and SFRO in the electronic and magnetic properties. SFWO is antiferromagnetic (AF) insulator and the Néel temperature is only 16$``$37K . Why is the W case so different from the Mo and Re cases despite the fact that W is the 5d analogue of Mo and next to Re in the row of the periodic table?

We will show that the stronger 2p(O)-5d(W) hybridization compared with 2p(O)-4d(Mo) hybridization is the main source of the difference between Mo and W pushing the 5d states higher in energy and passivating the FM stabilization mechanism in SFWO. Deeper 5d levels in the Re case compared with the W case cancel the effect of enhanced p-d hybridization and restores the FM stabilization mechanism. Note, however, that the standard LSDA (local spin-density approximation) or GGA (generalized gradient approximation)  cannot describe properly the ground state of SFWO. As Fe d states are strongly localized in these systems, the local Coulomb repulsion $`U_{\mathrm{eff}}`$, which is semi-empirically taken into account by the LDA+U method  in the present work, plays crucially important roles.

We adopt the plane-wave pseudopotential method. The 3d states of Fe, 4d states of Mo, 5d states of W and Re and 2p states of O are treated with the ultrasoft pseudopotentials  and the other states by the optimized norm-conserving pseudopotentials . The cut-off energy for describing the wave functions is 30 Ry, while that for the augmentation charge is 200 Ry. The crystal structures of all the three systems are cubic with the nearest Fe-M distance given as 3.945 Å, 3.975 Å and 3.945 Å for $`M`$=Mo, W and Re, respectively . For the FM state, the number of k-points used in the k-space integration is 19 in the irreducible Brillouin zone. For the AF state, two different configurations, AFI and AFII , are considered. In the AFI (AFII) configuration, the magnetic moments are aligned ferromagnetically within the (001) ((111)) plane and alternate along the  ()direction. The number of k-points in the AF configuration is chosen to be equivalent to that in the FM configuration. As for the electron-electron interaction, we adopt first the standard GGA  and then the semi-empirical LDA+U method . The details of the implementation of the LDA+U method in the pseudopotential scheme can be found in our previous publication .

Figure 1 shows a summary of the GGA calculations for three materials Sr2Fe$`M`$O6 with $`M`$=Mo, W and Re in both FM and AF states. (Note that only the results for AFII are shown here for the AF states.) The thin solid lines denote the local density of states (LDOS) for Fe 3d orbits and thick broken lines LDOS for 4d (Mo) or 5d (W, Re) states. The results for $`M`$=Mo and Re in FM state are basically the same as those shown in the previous works . The oxygen p bands extends from -8 eV to about -4 eV, the Fe majority spin t2g bands from about -4 eV to -2 eV followed by the majority spin eg bands extending up to near the Fermi level. In the majority spin state, the band just at and above the Fermi level is of t2g character of $`M`$. In the minority spin state, t2g states of Fe and $`M`$ coexist around the Fermi level. The formal valence of the combination of Fe$`M`$ is +8, meaning that the number of d electrons per Fe$`M`$ is 6 for $`M`$=Mo and W and 7 for $`M`$=Re. In both FM and AF states, the majority spin bands of Fe are completely filled with 5 electrons and the minority spin bands accommodate 1 electron for $`M`$=Mo and W and 2 electrons for $`M`$=Re. These materials are predicted to be metallic in both FM and AF orders with GGA calculations. Particularly, they are half metallic in the FM order and this half metallicity is preserved even in the LDA+U calculation as shown later. In Table I, the total energies for AF states with reference to those for FM states are given for the three materials. The calculation for the AFI state for $`M`$=W suggests that this magnetic order may have no chance of being realized in these materials. Therefore the AFI order will not be considered hereafter. Clearly, the FM state is significantly more stable than the AFII state in the GGA calculations for all the three systems. The results in GGA are qualitatively consistent with experimental facts for SFMO and SFRO but inconsistent for SFWO, which is antiferromagnetic and insulating experimentally . Deferring the discussion on the stability of ferromagnetism for SFMO and SFRO for a while, we first discuss the problems of SFWO and how to solve them.

The LDOS for the AF state of SFWO has a very sharp peak of the Fe t2g state origin just at the Fermi level. This suggests that the AFII state obtained in this stage may be unstable. Although the symmetry in the AFII state is reduced to D3d, lift of degeneracy in t2g orbits is not strong enough to split the t2g band. The situation is quite similar to FeO . In this case, the lattice is elongated along $`<111>`$ direction (even with the GGA level treatment) and furthermore the local Coulomb repulsion ($`U_{\mathrm{eff}}`$) strongly enhances the orbital polarization making the system insulating. On the analogy of FeO, we first studied effects of rhombohedral distortion of SFWO in GGA and found that such distortion either elongation or contraction along $`<111>`$ direction simply increases the total energy. The cubic lattice for SFWO even in the AFII state is actually observed experimentally. As these analyses suggest that there is little chance of stabilizing the AF state for SFWO with the GGA level calculation, we applied the LDA+U method to these materials. Although the LDA+U method is semi-empirical, it still provides us with some important insights into the problems. We set $`U_{\mathrm{eff}}`$ to be 4 eV and applied it only to the Fe d orbitals for the sake of simplicity. As was described in our previous paper, $`U_{\mathrm{eff}}`$ is nonzero in a rather limited region around the Fe nucleus and its actual value does not have definite meaning .

Figure 2 shows the LDA+U version of Fig.1, and Table I includes the corresponding total energies. The common characteristic feature in Fig.2 is the enhancement in the exchange splitting of Fe. Nevertheless, the electronic structure remains qualitatively the same for SFMO and SFRO in both FM and AF states except the fact that the weight of Mo and Re d states increased significantly around the Fermi level. On the other hand in SFWO, while the change in the FM state is minor, the AF state shows a dramatic change from Fig.1 to Fig.2. The t2g band of Fe splits due to orbital polarization induced by $`U_{\mathrm{eff}}`$ and the occupied state in the minority spin state just below the Fermi level is of a1g character. A band gap opens up and the AFII state becomes more stable than the FM state. The insulating nature of the ground state of SFWO is now correctly reproduced . The fact that the Néel temperature is only 16$``$37 K may suggest that the stabilization of the AFII state in the present calculation may be overestimated. However, the quantitative aspect can be tuned by $`U_{\mathrm{eff}}`$.

Having shown the calculated results which are qualitatively consistent with experimental facts, we start discussions on the underlying mechanisms in relation to the two fundamental questions raised at the beginning of the present paper. The first one concerns the mechanism of the strong stabilization of the FM state for SFMO and SFRO. Recently Sarma et al. proposed an interesting explanation to the origin of strong AF coupling between Fe and Mo in which they pointed out strong effective exchange enhancement at Mo due to the 3d (Fe)-4d(Mo) hybridization. Kanamori and Terakura  proposed a more general idea for the mechanism where a non-magnetic typical element located at the midpoint of neighboring high-spin 3d elements contributes to stabilization of the FM coupling of the 3d elements. Figure 3(a) is a schematic illustration explaining the mechanism. The states of the typical element located in between the majority and minority spin states of 3d elements are tentatively called p states. The key concept in this mechanism is the energy gain contributed by the negative spin polarization of the non-magnetic element induced by the p-d hybridization. Such spin polarization does not exist in the AF configuration and therefore there is no energy gain due to the spin-state relaxation at the typical element. In the present problem, the 4d states of Mo and 5d states of W and Re correspond to the p states in Fig.3(a). The analogy is obvious in the majority spin state in the FM order. In the minority spin state, as the 4d (or 5d) bands and the 3d bands are not well separated and the Fermi level lies in the 3d bands after including the hybridization, we need a careful analysis to distinguish the FM stabilization mechanism discussed above and the double exchange (DE) . Figure 3(b) illustrates the situation corresponding to SFMO and SFRO where the $`M`$ t2g bands are slightly below the Fe ones. We first treat the up and down spin states separately and then consider the electron transfer between two spin states. The standard DE mechanism takes account of the processes only in the first step. As for the hybridization between $`M`$ (=Mo, W, Re) bands and the majority spin Fe bands, the total energy change caused by band shift due to the 3d-4d(5d) hybridization does not depend on the relative spin direction between the neighboring Fe atoms up to the second order in the hybridization matrix element $`t`$. For example, the upward shift of the up spin $`M`$ bands by $`2t^2/\mathrm{\Delta }`$ in the FM state balances the upward shift of the both spin $`M`$ bands by $`t^2/\mathrm{\Delta }`$ in the AF state where $`\mathrm{\Delta }`$ denotes the energy separation. Subtle features exist in the minority spin state. Not only the band shift but also band broadening have to be considered. It is obvious that the width of the minority spin bands will be wider in the FM state than in the AF state. As the Fermi level lies in the minority spin bands, the band broadening contributes to the stability of the FM state like in the standard DE . In the present problem, we have an additional effect in the FM state coming from the electron transfer just like in Fig. 3(a). This electron transfer produces negative spin polarization at $`M`$ atoms and contributes to further stabilization of the FM state. In contrast to SFMO and SFRO, the $`M`$ t2g bands in SFWO are slightly above the Fe ones. In this case, the $`M`$ t2g bands are basically empty and the electron transfer will not occur. Therefore the FM-stabilization mechanism of Fig. 3(a) is passivated for SFWO, while the DE mechanism may still be effective. Although Table I still shows small negative spin polarization at the W atom, this is due to the stronger 5d-3d hybridization in the minority spin state than in the majority spin state. We also speculate that the main reason of considerable relative stability of the FM order for SFWO in the GGA calculation is the rather unstable electronic configuration in the AFII order. Because of this, a change in the electronic structure in the AFII state from GGA to LDA+U reduces the energy of the AFII state dramatically.

The second question concerns the origin of the different behavior of W from other two elements Mo and Re. It is clear from the above arguments that in order to answer this question, we have to clarify the origin of the difference in the energy position of the minority spin t2g bands. We assign the p-d hybridization between oxygen and $`M`$ to the main source of this difference. As the 5d orbital of W is more extended than the 4d orbital of Mo, the stronger 2p(O)-5d(W) hybridization pushes the 5d band, which is the p-d antibonding state, higher in energy. This mechanism is supported by the fact that the p-d bonding counter part is clearly deeper for SFWO than for SFMO (see Figs.1 and 2). As Re has deeper 5d level than W to accommodate one more d electron, the energy scheme for SFRO becomes similar to that for SFMO.

In summary, we showed that the electronic structures and magnetic ordering in the ground state of Sr2Fe$`M`$O6 ($`M`$=Mo, W and Re) are properly reproduced by the LDA+U method. A new mechanism was proposed to explain the high Curie temperature for $`M`$=Mo and Re cases. An explanation was also given to the sudden changes in the electronic and magnetic properties in the $`M`$=W case. The mechanism proposed by us  is very useful to predict qualitatively the change in the magnetic states by changing constituent elements.

We thank Tokura Group members in JRCAT for providing us with experimental information. Thanks are given also to Prof. N. Hamada and Prof. D. D. Sarma for valuable discussion. The present work was partly supported by NEDO.


text:
# The radio surroundings of the microquasar GRO J1655-40

## 1 Introduction

The soft X-ray transient GRO J1655$``$40 is one of the best black hole candidates in the Galaxy (Mirabel & Rodriguez 1999). It is a low mass X-ray binary whose compact object has a mass in the range $`5.57.9`$ $`M_{}`$ (Shahbaz et al. 1999). The system is located at 3.2 kpc from the Earth and presents frequent superluminal ejections at radio wavelengths (e.g. Tingay et al. 1995, Hjellming & Rupen 1995).

Recently, Israelian et al. (1999) found convincing evidence of a supernova origin for the black hole in GRO J1655$``$40. They detected a clear overabundance of $`\alpha `$-elements in the secondary star, whose internal temperatures are not high enough as to synthesize them. These elements seem to have been produced in a massive companion (with $`2540`$ $`M_{}`$) which exploded both contaminating the secondary star with matter enriched with nucleosynthetic products and creating the black hole. Israelian et al. (1999) estimated that the explosion should have occurred less than $`10^6`$ years ago.

In this Letter, we present the results of a study of the radio environment of GRO J1655$``$40. Our aim is to detect any existent signature of an explosive event in this region. A supernova explosion should have introduced modifications in the ISM around the binary system, displacing material and creating a low density bubble. Additionally, the expanding shock front should accelerate electrons up to relativistic energies, producing a shell-type source in the radio continuum. If the age of the system is less than $`10^6`$ years, these observational features might be detectable.

## 2 Data analysis and new observations

### 2.1 Radio continuum

We have used data from the 4.85-GHz PMN survey (Condon et al. 1993) in order to produce a large-scale map ($`1^{}\times 1^{}`$) of the surroundings of GRO J1655$``$40. Although these data are optimized for sizescales of less than 40 arcminutes, the good resolution and sensitivity of the images render them a valuable tool to search for SNR candidates, especially if they are not of large angular size (Duncan et al. 1997). For this work, we apply to the 4.85-GHz data an additional filtering process in order to remove the galactic diffuse emission on scales larger than 30 arcminutes (see Combi et al. 1998 for details of the Gaussian filtering method). The resulting map is shown in Figure 1, along with a 1.4-GHz image of the small-scale emission of the inner region obtained from the NVSS Sky Survey (Condon et al. 1998) at a better resolution. These latter observations were obtained with the Very Large Array (VLA) in the compact D and DnC configurations.

### 2.2 HI observations

We have performed HI observations towards the area of interest with a 30-m single dish telescope located at the Instituto Argentino de Radioastronomía (IAR), Villa Elisa, Argentina. The observations were carried out during four consecutive sessions on November 20-23, 2000. The receiver is a helium-cooled HEMT amplifier with a 1008-channel autocorrelator at the backend. System parameters and additional details of the observational technique can be found in Combi et al. (1998). The HI line was observed in hybrid total power mode and the sky was sampled on a $`0.35^{}`$ rectangular grid. Each grid position was observed during 60 s with a velocity resolution of $`1`$ km s-1 and a coverage of $`\pm 450`$ km s-1. A set of HI brightness temperature maps ($`\mathrm{\Delta }T_{\text{rms}}0.2`$ K) were made for the velocity interval $`\mathrm{\Delta }v=`$($`80`$ km s-1, $`+20`$ km s-1). Those maps for the interval $`42`$ to $`30`$ km s-1, corresponding to the HI distribution at a distance centered at $`3.5`$ kpc, are shown in Figure 2.

## 3 Main results

Figure 1 (upper panel) shows the radio continuum image at 4.85 GHz of the surroundings of GRO J1655$``$40. A small ($`24^{}\times 24^{}`$) semi-circular radio source, which very much resembles a SNR, can be clearly seen in the center of the image. This extended, partial shell source is centered at $`(l,b)(345.14^{},+2.48^{})`$. It has an integrated flux density of $`0.5`$ Jy at 4.85 GHz, although uncertainties due to the background filtering make difficult to estimate a realistic value. The microquasar is located at a distance of $`0.15`$ deg from the center of the radio source.

Figure 1 (lower panel) shows the VLA image at 1.4 GHz (Condon et al. 1998) of the vicinity of GRO J1655$``$40. The microquasar was undetectable at the epoch of the observations and its position is marked with a cross. There are seven point sources near the microquasar with fluxes above 5 mJy. They have been labeled from S1 to S7 in the image (notice that the SNR candidate is not visible in this map because these data are optimized for point source detection). The bright, compact source located at $`(l,b)345.1^{}`$, $`+2.37^{}`$ is identified in the PMN catalogue as PMN J1654-3949. The sources S1 to S6 have no entry in any point source catalog at present. The measured characteristics of all these sources are summarized in Table 1.

Figure 2 clearly shows the development of a cavity around GRO J1655$``$40 in the HI maps from $`42`$ to $`30`$ km s-1. Although the cavity do not form a shell due to the steep density gradient towards the galactic plane, the existence of a local minimum in the HI distribution at the position of the microquasar is evident. This kind of “horse-shoe” morphology has been observed around early-type stars with very strong stellar winds (e.g. Benaglia & Cappa 1999), where a strong energy release in a highly inhomogeneous medium has occurred.

In Figure 3 we present the integrated column density map for the velocity interval from $`40`$ to $`30`$ km s-1. Standard galactic rotation models (Fich et al. 1989) indicate that this HI hole is located at a kinematic distance of 3 to 4 kpc, in good accordance with the estimated distance to GRO J1655$``$40.

## 4 Discussion

If we assume that the extended radio structure shown in Fig. 1, upper panel, is a supernova remnant located at the same distance than the microquasar, we get that its radius should be $`15`$ pc. The surface brightness of the source results to be abnormally low for such a remnant ($`\mathrm{\Sigma }_{4.8\mathrm{GHz}}1.3\times 10^{22}`$ W m-2 Hz-1 sr-1). On the other hand, the microquasar is displaced about $`0.15^{}`$ from the center of the semi-shell source. Although the transverse velocity of GRO J1655$``$40 is not known, it is constrained to be less than 300 km s-1 (Brandt et al. 1995). At this velocity, it would have taken more than $`2.7\times 10^4`$ yr to reach its current position. And for a velocity similar to the observed radial velocity ($`114`$ km s-1), the age of the system would be $`\mathrm{82\hspace{0.17em}500}`$ yr. These timescales are incompatible with the size of the remnant and our knowledge of the ambient density in this region (see below), and lead us, consequently, to discard the extended radio source in Fig. 1 as possibly associated with GRO J1655$``$40. This source could be a background SNR or, perhaps, an artifact formed by several close and unresolved point sources.

At larger sizescales there are not other clear continuum structures that could be associated with a SNR. We have inspected larger maps, of $`6^{}\times 6^{}`$, in search of additional signatures, but nothing was found. This is a very confused region located nearly towards the Galactic Center and the detection of weak continuum features with low surface brightness is extremely difficult. Even at small scales there is a significant number of sources, as it is shown in Fig. 1, lower panel. Here, the two closest sources to the microquasar, S3 and S6, seem to be at approximate equal angular distances, in opposite directions, of the X-ray binary. However, it is unlikely that they could be old emission knots ejected by the central object, because they are on an axis displaced $`30^{}`$ from the axis determined by the known ejections (Hellming & Rupen 1995), which is also shown in the figure.

Contrary to the often misleading continuum data, the HI observations can provide more clear signatures of an explosive event in this region, because the differential rotation motion of the material allows a glimpse of the matter distribution along the line of sight. The channel maps in Fig. 2 show the development in the velocity distribution of a cavity in the HI with a “horse-shoe” shape and an average radius of $`75`$ pc at a distance of 3.2 kpc. The integrated column densities shown in Fig. 3 can be used to estimate the mass displaced to form the cavity. Assuming that the background density is given by the density at the minimum, we estimate that at least 25 000 $`M_{}`$ have been removed from this region. Since the average expansion velocity of the shell is $`15`$ km s-1, we obtain that a lower limit to the original energy release which started the mass motion at the center of the cavity is $`6\times 10^{49}`$ erg. The expansion velocity should have been considerably higher in the past, in such a way that the real value is probably in the range $`10^{50}10^{51}`$ erg, what is quite consistent with the supernova hypothesis. This is reinforced by the fact that there are not known early-type stars with stellar winds strong enough to create such a bubble in this sector of the sky. It is noteworthy that even higher values for the initial energy release are possible if the explosion was caused by a hypernovae event as recently proposed by Brown et al. 2000.

The mean number density in this region can be obtained averaging the column densities within the bubble. This provides a rough estimate of the average density of the ISM before the explosion. We obtain a value $`n_00.4`$ cm-3. This density and the size of the cavity ($`R75`$ pc) indicate that the SNR should be well within the radiative phase. An exact estimate of the time elapsed since the explosion is not possible without a knowledge of the detailed structure of the ISM in this region, but simulations with an original energy release of few times $`10^{50}`$ erg and some reasonable conditions for the cloud properties suggest ages in the range $`10^510^{5.5}`$ yr (Cowie et al. 1981). The continuum emission of the remnant, with sizescales of several degrees, should have been removed during the filtering of the background radiation in the PMN survey. New observations with good sensitivity and single dish telescopes over a large area are necessary to detect the SNR radio structure and help to make a more accurate estimate of its age, if the remnant is not already dissipated.

In summary, we have presented radio continuum and HI maps of the surroundings of the black hole candidate GRO J1655$``$40. A bubble centered at the microquasar has been detected in our HI data. The energy necessary to create this bubble and its approximate age are consistent with the hypothesis that it is the trace left in the ISM by the supernova explosion occurred when the black hole was formed.

###### Acknowledgements.

This work was supported by CONICET (PIP 0430/98), ANPCT (PICT 03-04881) and Fundación Antorchas.


text:
# Measurement of Lagrangian velocity in fully developed turbulence

## Abstract

We have developed a new experimental technique to measure the Lagrangian velocity of tracer particles in a turbulent flow, based on ultrasonic Doppler tracking. This method yields a direct access to the velocity of a single particle at a turbulent Reynolds number $`R_\lambda =740`$. Its dynamics is analyzed with two decades of time resolution, below the Lagrangian correlation time. We observe that the Lagrangian velocity spectrum has a Lorentzian form $`E^L(\omega )=u_{\mathrm{rms}}^2T_L/(1+(T_L\omega )^2)`$, in agreement with a Kolmogorov-like scaling in the inertial range. The probability density function (PDF) of the velocity time increments displays a change of shape from quasi-Gaussian a integral time scale to stretched exponential tails at the smallest time increments. This intermittency, when measured from relative scaling exponents of structure functions, is more pronounced than in the Eulerian framework.

Lagrangian characteristics of fluid motion are of fundamental importance in the understanding of transport and mixing. It is a natural approach for reacting flows or pollutant contamination problems to analyze the motion of individual fluid particles Pope . Another characteristic of mixing flows is their high degree of turbulence. For practical reasons, most of the experimental work concerning high Reynolds number flows has been obtained in the Eulerian framework. Lagrangian measurements are challenging because they involve the tracking of particle trajectories: enough time resolution, both at small and large scales, is required to describe the turbulent fluctuations.

Early Lagrangian information have been extracted from the dispersion of particles, following Taylor’s approach. Recently numerical and experimental studies have focused on resolving the motion of individual fluid or tracer particles. The emerging picture is as follows. The one-component velocity auto-correlation function is quasi-exponential with a characteristic time of the order of the energy injection scale VirantDracos97 ; Sato ; YeungPope89 . The velocity power spectrum is expected to have a scaling $`E^L(\omega )\omega ^2`$, as recently reported Lien ; Yeung2001  and expected from a Kolmogorov similarity arguments. In the same spirit, the second order structure function should scale as $`D_2^L(\tau )=C_0ϵ\tau `$, where $`ϵ`$ is the the power dissipation. Measurements of atmospheric balloons Hanna  have given $`C_0=4\pm 2`$ and a limit $`C_07`$ has been suggested in stochastic models Sawford . Recent experiments VothBodenschatz98 ; PortaBodenschatz01  using high speed optical techniques have shown that the statistics of the Lagrangian acceleration are strongly non-Gaussian.

We have developed a new experimental method, based on sonar techniques MordantJASA , in order to study in a laboratory experiment the Lagrangian velocity across the inertial range of time scales. We obtain the first measurement of single particle velocity for times up to the flow large scale turnover time, at high Reynolds number. In this Letter, we report the results of this measurements and compare with previous observations and numerical predictions.

Our technique is based on the principle of a continuous Doppler sonar. A small (2mm$`\times `$2mm) emitter continuously insonifies the flow with a pure sine wave, at frequency $`f_0=2.5`$ MHz (in water). The moving particle backscatters the ultrasound towards an array of receiving transducers, with a Doppler frequency shift related to the velocity of the particle: $`2\pi \mathrm{\Delta }f=𝐪.𝐯`$. The scattering wavevector $`𝐪`$ is equal to the difference between the incident and scattered directions. A numerical demodulation of the time evolution of the Doppler shift gives the component of the particle velocity along the scattering wavevector $`𝐪`$. It is performed using a high resolution parametric method which relies on an Approximated Maximum Likelihood scheme coupled with a generalized Kalman filter MordantJASA . The study reported here is made with a single array of transducers so that only one Lagrangian velocity component is measured.

The turbulent flow is produced in the gap between two counter-rotating discs mor97 . This setup has the advantage to generate a strong turbulence in a compact region of space, with no mean advection. In this way, particles can be tracked during times comparable to the large eddy turnover time. Discs of radius $`R=9.5`$ cm are used to set water into motion inside a cylindrical vessel of height $`H=18`$ cm. To ensure inertial entrainment, the discs are fitted with 8 blades with height $`h_b=5`$ mm. In the measurement reported here, the power input is $`ϵ=25`$ W/kg. It is measured on the experiment cooling system, from the injection-dissipation balance. The integral Reynolds number is $`Re=R^2\mathrm{\Omega }/\nu =\mathrm{6.5\hspace{0.33em}10}^4`$, where $`\mathrm{\Omega }`$ is the rotation frequency of the discs (7.2 Hz), and $`\nu =10^6`$ m2/s is the kinematic viscosity of water. A conventional turbulent Reynolds number can be computed from the measured $`rms`$ amplitude of velocity fluctuations ($`u_{rms}=0.98`$ m/s) and an estimate of the Taylor microscale ($`\lambda =\sqrt{15\nu u_{rms}^2/ϵ}=0.88`$mm); we obtain $`R_\lambda =740`$. This value is consistent with earlier studies in the same geometry; it corresponds to the range of turbulent Reynolds numbers where measurements of particle acceleration have been reported PortaBodenschatz01 .

The flow is seeded with a small number of neutrally buoyant (density 1.06) polystyrene spheres with diameter $`d=250\mu `$m . It is expected that the particles follow the fluid motion up to characteristic times of the order of the turbulence eddy turnover time, at a scale corresponding to their diameter, i.e. $`\tau _{\mathrm{min}}d/u_dϵ^{1/3}d^{2/3}`$, using standard Kolmogorov phenomenology. For beads of diameter 250 $`\mu `$m, one estimates $`\tau _{\mathrm{min}}1.3`$ ms. This value is within the resolution of the demodulation algorithm whose cut-off frequency is at 3 kHz. Note that the Kolmogorov dissipative time ($`\tau _\eta =\sqrt{\nu /ϵ}=0.2`$ ms) is smaller, so that we do not expect to resolve the dissipative region. The statistical quantities are calculated from $`3\times 10^6`$ velocity data points, taken at a sampling frequency equal to 6500 Hz. The acoustic measurement zone is in central region of the flow, 10 cm thick in the axial direction and almost spanning the cylinder cross-section. In this region the flow is a good approximation to isotropic and homogeneous conditions: at all points, the mean velocity is non zero, but equal to about one tenth of its $`rms`$ value.

We first consider the Lagrangian velocity auto-correlation function:

$$R^L(\tau )=\frac{v(t)v(t+\tau )_t}{v^2}.$$
(1)
We observe – Fig.1a – that it has a slow decrease which can be modeled by an exponential function $`\rho _v(\tau )e^{\tau /T_L}`$. This expression defines an integral Lagrangian time scale $`T_L=22`$ ms. For comparison, the period of rotation of the discs is 140 ms and the sweeping period of the blades is 17 ms. The measured Lagrangian time scale thus appears as a time characteristic of the energy injection. The exponential reproduces extremely well the variation of the auto-correlation function, from about $`5\tau _\eta `$ at small scales to $`4T_L`$. These limits coincide with the upper and lower resolution of the technique, so that we observe an exponential decay over the entire range of our measurement. However, as the variance of the acceleration must be finite PortaBodenschatz01  there has to be some lower cut-off to this behavior, at times of order $`\tau _\eta `$. These observations extend and confirm previous numerical and experimental studies at moderate Reynolds numbers VirantDracos97 ; Sato ; Yeung2001 . Note that the exponential decay of the Lagrangian velocity auto-correlation is a key feature of stochastic models of dispersion since it appears as a linear drift term in a Langevin model of particle dynamics Pope ; Sawford .

We show in Fig.1b the velocity power spectrum, computed both from the data and as the Fourier transform of the exponential decay of the auto-correlation function:

$$E_{\mathrm{fit}}^L(\omega )=\frac{u_{rms}^2T_L}{1+(T_L\omega )^2}.$$
(2)
We observe a clear range of power law scaling $`E^L(\omega )\omega ^2`$. This is in agreement with a Kolmogorov K41 picture in which the spectral density at a frequency $`\omega `$ is a dimensional function of $`\omega `$ and $`ϵ`$: $`E^L(\omega )ϵ\omega ^2`$. To our knowledge, this is the first time that it is directly observed at high Reynolds number and in a laboratory experiment, although it has been reported in oceanic studies Lien  and in lower Reynolds number direct numerical simulations Yeung2001 . Departure from the Kolmogorov behavior is observed at low frequencies in agreement with the exponential decay of the auto-correlation. At high frequencies, the spectrum deviates from the Lorentzian form due to the particle response. Note in Fig.1b that the measurement is made over a dynamical range of about 60 dB.

We now consider the second order structure function of the velocity increment

$$D_2^L(\tau )=(v(t+\tau )v(t))^2_t=(\mathrm{\Delta }_\tau v)^2.$$
(3)
We emphasize that these are time increments, and not space increments as in the Eulerian studies. The profile $`D_2^L(\tau )`$ is shown in the inset of Fig.2. It is linked to the auto-correlation by $`D_2^L(\tau )=2u_{\mathrm{rms}}^2\left(1R^L(\tau )\right)`$: at small times one observes the trivial scaling $`D_2^L(\tau )\tau ^2`$ and at large times $`D_2^L(\tau )`$ saturates at $`2u_{\mathrm{rms}}^2`$ (as $`v(t)`$ and $`v(t+\tau )`$ are uncorrelated).

In between these two limits, one expects an inertial range of scales with a Kolmogorov-like scaling

$$D_2^L(\tau )=C_0ϵ\tau ,$$
(4)
where $`C_0`$ is a ‘universal’ constant. Such a behavior is consistent with dimensional analysis and with an $`\omega ^2`$ scaling range in the velocity power spectrum. Fig.2 shows $`D_2^L(\tau )/ϵ\tau `$; a plateau with a constant $`C_0`$ is not observed. Note that this also the case in Eulerian measurements when the third order structure function is represented in linear coordinates mal98 . The function reaches a maximum at $`20\tau _\eta `$, for which $`C_02.9`$. This value is in agreement with the estimation $`C_0=4\pm 2`$ in Hanna  and in the range of values (between 3 and 7) used in stochastic models for particle dispersion Du . In our case there may also be a bias at small times due to particle effects. However if we assume the exponential fit for the velocity autocorrelation function to be valid down to the smallest scales, we obtain a value $`C_0=3.5`$ as an upper bound for the maximum of $`D_2^L(\tau )/ϵ\tau `$. In our set of measurements between $`R_\lambda =100`$ and $`R_\lambda =1100`$, we have observed an increase of $`C_0`$ (defined in the same way) from 0.5 to 4. We point out that in the absence of an equivalent of the Kármán-Howarth relationship for the Lagrangian time increments, a limit value of $`C_0`$ is not a priori fixed. Dimensional analysis yields $`D_2^L(\tau )=C_0(Re)ϵ\tau `$ and similarity arguments give $`C_0(Re)\mathrm{const}.`$ or $`C_0(Re)Re^\alpha `$ in the limit of infinite Reynolds numbers.

To further describe the statistics of the Lagrangian velocity fluctuations, we have analyzed the statistics of the velocity increments $`\mathrm{\Delta }_\tau v`$. Their PDF $`\mathrm{\Pi }_\tau `$ for $`\tau `$ covering the accessible range of time scales is shown in Fig.3.

To emphasize the functional form, the velocity increments have been normalized by their standard deviation so that all PDFs have unit variance. A first observation is that the PDFs are symmetric, in agreement with the local symmetries this flow. Another is that the PDFs almost Gaussian at integral time scales and progressively develop stretched exponential tails for small time increments. At the smallest increment, the stretched exponential shape is in agreement with measurements of the PDF of Lagrangian acceleration at identical Reynolds numbers PortaBodenschatz01 . In our case, the limit form of the velocity increments PDF is not as wide as that of the acceleration because the Kolmogorov scale is not resolved. Note that in regards of the evolution of the PDF, the intermittency is at least as developed in the Lagrangian frame as it is in the Eulerian one Gagne .

The continuous evolution with scale can be quantified using the flatness factor. We show in Fig.4 the variation excess kurtosis $`K(\tau )=(\mathrm{\Delta }_\tau v)^4/(\mathrm{\Delta }_\tau v)^2^23`$. It is null at integral scale as expected from the Gaussian shape of the PDF and increases steeply at small scales. Below about $`5\tau _\eta `$, the increase is limited by the cut-off of the particle; an extrapolation of the trend to $`\tau _\eta `$ yields $`K(\tau _\eta )40`$ in agreement with acceleration measurements in PortaBodenschatz01 .

More generally, one can choose to describe the evolution of the PDFs by the behavior of their moments (or ‘structure functions’) $`D_q^L(\tau )=|\delta _\tau v|^q`$. Indeed, a consequence of the change of shape of the PDFs with scale is that their moments, as the flatness factor above, vary with scale. Classically in the Eulerian picture, one expects scaling in the inertial range, $`D_q^E(r)r^{\zeta _q}`$, at least in the limit of very large Reynolds numbers. At the finite Reynolds number where most experiments are made, the lack of a true inertial range is usually compensated by studying the relative scaling of the structure functions – the ESS ansatz Benzi . We use the second order structure function as a reference. Indeed the dimensional estimation of $`D_2^L`$ (as that of $`D_3^E`$) depends linearly on the increment and on the dissipation. Fig.5 shows that, as in the Eulerian frame, a relative scaling is observed for the Lagrangian structure functions of orders 1 to 5, $`D_q^L(\tau )D_2^L(\tau )^{\xi _q}`$. We observe that the relative exponents follow a sequence close to, but more intermittent than the corresponding Eulerian quantity. Indeed, we obtain: $`\xi _1^L/\xi _3^L=0.42`$, $`\xi _2^L/\xi _3^L=0.75`$, $`\xi _4^L/\xi _3^L=1.17`$, $`\xi _5^L/\xi _3^L=1.28`$ to be compared to the commonly accepted Eulerian values Arneodo96  $`\xi _1^E/\xi _3^E=0.36`$, $`\xi _2^E/\xi _3^E=0.70`$, $`\xi _4^E/\xi _3^E=1.28`$, $`\xi _5^E/\xi _3^E=1.53`$.

In conclusion, using a new experimental technique, we have obtained a Lagrangian velocity measurement that covers the inertial range of scales. Our results are consistent with Kolmogorov-like dimensional predictions for second order statistical quantities. At higher orders, the observed intermittency is very strong. How the Lagrangian intermittency is related to the statistical properties of the energy transfers is an open question. From a dynamical point of view, the Navier-Stokes equation in Lagrangian coordinates could be modeled using stochastic equations. Work is currently underway to compare the dynamics of the Lagrangian velocity to predictions of Langevin-like models.

acknowledgements: We thank Bernard Castaing for interesting discussions and Vermon Corporation for the design of the ultrasonic transducers. This work is supported by grant ACI No.2226 from the French Ministère de la Recherche.