© ESO, 2012

1. Introduction

The distortions in the shape of binary stars due to proximity effects cause a slow rotation of the line of the apses, known as apsidal-motion, that provides a robust tool for searching into the interior of the stars. These distortions can be calculated as a function of the internal structure of the stars (Kopal 1978). An accurate measurement of this phenomenon can be obtained by monitoring the primary and secondary minima of eclipsing binaries. The resulting empirical apsidal-motion rates can be compared with the theoretical ones, derived from evolutionary stellar models. For an update of the comparison theory-observations, we refer to Claret & Giménez (2010) and, for relativistic systems, sto Claret (1997) and more recently Wolf et al. (2010).

Tidal evolution can also be used to investigate the interior of the stars. The available approximations to the differential equations that govern the orbital elements depend on the moment of inertia of each star and in some formalisms, also the apsidal-motion constants (see for example, Hut 1981; Eggleton & Kiseleva-Eggleton 2002). One of the first calculations of a stellar moment of inertia was carried out by Motz (1952). However, the adopted models at that time were too simple. The first attempt to produce a systematic calculation of the fractional radius of gyration was that of Claret & Giménez (1998), who tabulated the fractional radius of gyration for a grid of hydrogen-burning stars for use in stellar dynamics. Later, in Claret & Giménez (1991), the same authors investigated the influence of core overshooting and mass loss on the internal structure constants.

