skip to main content

gaia data release 3 documentation

7.6 Eclipsing binaries

7.6.3 Physical model

The EB simulator used for fitting the observational data (as explained in Section 7.6.4), generates synthetic multicolour light and radial velocity curves from a set of physical parameters. The radial velocity curves are required for the combined solutions described in Section 7.7. The simulator offers capabilities similar to those of the well-established Wilson-Devinney LC code (Wilson and Devinney 1971; Wilson 1979, 1990). It models stellar surfaces as equipotentials of the Roche potential, discretised using a mesh, and it can take into account the effects of limb darkening, gravity brightening, mutual irradiation, asynchronous rotation, the presence of a third light source, and the presence of spots (although the last two have not been enabled in Gaia DR3).

The Roche potential

The underlying physical model treats the surfaces of the two stellar components as equipotentials of an “effective” Roche potential which takes into account asynchronous rotation (Wilson 1979). The Roche geometry was favoured over simpler geometries because a significant fraction of the EBs detected by Gaia were expected to have short periods where proximity effects can be important. Furthermore, using the asynchronous version of the potential allows for a more general and precise treatment of stellar geometry in such systems. In particular, if the two components are i and j, the potential for component i is given by

Pi(x,y,z)=1r+qi[1d2-2dx+r2-xd2]+12d3fi2(1+e)(qi+1)(x2+y2). (7.43)

where r2=x2+y2+z2, qi=Mj/Mi is the mass ratio of the two components, d is their distance in units of the periastron distance, e is the orbital eccentricity, and fi=ωi/ω is the rotational parameter, i.e., the ratio of the angular speed ωi, assumed constant, of component i to the mean orbital angular speed ω. Stellar and orbital angular velocity vectors are presumed parallel.

It should be noted that in the case of an eccentric orbit, the distance between the components is not constant and the potential becomes time-dependent, and therefore non-conservative. Equation 7.43 is only valid under the assumption that the components can readjust to equilibrium on time scales that are short compared to the orbital period (Wilson 1979; Sepinsky et al. 2007). In that case, rotation will tend to synchronise to the angular rotation rate around periastron (Hut 1981), which leads to the following periastron-synchronised rotational parameter:

f2=1+e(1-e)3. (7.44)

When e0, the simulator uses this equation to set both f1 and f2 in Equation 7.43.

The orbit itself is calculated in the usual way, requiring as input parameters the eccentricity and the argument of periastron (ω). For a circular orbit (e=0), ω is set to zero. The projection of the orbit on the plane of the sky is determined by the orbital inclination (i), where i=90 corresponds to an orbital plane seen edge-on (central eclipse).

The simulator generates a synthetic LC by calculating the flux emitted by the EB in the direction of the observer at each of a requested (input) list of phases. The procedure, along with the requisite input parameters (in addition to the mass ratio and the rotational parameters mentioned above), is briefly described below. The specific values selected for those parameters for Gaia DR3 are given in Section 7.6.5.

The fill factors 𝐬1 and 𝐬2

The simulator takes as input the fill factor, s, of each component and determines the polar radius, rpole, of each component via the relation


where rcrit is the intersection of the polar axis of the component with the critical lobe of the system, i.e., with the equipotential surface that passes through the first Lagrange point, L1. According to this definition, a binary system is detached (the stars are separated) when the fill factors of both stars, s1 and s2, are <1; semi-detached (the stars separated but with one component filling its critical lobe) when either s1=1 and s2<1, or s1<1 and s2=1; or overcontact (the stars share a common envelope) when s1>1 or s2>1, in which case the fill factor of only one component may be specified, the fill factor of the other component being calculated using the value of rcrit that corresponds to the common equipotential surface of the system. The simulator checks that the fill factors do not exceed a value, dependent on the mass ratio, beyond which stellar mass is lost through the Lagrange point L2 and/or L3.

In the case of eccentric motion, the potential, and hence the fill factors, are time-dependent. In that case, the values of the fill factors are taken to refer to the time of periastron.

In the case of overcontact systems, circularity and synchronous rotation are enforced, i.e., e=0 and f1=f2=1.

The stellar surface mesh

The surface of each component is discretised using a mesh of quadrilateral surface elements of roughly equal areas, with special care taken in the representation of the tidally distorted regions of the stars. The mesh resolution of the larger component is determined by the the input parameter N, the number of requested surface elements. The number of surface elements for the smaller component is then determined automatically as (rs/rb)2N, where rs and rb are, respectively, the polar radii of the smaller and of the larger component. In this way, the area of the surface elements remains similar in each component. However, the number of surface elements of the smaller component is not allowed to become smaller than the value of input parameter Nmin. During eclipse calculations, the elements of the eclipsed component located near the limb of the eclipsing component contribute light according to their partial visibility, i.e., the visible fraction of their surface, in order to improve LC accuracy and reduce discretisation effects. In the case of eccentric motion, the surface meshes are recalculated at each requested phase so that the volume of each star is conserved along its orbit.

Flux calculation

Figure 7.52: The eclipsing binary simulator represents the surfaces of the two stellar components using meshes that trace the corresponding equipotential surfaces of the Roche model. The top row uses sequential colour maps of varying lightness to represent the flux emitted from the stellar surfaces at two orbital phases (a) and (b). The bottom row uses rainbow colour maps to depict the departures of the local temperature from the effective temperature of each component due to gravity brightening and mutual irradiation effects. The stellar surface meshes are also visible in the bottom row. For the larger component, Teff,1=7000K and s1=1, while for the smaller component Teff,2=5000K and s2=0.98. The mass ratio q=M2/M1=0.5, and the number of surface elements of the larger component is N=3000.

The total flux emitted by the system at a given photometric passband is the sum of the fluxes emitted in the direction of the observer by the visible surface elements. The calculation of the flux of each element depends on the choice of a stellar atmosphere model. The simulator calculates fluxes either for a blackbody atmosphere, where the flux depends only on the effective surface temperature (Teff), or for an arbitrary stellar atmosphere model provided to the simulator in a tabular form listing the flux as a function of temperature, metallicity and surface gravity (log g). In the latter case, a blackbody law is assumed for temperatures outside the temperature range of the tabulated data, and a scaling is applied so that, at each end of the range, the blackbody fluxes are identical to those provided by the atmosphere models. For the Gaia DR3 run, a blackbody law has been used, as explained in Section 7.6.4.

The simulator provides a choice of limb darkening laws, namely linear, quadratic and square-root. When enabled, gravity brightening is calculated as a function of the effective temperature if it is below the radiative limit (set to 7 700K), or set to 1 above it. The calculation of mutual irradiation between the two components, if enabled, implements a correction to the effective temperature of the surface elements concerned. The simulator offers the choice between a simple and a more accurate but slower treatment of this effect. In the simple treatment, the temperature correction is computed independently for each component without taking into account the temperature increase that it induces on the other star. By contrast, the accurate treatment takes into account the fact that the temperature increase of the irradiated star leads to an increased flux which, in return, will increase the temperature of certain surface elements of the irradiating star. This calculation is iterated until the relative difference of the average temperature over all the affected elements is below a threshold value (set to 0.005). Usually convergence is achieved after a few iterations. Figure 7.52 is a visualisation of the Roche model as implemented by the simulator for a particular choice of physical parameters.