7.6.3 Physical model

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

where $r^2=x^2+y^2+z^2$, $q_i=M_j/M_i$ is the mass ratio of the two components, $d$ is their distance in units of the periastron distance, $e$ is the orbital eccentricity, and $f_i={\omega }_i/\omega$ is the rotational parameter, i.e., the ratio of the angular speed ${\omega }_i$, assumed constant, of component $i$ to the mean orbital angular speed $\omega$. 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:

When $e\ne 0$, the simulator uses this equation to set both $f_1$ and $f_2$ in Equation 7.43.

The orbit itself is calculated in the usual way, requiring as input parameters the eccentricity and the argument of periastron ($\omega$). For a circular orbit ($e=0$), $\omega$ 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}^{\circ }$ 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 simulator takes as input the fill factor, $s$, of each component and determines the polar radius, $r_{\mathrm{pole}}$, of each component via the relation

where $r_{\mathrm{crit}}$ 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, $L_1$. According to this definition, a binary system is detached (the stars are separated) when the fill factors of both stars, $s_1$ and $s_2$, are ; semi-detached (the stars separated but with one component filling its critical lobe) when either $s_1=1$ and , or and $s_2=1$; or overcontact (the stars share a common envelope) when $s_1>1$ or $s_2>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 $r_{\mathrm{crit}}$ 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 $L_2$ and/or $L_3$.

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 $f_1=f_2=1$.

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 ${(r_s/r_b)}^{2}N$, where $r_s$ and $r_b$ 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 $N_{\min }$. 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.

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 ($T_{\mathrm{eff}}$), 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 $\mathrm{7\hspace{0.17em}700}K$), 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.