In relation still to stellar distortions, the computation of the theoretical light curves of eclipsing binaries, planetary transits, and single rapidly rotating stars requires the knowledge of how the flux is distributed over the distorted stellar surfaces. There are two important historical papers that provide the values of the gravity-darkening exponents (GDE) β₁: that of von Zeipel (1924), which treats stars with radiative envelopes (β₁ = 1.0), and that of Lucy (1967), who used convective envelopes, which are appropriate for late-type stars, to derive an average of β₁ = 0.32. These values were used over several decades until the paper of Claret (1998), who introduced a numerical method based on the triangles strategy of Kippenhahn et al. (1967, which provides values of the gravity-darkening exponents as a function of age, mass, and chemical composition. A smooth transition in the values of β₁ is achieved between the convective/radiative transport of energy.

All the aforementioned parameters are presently available for several masses and chemical compositions covering evolutionary phases such as the main sequence (MS) and red giant branch but for early phases the calculations are very scarce or even lacking. In the present paper, we present the calculations of the apsidal-motion constants, the fractional radius of gyration, and the gravity-darkening exponents for an extensive grid of PMS. These calculations can be used to investigate the apsidal-motion of double-lined eclipsing binaries, the tidal evolution, and the brightness distribution of very young stars.

2. The stellar models

The code used to generate the PMS models is essentially the same as that described in Claret (2004). For completeness, we give some of its main characteristics describing PMS evolution. We adopt the equation of state OPAL (Rogers et al. 1996). The code uses the opacity tables provided by Iglesias & Rogers (1996), which are completed by the calculations by Ferguson et al. (2005) for lower temperatures. We adopt X = 0.70 and Z = 0.02, which is representative of a typical chemical composition. Two additional compositions are available of (X,Z) = (0.757, 0.001) and (0.64, 0.04).

Convection was treated by using the mixing-length theory (α_MLT = 1.5) and we adopt the Schwarzschild criterion to delineate the convective boundaries. To get more accurate computations we decreased the triangle size used to define an envelope in Herzsprung-Russel (HR) diagram to Δlog T_eff = 0.001 and Δlog L = 0.004. This is particularly important to determinate the depth of the convective envelope. The calculations cover the mass range 0.05–30 M_⊙ (Table 1). All models were computed using our stellar evolution code (Claret 2004), except those with masses 0.05, 0.1, 0.2, and 0.4 M_⊙ for which we used the MESA code (Paxton et al. 2011). Figure 1 shows the HR diagram for the PMS models with (X,Z, α_MLT) = (0.70, 0.02, 1.5).

3. Gravity-darkening exponents, apsidal-motion constants, radius of gyration and potential energy

Fig. 1 Evolutionary tracks for PMS models. Attached numbers denote the masses of the model in solar units.

We introduce the equations and numerical methods necessary to calculate the gravity-darkening exponents, the radius of gyration, the potential energy, and the internal structure constants. The gravity-darkening exponent is an important parameter that is often used in the analysis of the light curves of eclipsing binaries and the modelling of rapidly rotating stars. For stars in strict radiative equilibrium (pseudo-barotrope), von Zeipel (1924) showed that the variation in the brightness over the surface is proportional to the effective gravity: $\mathit{F}\mathrm{=}\mathrm{-}\frac{\mathrm{4}\mathit{ac}{\mathit{T}}^{\mathrm{3}}}{\mathrm{3}\mathit{\kappa \rho }}\frac{\mathrm{d}\mathit{T}}{\mathrm{d\Psi }}{\mathit{g}}^{{\mathit{\beta }}_{\mathrm{1}}}\mathit{,}$(1)or ${{\mathit{T}}_{\mathrm{eff}}}^{\mathrm{4}}\mathrm{\propto }{\mathit{g}}^{{\mathit{\beta }}_{\mathrm{1}}}\mathit{,}$(2)where β₁ = 1.0, g is the local gravity, Ψ is the total potential, and the remaining variables have their usual meaning in studies of stellar structure. Lucy (1967) computed values of β₁ for convective envelopes, finding an average value of 0.32. As mentioned in the introduction, we have introduced a numerical method to compute β₁, that allows us to calculate the gravity-darkening exponents as a function of age, mass, chemical composition, etc. Such a method provides a set of β₁ that supersedes the old values of 0.32 and 1.0 for convective and radiative envelopes respectively and a smooth transition is achieved between both energy transport mechanisms. The aforementioned smooth transition zone was confirmed by the observations of Che et al. (2011) for very rapid rotators.

The numerical method is based on the triangles strategy (Claret 1998). Given that the external boundary conditions are unchanged during the evolution, the integrations in the outer layers will be the same, provided that the corresponding triangles in the HR diagram are small. By generalizing to several points, we can then simulate distorted configurations with different effective temperature distributions over the surface. The numerical differentiation of these envelopes allows us, finally, to obtain β₁. A variation in this method was presented in Claret (2012), where the GDE were computed as a function of the optical depth. In that paper, we showed that the classical von Zeipel’s theorem (β₁ = 1.0) is not strictly valid in the upper layers of hot rotating stars.

We derive the theoretical apsidal-motion constants k_j by integrating the differential equations of Radau $\frac{\mathit{a}\mathrm{d}{\mathit{\eta }}_{\mathrm{j}}}{\mathrm{d}\mathit{a}}\mathrm{+}\frac{\mathrm{6}\mathit{\rho }\mathrm{\left(}\mathit{a}\mathrm{\right)}}{\overline{\mathit{\rho }}\mathrm{\left(}\mathit{a}\mathrm{\right)}}\mathrm{(}{\mathit{\eta }}_{\mathrm{j}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{+}{\mathit{\eta }}_{\mathrm{j}}\mathrm{(}{\mathit{\eta }}_{\mathrm{j}}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{=}\mathit{j}\mathrm{(}\mathit{j}\mathrm{+}\mathrm{1}\mathrm{)}\mathit{,} \mathit{j}\mathrm{=}\mathrm{2}\mathit{,}\mathrm{3}\mathit{,}\mathrm{4}\mathit{,}$(3)where ${\mathit{\eta }}_{\mathrm{j}}\mathrm{\equiv }\frac{\mathit{a}}{{\mathit{ϵ}}_{\mathrm{j}}}\frac{\mathrm{d}{\mathit{ϵ}}_{\mathrm{j}}}{\mathrm{d}\mathit{a}}$(4)and a is the mean radius of the configuration, ϵ_j is the deviation from sphericity, ρ(a) is the mass density at the distance a from the centre, and $\overline{\mathit{\rho }}\mathrm{\left(}\mathit{a}\mathrm{\right)}$ is the mean mass density within a sphere of radius a. The k_j are computed using ${\mathit{k}}_{\mathrm{j}}\mathrm{=}\frac{\mathit{j}\mathrm{+}\mathrm{1}\mathrm{-}{\mathit{\eta }}_{\mathrm{j}}\mathrm{\left(}\mathit{R}\mathrm{\right)}}{\mathrm{2}\mathrm{(}\mathit{j}\mathrm{+}{\mathit{\eta }}_{\mathrm{j}}\mathrm{(}\mathit{R}{\mathrm{\left)}}^{\mathrm{\right)}}}\mathit{,}$(5)where R refers the values of η_j at the stellar surface.

The fractional radius of gyration β is obtained for each stellar configurations using $\mathit{I}\mathrm{=}\frac{\mathrm{8}\mathit{\pi }}{\mathrm{3}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{R}}\mathit{\rho }\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathit{r}}^{\mathrm{4}}\mathrm{d}\mathit{r}\mathrm{=}\mathit{M}\mathrm{(}\mathit{\beta R}{\mathrm{)}}^{\mathrm{2}}$(6)and the potential energy $\mathrm{\Omega }\mathrm{=}\mathrm{-}\mathit{G}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{M}}\frac{\mathit{m}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathit{r}}\mathrm{d}\mathit{m}\mathrm{=}\mathrm{-}\frac{\mathit{\alpha G}{\mathit{M}}^{\mathrm{2}}}{\mathit{R}}\mathit{,}$(7)where the form-factors β and α describe the moment of inertia and the potential energy, respectively. Equations (3), (6), and (7) were integrated using a fourth order Runge-Kutta method and the values of k_j, β, and α were finally derived.

Fig. 2 The gravity-darkening exponent evolution. The continuous line represents a 2 M⊙ during the PMS phase, while the dashed line denotes the same model after the PMS stage. The two arrows indicate the direction of the evolution in both cases.

Fig. 3 The apsidal-motion constant log k2 as a function of mass and log g for some selected PMS models. The attached numbers indicate the model mass in solar units.

A contracting model with 2 M_⊙ is very useful to illustrate the evolution of β₁ (Fig. 2). The continuous line represents the PMS-ZAMS phase and the dashed one denotes the post ZAMS models. The arrows indicate the time evolution direction. In the first stages of PMS, the depth of the convective envelope is large and decreases progressively as the star approaches the ZAMS. The corresponding β₁ is almost constant in this range. Outside this interval, this value increases, reaching a local maximum around log T_eff = 3.70. This behaviour is due to the onset of the convection: the convection is present even in thin layers. For effective temperatures lower than 5000 K, there is a decrease in the optical depth where convection appears. The hydrogen recombination is responsible for this and the efficiency of the convection is higher. For higher effective temperatures, convection loses importance and the corresponding β₁ decreases until a local minimum at log T_eff ≈ 3.84 (that value depends on the envelope input physics). In addition, the von Zeipel value of β₁ is progressively restored towards the ZAMS point.

During the PMS, β₁ is systematically larger than the computed ones for later phases of log T_eff < 3.84. Although the PMS and post-ZAMS models show similar effective temperatures, the respective physical conditions in the envelopes are different. When the depth of the convective envelopes are of the same order, β₁ converges to similar values. This can be seen in particular at the beginning of the PMS or for models near the ZAMS.

Our codes which are designed to compute the GDE, do not include ingredients such as the evolution of the chemical abundance distribution caused by gravitational settling, chemical/thermal diffusion of the elements, or a more specific equation of state for very compact objects and therefore, for consistency, we do not calculate the gravity-darkening exponents for very low mass stars. These objects, owing to their peculiar characteristics, will be studied in a separate paper.

Figure 3 shows the evolution of log k₂ as a function of log g for some selected PMS sequences. The masses (in solar units) are indicated in the figure. Two interesting features can be seen in Fig. 3. The first one is related to the first stages of the PMS tracks: the models have similar mass distributions which correspond to configurations with an effective polytropic index 1.5–2.0.

As the models evolve, they become more centrally condensed and, consequently, the apsidal-motion constant decreases. The second point refers to the double maxima in log k₂, which appears for models with initial masses higher than 1.2 M_⊙. These changes in the internal structure of the models when they are close to the ZAMS are due to the reduction in the original C¹² through the nuclear reactions C¹²(p, γ)N¹³(β ⁺ , ν)C¹³(p, γ)N¹⁴. The models have a nuclear energy source, though the gravitational and thermal energies also contribute. The mass concentration either reduces or increases depending on the predominant source of energy. A similar characteristic is also detected in later evolutionary phases (for m > 1.2   M_⊙) and is due to a change in the dominant nuclear source of energy from the proton-proton chain to CNO cycle.

Fig. 4 Comparison between apsidal-motion constants log k2: with our stellar evolutionary code (continuous line) and that computed by using MESA (dashed line). The figure in the upper right corner shows a comparison of the corresponding HR diagram for the same models. The characteristics of the models are 3 M⊙, (X,Z, αMLT) = (0.70, 0.02, 1.5).

Figure 4 compares the log k₂ for a 3 M_⊙ model computed using our stellar evolution code (Claret 2004, continuous line) with the evolution of log k₂ computed with the MESA code (version 3290) recently developed by Paxton et al. (2011) denoted

by the dashed line. The agreement between the two calculations is very good in all phases of the PMS evolution. We also found very good agreement for later phases such as the MS or red giant branch. The figure at the upper right corner of Fig. 4 shows the corresponding HR diagram for the same models. The agreement can be considered to be good; the maximum differences in T_eff are of the order of 1–3% in the either early stages or near ZAMS.

All calculations presented here for α and β have been performed for spherical configurations. However, in the case of elliptical symmetry and an ellipsoidal mass distribution the corresponding α and β can also be computed by writing them in terms of both the spherical structure and of a correction factor that depends on the geometry of the perturbation. We return to this specific subject in future investigations.

Finally, we note that the PMS tables contain 11 columns corresponding to the age (in years), the mass (solar units), log L (in solar units), log T_eff, log g(cgs), log k₂, log k₃, log k₄, α, β, and β₁.

Acknowledgments

B. Paxton and A. Dotter are particularly acknowledged for their help in the implementating MESA. The Spanish MEC (AYA2009-10394, AYA2009-14000-C03-01) is gratefully acknowledged for its support during the development of this work. This research has made use of the SIMBAD database, operated at the CDS, Strasbourg, France, and of NASA’s Astrophysics Data System Abstract Service.

References

All Tables

All Figures

	Fig. 1 Evolutionary tracks for PMS models. Attached numbers denote the masses of the model in solar units.
In the text

	Fig. 2 The gravity-darkening exponent evolution. The continuous line represents a 2 M⊙ during the PMS phase, while the dashed line denotes the same model after the PMS stage. The two arrows indicate the direction of the evolution in both cases.
In the text

	Fig. 3 The apsidal-motion constant log k2 as a function of mass and log g for some selected PMS models. The attached numbers indicate the model mass in solar units.
In the text

	Fig. 4 Comparison between apsidal-motion constants log k2: with our stellar evolutionary code (continuous line) and that computed by using MESA (dashed line). The figure in the upper right corner shows a comparison of the corresponding HR diagram for the same models. The characteristics of the models are 3 M⊙, (X,Z, αMLT) = (0.70, 0.02, 1.5).
In the text