Radiative Properties of Binary Stars

The distribution of intensity over the stellar disk depends on the dominant energy transfer mechanism in a stellar envelope. The mechanism is primarily determined by the effective temperature of the star: for ${T}_{\mathrm{eff}}\gt 8000$ K the envelope is in a predominantly radiative equilibrium, whereas for ${T}_{\mathrm{eff}}\lt 5000$ K the envelope is in a predominantly convective equilibrium (Claret 1999); in between, both mechanisms significantly contribute to the energy transfer.

For a rotating star in a strict radiative equilibrium, a pioneering work by von Zeipel (1924) has shown that the flux distribution over the surface is proportional to the effective gravity:

$$\begin{eqnarray}&&{F}_{\lambda }=-\frac{16\sigma {T}^{3}}{3\bar{\kappa }\rho }\frac{{dT}}{d{\rm{\Omega }}}{g}^{\beta },\end{eqnarray} \tag{ 5.1 }$$

where $\sigma =5.67\times {10}^{-8}\,{\rm{W}}\,{{\rm{m}}}^{-2}\,{{\rm{K}}}^{-4}$ is the Stefan–Boltzmann constant, T is the local surface temperature, $\bar{\kappa }$ is the Rosseland mean opacity, ρ is the density of the gas, g is the local gravity acceleration, and β is the gravity-darkening coefficient, for which von Zeipel (1924) rigorously showed that $\beta =1$ for strictly radiative envelopes. Equation (5.1) implies that the equatorial regions are darker than the polar ones, hence the name "gravity darkening." The local temperature may also be coupled to the gravitational acceleration:

$$\begin{eqnarray}&&T\propto {F}^{1/4}\Rightarrow T\propto {g}^{\beta /4}.\end{eqnarray} \tag{ 5.2 }$$

The generalized formalism of computing the value of the gravity-darkening coefficient β for radiative envelopes was presented later by Kippenhahn (1977).

Prediction for the gravity-darkening exponent β of purely convective envelopes was first done by Lucy (1967), who used convective stellar envelope models to compute the average value $\beta =0.32$ for solar-type stars. The driving idea of Lucy's study was that, within the fully convective zone, the gradient of specific entropy is 0. His predictions for the gravity-darkening coefficient were later confirmed observationally by Rafert & Twigg (1980), who obtained the average value $\beta =0.31$ for the sample of 28 stars.

Putting the two results for radiative and convective equilibrium together, gravity darkening is given by:

$$\begin{eqnarray}F={F}_{\mathrm{pole}}\,{\left(\frac{g}{{g}_{\mathrm{pole}}}\right)}^{\beta },\quad \beta =\left\{\begin{array}{ll}1.00 & {\rm{for\; radiative\; envelopes}}\\ 0.32 & {\rm{for\; convective\; envelopes}}.\end{array}\right.\end{eqnarray} \tag{ 5.3 }$$

Factor

$$\begin{eqnarray}&&{ \mathcal G }=\frac{F}{{F}_{\mathrm{pole}}}={\left(\frac{g}{{g}_{\mathrm{pole}}}\right)}^{\beta }={T}_{\mathrm{pole}}\,{\left(\frac{g}{{g}_{\mathrm{pole}}}\right)}^{\beta /4}\end{eqnarray} \tag{ 5.4 }$$

is called the gravity-darkening correction.

Two later studies by Alencar & Vaz (1997) and Alencar et al. (1999) followed Lucy's approach of calculating gravity darkening for convective stellar atmospheres; they used nonilluminated and illuminated Uppsala model atmospheres (Gustafsson 1983), respectively, deriving an empirical formula for $\beta ({T}_{\mathrm{eff}})$, depicted in Figure 5.1 (left panel). The results indicate agreement with Lucy's estimate only at temperatures $\sim \,6500\,{\rm{K}}$.

Zoom In Zoom Out Reset image size

Since the gravity-darkening effect depends on the circumstances in stellar envelopes, it is related not only to the atmospheric parameters (${T}_{\mathrm{eff}}$ and $\mathrm{log}\,g$) but also to the internal stellar structure and rotation properties. This was pointed out by Claret (1999), who introduced the computation of gravity-darkening coefficients into the stellar evolution code. Statistical studies based on this new formalism (Claret 2000, Figure 1 therein) yielded a more accurate theoretical prediction for stellar envelopes where both radiation and convection contribute significantly to energy transfer.

The gravity-darkening correction given by Equation (5.4) is bolometric. Thus, to apply gravity darkening appropriately to passband data, a bolometric correction needs to be made. Such a correction has been first proposed by Martynov (1973) and later expanded by Bloemen et al. (2011) to the following form:

$$\begin{eqnarray}&&\beta (\lambda ;{T}_{\mathrm{eff}},\mathrm{log}\,g,\ldots )={\left(\frac{\partial \mathrm{log}I(\lambda )}{\partial \mathrm{log}g}\right)}_{{T}_{\mathrm{eff}}}+\left(\frac{\partial \,\mathrm{log}\,{T}_{\mathrm{eff}}}{\partial \,\mathrm{log}\,g}\right){\left(\frac{\partial \mathrm{log}I(\lambda )}{\partial {{\rm{T}}}_{\mathrm{eff}}}\right)}_{\mathrm{log}g}.\end{eqnarray} \tag{ 5.5 }$$

The right panel of Figure 5.1 depicts the dependence of β on temperature for the Kepler passband (Claret & Bloemen 2011).

Espinosa Lara & Rieutord (2011) also pointed out the failure of von Zeipel's theorem for rapidly rotating stars. They proposed a model that describes latitudinal variation of gravity darkening as a consequence of rotation. In their subsequent paper, Espinosa Lara & Rieutord (2012) further developed the formalism for gravity darkening in binary stars and advocated for the abandonment of a single-parameter description and replacing it with the following expression:

$$\begin{eqnarray}&&f(x,y,z)={f}_{0}\,\exp \,\left({\int }_{0}^{{r}_{\mathrm{point}}}\frac{2(1+q)}{{g}_{\mathrm{eff}}^{2}(x,y,z)}{g}_{\mathrm{eff}}(x,0,0){dx}\right),\end{eqnarray} \tag{ 5.6 }$$

where f(x, y, z) is surface flux at point (x, y, z), ${f}_{0}={\mathrm{lim}}_{r\to 0}\,f=L/4\pi {GM}$, ${{\boldsymbol{g}}}_{\mathrm{eff}}=-{\rm{\nabla }}\psi$ is the effective gravitational acceleration, and g_eff is the magnitude of ${{\boldsymbol{g}}}_{\mathrm{eff}}$. This integral can be solved numerically to obtain a relationship between the flux and gravitational acceleration.

When components in a binary system are close, i.e., when their radii are ∼ 15%–20% of their separation or more (Wilson 1990), surface irradiation of one component by the other becomes important. A part of the incident flux heats the surface of the irradiated star, increasing the temperature, and the other part of the incident flux is reflected off the irradiated star. This effect is called the reflection effect.

Before we develop the formalism to account for the reflection effect, we should identify the cases where reflection contributes significantly to the overall flux of the irradiated component. There are indeed quite a number of such cases, ranging from stars similar in size and effective temperature (Hilditch et al. 1996), to cataclysmic and pre-cataclysmic binaries (Pigulski & Michalska 2002), to rapidly rotating Algols (Olson & Etzel 1994), to elliptical ellipsoidal variables known as heartbeat stars (Welsh et al. 2011).

The rigorous treatment of the reflection effect in binary stars was first developed by Wilson (1990); before that work only approximate solutions existed—see a thorough review by Vaz (1985) on this topic.

Consider a close binary as depicted in Figure 5.2. To properly account for the heating of a given surface element dA₁ on the primary star, we need to sum all contributions to irradiation over all surface elements dA₂ on the secondary star. Following the energy conservation principle, recall that heating is a bolometric process, whereas radiation is wavelength dependent, described by the spectral energy distribution (SED) function (blackbody, model atmospheres, ...). Bolometric emergent intensity ${I}_{2}^{\mathrm{bol}}$ emitted from surface element dA₂ into the direction of dA₁ may be expressed as the limb-darkened normal bolometric emergent intensity ${I}_{2\perp }^{\mathrm{bol}}$:

$$\begin{eqnarray}&&{I}_{2}^{\mathrm{bol}}={I}_{2\perp }^{\mathrm{bol}}{{ \mathcal L }}_{2}^{\mathrm{bol}}({\mu }_{2}),\end{eqnarray} \tag{ 5.7 }$$

where ${\mu }_{2}=\cos \,{\psi }_{2}$ is the cosine of the angle between the surface normal to dA₂ and the direction ${\boldsymbol{s}}$ that connects the two surface elements. The incident intrinsic flux from dA₂ to dA₁ can then be written per Equation (4.15) as:

$$\begin{eqnarray}&&{{dF}}_{2\to 1}^{\mathrm{bol}}=\frac{{\mu }_{1}{\mu }_{2}}{{s}^{2}}{I}_{2\perp }^{\mathrm{bol}}{{ \mathcal L }}_{2}^{\mathrm{bol}}\,({\mu }_{2}){{dA}}_{2},\end{eqnarray} \tag{ 5.8 }$$

where ${\mu }_{1}=\cos \,{\psi }_{1}$ is the cosine of the incident angle between the surface normal to dA₁ and the direction $-{\boldsymbol{s}}$.

Zoom In Zoom Out Reset image size

We now introduce a reflection factor ${ \mathcal R }$ as the ratio of the total (intrinsic + reflected) flux to intrinsic flux. If there is no reflection, ${ \mathcal R }\equiv 1$. The total incident flux from dA₂ to dA₁ can then be computed by multiplying the intrinsic incident flux, given by Equation (5.8), by ${ \mathcal R }$ ₂:

$$\begin{eqnarray}&&{F}_{2\to 1}^{\mathrm{bol}}={\int }_{\partial V}{{ \mathcal R }}_{2}{ \mathcal V }({\boldsymbol{r}})\frac{{\mu }_{1}{\mu }_{2}}{{s}^{2}}{I}_{2\perp }^{\mathrm{bol}}{{ \mathcal L }}_{2}^{\mathrm{bol}}({\mu }_{2})\,{{dA}}_{2},\end{eqnarray} \tag{ 5.9 }$$

where we used the visibility function ${ \mathcal V }({\boldsymbol{r}})$ defined by Equation (4.13), this time pertaining to the visibility of surface elements along the line s.

A fraction of the total incident flux will be used for heating the irradiated surface, and a fraction will be reflected. The two fractions need to sum up to 1, and they are determined by the bolometric albedo of the star, ${ \mathcal A }$. It is the ratio of the reradiated flux to irradiated flux. Thus, the fraction of reflected flux is given by:

$$\begin{eqnarray}&&{F}_{1,\mathrm{refl}}^{\mathrm{bol}}={{ \mathcal A }}_{1}{\int }_{\partial V}{{ \mathcal R }}_{2}{ \mathcal V }({\boldsymbol{r}})\frac{{\mu }_{1}{\mu }_{2}}{{s}^{2}}{I}_{2\perp }^{\mathrm{bol}}{{ \mathcal L }}_{2}^{\mathrm{bol}}({\mu }_{2})\,{{dA}}_{2}.\end{eqnarray} \tag{ 5.10 }$$

The intrinsic bolometric normal emergent intensity ${I}_{1\perp }^{\mathrm{bol}}$ emitted from the surface element dA₁ is computed from the polar emergent intensity ${I}_{\mathrm{polar}}^{\mathrm{bol}}$ corrected for gravity darkening via Equation (5.4). Note that Equation (5.4) relates fluxes and not intensities, but as these are bolometric intensities (integrated over all wavelengths), it means that ${F}_{\mathrm{bol}}=\pi {I}_{\mathrm{bol}}$ per Stefan–Boltzmann's law, irrespective of the actual distribution of I(λ). In other words, all models for I(λ) (blackbody, Kurucz, ...), when integrated over all λ, must yield $\sigma {T}_{\mathrm{eff}}^{4}/\pi$.

The intrinsic bolometric flux ${F}_{1,\mathrm{intr}}^{\mathrm{bol}}$ is computed per Equation (4.5):

$$\begin{eqnarray}&&{F}_{1,\mathrm{intr}}^{\mathrm{bol}}=2\pi {I}_{1\perp }^{\mathrm{bol}}\int {{ \mathcal L }}_{1}^{\mathrm{bol}}({\mu }_{1})\,{\mu }_{1}\,d{\mu }_{1}.\end{eqnarray} \tag{ 5.11 }$$

The ratio ${F}_{1,\mathrm{refl}}^{\mathrm{bol}}/{F}_{1,\mathrm{intr}}^{\mathrm{bol}}$ is easily related to the reflection factor ${{ \mathcal R }}_{1}$:

$$\begin{eqnarray}&&{{ \mathcal R }}_{1}=\frac{{F}_{1,\mathrm{intr}}^{\mathrm{bol}}+{F}_{1,\mathrm{refl}}^{\mathrm{bol}}}{{F}_{1,\mathrm{intr}}^{\mathrm{bol}}}=1+\frac{{F}_{1,\mathrm{refl}}^{\mathrm{bol}}}{{F}_{1,\mathrm{intr}}^{\mathrm{bol}}}.\end{eqnarray} \tag{ 5.12 }$$

Similarly, for the other star,

$$\begin{eqnarray}&&{{ \mathcal R }}_{2}=\frac{{F}_{2,\mathrm{intr}}^{\mathrm{bol}}+{F}_{2,\mathrm{refl}}^{\mathrm{bol}}}{{F}_{2,\mathrm{intr}}^{\mathrm{bol}}}=1+\frac{{F}_{2,\mathrm{refl}}^{\mathrm{bol}}}{{F}_{2,\mathrm{intr}}^{\mathrm{bol}}}.\end{eqnarray} \tag{ 5.13 }$$

Equations (5.9), (5.10), and (5.12) thus form a recursive set to compute ${{ \mathcal R }}_{1}$, and their counterparts for the secondary star form a recursive set to compute ${{ \mathcal R }}_{2}$. In practice, this means that we start from some initial values for ${{ \mathcal R }}_{1}$ and ${{ \mathcal R }}_{2}$, say, unity, and then iterate the recursion set until we reach a desired accuracy.

Once the converged values of ${{ \mathcal R }}_{1}$ and ${{ \mathcal R }}_{2}$ are computed, we can calculate a new effective temperature of the surface element that encompasses the radiation of both intrinsic and reflected energy:

$$\begin{eqnarray}&&{{T}_{\mathrm{eff}}}_{1,\mathrm{new}}={{ \mathcal R }}_{1}^{1/4}{{T}_{\mathrm{eff}}}_{1,\mathrm{old}},\quad \mathrm{and}\quad {{T}_{\mathrm{eff}}}_{2,\mathrm{new}}={{ \mathcal R }}_{2}^{1/4}{{T}_{\mathrm{eff}}}_{2,\mathrm{old}}.\end{eqnarray} \tag{ 5.14 }$$

We can do this only because all quantities that partake in reflection are bolometric.

The above procedure may be repeated iteratively any number of times to accommodate for multiple reflections: the total emergent flux from the primary star now heats the secondary star by an amount greater than what the intrinsic emergent flux by itself would cause.

A keen eye will catch that this description neglects the fraction of the incident flux that is not reflected. This may or may not be a problem: if we assume that the absorbed energy is redistributed across the entire star, the effect will not be significant. However, if we assume that the absorbed energy is localized to the immediate vicinity of the irradiated area, then we expect an additional local temperature increase. This was first worked out by Prša et al. (2016). The basic difference in the treatment is the energy balance equation (see Figure 5.3):

$$\begin{eqnarray}&&{F}_{\mathrm{in}}({\boldsymbol{r}})={\hat{{ \mathcal L }}}_{\mathrm{LD}}{F}_{0}({\boldsymbol{r}})+{\hat{{ \mathcal L }}}_{{\rm{L}}}(\rho {F}_{\mathrm{in}})({\boldsymbol{r}}),\end{eqnarray} \tag{ 5.15 }$$

$$\begin{eqnarray}&&{F}_{\mathrm{out}}({\boldsymbol{r}})={F}_{0}({\boldsymbol{r}})+\rho ({\boldsymbol{r}}){F}_{\mathrm{in}}({\boldsymbol{r}}),\end{eqnarray} \tag{ 5.16 }$$

where $\hat{{ \mathcal L }}$ is the radiosity operator that describes the physical process of heating and reflection/scattering. The two simplest processes are Lambertian, described by ${\hat{{ \mathcal L }}}_{{\rm{L}}}={ \mathcal R }/\pi$, and limb-darkened distribution, described by ${\hat{{ \mathcal L }}}_{\mathrm{LD}}={ \mathcal R }\,\circ \,{ \mathcal L }$. An expanded reflection theory has been developed by Horvat et al. (2018), which includes different energy redistribution modes: local, latitudinal, and global (Figure 5.4). Their approach further generalizes the treatment of reflection, but it comes at a significantly higher computational cost.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The discovery of chromospheric activity of the Sun dates back to prehistoric times (ancient China and ancient Greece), but it was not until Galileo & Harriot observed them in 1610 and until Fabricius published his observations in 1611 that the world truly embraced the concept of sunspots, challenging the perfect, unchanging nature of the universe. Proposing spots on other stars followed shortly, by Boulliau (1667), but we had to wait for actual observations of spots until Kron (1947), who observed eclipsing binary AR Lac and invoked starspots to explain aperiodic out-of-eclipse variability. Over the next several decades we learned that starspots are connected to stellar magnetic fields (Spruit & Weiss 1986), so we expect most (if not all) convective envelope stars to feature spots. Another type of spot is related to mass transfer and accretion in cataclysmic variables (Warner 2003), which are hotter and brighter than the underlying photosphere. A comprehensive review of our current understanding of starspots is given by Strassmeier (2009).

Modern observations of starspots abound; Figure 5.5 depicts a light curve of KIC 2438502, an eclipsing binary with an 8.36-day period that features an RS CVn-type component that rotates with a period of ∼8.23 days. Most out-of-eclipse variation in this light curve originates in chromospheric activity, which is not periodic; hence, it results in the smearing of the phased light curve. Starspot evolution of two persistent active longitudes for this object has been reported and studied by Roettenbacher et al. (2016).

Zoom In Zoom Out Reset image size

In terms of modeling, spots may be described by four parameters:

• SPOT COLATITUDE ${\theta }_{\mathrm{sp}}$: the colatitude of the spot on the stellar surface, described in the same coordinate system as the star itself.
• SPOT LONGITUDE ${\phi }_{\mathrm{sp}}$: the longitude of the spot on the stellar surface, described in the same coordinate system as the star itself.
• SPOT ANGULAR RADIUS ${\varrho }_{\mathrm{sp}}$: angular size of the spot, usually given in radians.
• SPOT TEMPERATURE FACTOR ${\tau }_{\mathrm{sp}}$: the ratio between the temperature of the spot and the local temperature of the underlying photosphere. Values ${\tau }_{\mathrm{sp}}\gt 1$ correspond to hot spots, while values ${\tau }_{\mathrm{sp}}\lt 1$ correspond to cool spots.

The key assumption when modeling the area covered by a spot is that its local temperature is modified by the temperature factor and the area then radiates as if it were a normal photosphere at that modified temperature.

We can justly question the adequacy of this simplistic description: to properly model spots, advanced and reliable methods of data acquisition (polarimetry, Doppler tomography) are required (see chapter 2). It is not straightforward to use the data in the sophisticated treatment of spots (e.g., nonuniform spot brightness, magnetohydrodynamic effects, etc.), especially as introducing spots without clear support by the data may easily reproduce any physical effect—the degeneracy is typically severe. Even such a simplified treatment, along with any temporal evolution of spot parameters, is thus still ahead of the majority of our current observational capabilities.

The amount of light that reaches the observer is affected by the kinematic properties of the radiating body. There are three important velocity-induced changes to the received flux: (1) Doppler shift of the spectrum causes the passband-weighted integral to change (Equations (4.27) and (4.28)); (2) Doppler shift of the frequency causes the arrival rate of the photons to change because of time dilation; and (3) relativistic beaming, where radiation is no longer isotropic but becomes direction-dependent because of light aberration (Rybicki & Lightman 1979, chap. 4.8). Loeb & Gaudi (2003) showed that the boosting signal in exoplanet-hosting stars is larger than variability due to reflected light from the planet; Zucker et al. (2007) performed a similar study for binary stars and showed that ellipsoidal variability, reflection, and boosting are same-order-of-magnitude effects in binaries with a $\sim \,10$-day orbital period, while boosting dominates in binaries with a $\sim \,100$-day orbital period. Studies by van Kerkwijk et al. (2010) and Bloemen et al. (2011, 2012) apply this correction to Kepler objects KOI-74, KOI-81, and KPD 1946+4340 and demonstrate that the contributions of boosting are indeed crucial to model the light curves correctly.

The combined Doppler boosting signal can be written as:

$$\begin{eqnarray}&&{I}_{\lambda }={I}_{0,\lambda }\left(1-{ \mathcal B }(\lambda )\frac{{v}_{r}}{c}\right),\end{eqnarray} \tag{ 5.17 }$$

where I_λ is the boosted passband intensity, I_{0, λ} is the initial passband intensity, v_r is the radial velocity, c is the speed of light, and ${ \mathcal B }(\lambda )$ is the boosting index,

$$\begin{eqnarray}&&{ \mathcal B }(\lambda )=5+\alpha \equiv 5+\frac{d\,\mathrm{ln}\,I}{d\,\mathrm{ln}\,\lambda },\end{eqnarray} \tag{ 5.18 }$$

where $\alpha \equiv d(\mathrm{ln}\,I)/d(\mathrm{ln}\,\lambda )$ is the spectral index, and 5 comes from the Lorenz invariance of Iλ⁵ (Loeb & Gaudi 2003). Figure 5.6 depicts the boosting index for the broad Johnson V passband, for a series of ${T}_{\mathrm{eff}}=6000$ K, $\mathrm{log}\,g=4.0$, [M/H] = 0.0 spectra where μ ranges from 0.001 to 1.0. The dependence of the boosting index on wavelength is notable, but in the simplified approach we approximate the monochromatic boosting factor ${ \mathcal B }(\lambda )$ with the passband-averaged value:

$$\begin{eqnarray}&&{{ \mathcal B }}_{\mathrm{pb}}=\frac{{\int }_{\lambda }{ \mathcal P }(\lambda ){ \mathcal S }(\lambda ){ \mathcal B }(\lambda )d\lambda }{{\int }_{\lambda }{ \mathcal P }(\lambda ){ \mathcal S }(\lambda )d\lambda }\end{eqnarray} \tag{ 5.19 }$$

(or its photon-weighted counterpart per Equation (4.28)). This is done by fitting a Legendre polynomial of the fifth order to the $\mathrm{ln}\,{I}_{\mu }(\mathrm{ln}\,\lambda )$ data; the order is determined by evaluating the rank of the coefficient matrix in the least-squares fit and seeing where the values become susceptible to numerical noise. The series was fit iteratively, as data points that lie below the −1σ threshold were discarded by sigma-clipping the data set. The iteration stopped when no further points are removed. Another benefit of Legendre polynomials is their analytical derivative; this analytic form is used to derive the average boosting index for every model atmosphere. Figure 5.6 also shows that the dependence of the boosting index on μ is significant, so the traditional approach of using integrated flux SEDs to estimate ${{ \mathcal B }}_{\mathrm{pb}}$ is less precise than the treatment presented here.

Zoom In Zoom Out Reset image size

Boosting is a significant effect, and it becomes a dominant effect over ellipsoidal variation and reflection for longer-period systems. Figure 5.7 shows a comparison of the amplitude of ellipsoidal variation (blue), reflection (green), and Doppler boosting (red) as a function of orbital period for an A0–K0 main-sequence pair in the Johnson V passband (Prša et al. 2016). Ellipsoidal variation and reflection dominate the short-period end, with boosting taking over at around 8 days. Figure 5.8 depicts a light curve with no effects computed (black), with only ellipsoidal variations (blue), with ellipsoidal variations and reflection (green), and with ellipsoidal variations, reflection, and boosting (red).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Interstellar extinction is an extrinsic effect that causes a wavelength-dependent reduction of light that reaches the observer. It was first reported by Trumpler (1930) and connected with dust grains that form the interstellar medium. Extinction is also known as reddening because it affects the blue part of the spectrum more than the red, making objects appear redder than they are intrinsically (Draine 2003). It is usually described by the color excess:

$$\begin{eqnarray}&&E{(B-V)={(B-V)}_{\mathrm{observed}}-(B-V)}_{\mathrm{intrinsic}}.\end{eqnarray} \tag{ 5.20 }$$

It provides a measure of reddening in terms of the B − V color index change. The slope of the extinction curve in the optical part of the spectrum is given by:

$$\begin{eqnarray}&&{R}_{V}\equiv \frac{{A}_{V}}{{A}_{B}-{A}_{V}},\end{eqnarray} \tag{ 5.21 }$$

where A_B and A_V are extinctions (in mag) in the B and V passbands, respectively. The value of R_V depends on the grain size. For example, Rayleigh scattering similar to that in our atmosphere ($A(\lambda )\propto {\lambda }^{-4}$) would produce a slope of R_V ∼ 1.2; extinction in our Galaxy ranges between ∼2 and ∼6, with an average around R_V ∼ 3.1. The larger the grain size, the larger the R_V ; in the limit of very large grain sizes, ${R}_{V}\to \infty$ (Draine 2003). This regime is called gray extinction, i.e., extinction does not depend on the wavelength.

The most frequently used model to describe wavelength dependence of interstellar extinction is that of Cardelli et al. (1989). It is an empirical formula with the assumed extinction slope R_V = 3.1. The formula is split into three wavelength regions:

Region	λ-interval	$x=(1/\lambda )$-interval
UV	$170\,\mathrm{nm}\leqslant \lambda \lt 303\,\mathrm{nm}$	$3.3\,\mu {{\rm{m}}}^{-1}\leqslant x\lt 5.9\,\mu {{\rm{m}}}^{-1}$
Optical/NIR	$303\,\mathrm{nm}\leqslant \lambda \lt 910\,\mathrm{nm}$	$1.1\,\mu {{\rm{m}}}^{-1}\leqslant x\lt 3.3\,\mu {{\rm{m}}}^{-1}$
IR	$910\,\mathrm{nm}\leqslant \lambda \lt 3330\,\mathrm{nm}$	$0.3\,\mu {{\rm{m}}}^{-1}\leqslant x\lt 1.1\,\mu {{\rm{m}}}^{-1}$

The amount of interstellar extinction at a given wavelength (see Figure 5.9) is described by the interpolation formula:

$$\begin{eqnarray}&&\left\langle \frac{A(\lambda )}{{A}_{V}}\right\rangle =a(x)+b(x)/{{ \mathcal R }}_{{\rm{V}}},\end{eqnarray} \tag{ 5.22 }$$

where A(λ) is extinction in mag at the given wavelength, expressed relative to the visual extinction A_V ; x = 1/λ is given in $\mu {{\rm{m}}}^{-1}$. Coefficients a(x) and b(x) are computed by the following expressions:

Region
UV	$a{(x)=1.752-0.316x-0.104/[(x-4.67)}^{2}+0.341]-0.04473{(x-5.9)}^{2}\,-\,0.009779{(x-5.9)}^{3}$
	$b{(x)=-3.090+1.825x+1.206/[(x-4.62)}^{2}+0.263]+0.2130{(x-5.9)}^{2}\,+\,0.1207{(x-5.9)}^{3}$
Optical/NIR	$\,y:= (x-1.82)$
	$a(x)=1+0.17699y-0.50447{y}^{2}-0.02427{y}^{3}+0.72085{y}^{4}+0.01979{y}^{5}-0.77530{y}^{6}+0.32999{y}^{7}$
	$b(x)=1.41338y+2.28305{y}^{2}+1.07233{y}^{3}-5.38434{y}^{4}-0.62251{y}^{5}+5.30260{y}^{6}-2.09002{y}^{7}$
IR	$a(x)=0.574{x}^{1.61}$
	$b(x)=-0.527{x}^{1.61}$

Zoom In Zoom Out Reset image size

A standard approach used in the literature to deredden the data is to estimate the amount of reddening from extinction tables (e.g., Schlegel et al. 1998; Gonzalez et al. 2013; Green et al. 2018) given an object's equatorial coordinates. The tables provide color excess along the line of sight, which, given R_V and distance, can be used to calculate passband extinction A(λ). This extinction is then subtracted from photometric observations. These tables are either 2D (integrated along the line of sight) or 3D (provide color excess as a function of distance). The "ultimate" extinction table is expected to come from ESA's Gaia mission (Schultheis et al. 2015).

There is a significant problem with this standard approach when applied to eclipsing binary stars: as there are two stars that contribute to total flux and their relative contributions change with phase (i.e., ellipsoidal variations, eclipses, spots, etc.), the amount of reddening also changes, and subtracting a single constant introduces systematic trends as a function of phase. The magnitude of these trends increases with the color difference between the components and with the amount of extinction (Prša & Zwitter 2005).

A much more suitable approach to account for interstellar extinction is to multiply the underlying model atmospheres by monochromatic extinction A(λ) and integrate over wavelength:

$$\begin{eqnarray}&&\langle {I}_{\mathrm{pb},\mathrm{norm}}\rangle =\frac{{\int }_{\lambda }{ \mathcal P }(\lambda )A(\lambda ){I}_{\lambda ,\mathrm{norm}}\,(\lambda )\,d\lambda }{{\int }_{\lambda }{ \mathcal P }(\lambda )A(\lambda )\,d\lambda },\end{eqnarray} \tag{ 5.23 }$$

so Equation (4.27) becomes:

$$\begin{eqnarray}&&{F}_{\mathrm{pb}}=\frac{\mu ^{\prime} }{{D}^{2}}\left({\int }_{\lambda }{ \mathcal P }(\lambda )A(\lambda )d\lambda \right){\int }_{\partial V}\langle {I}_{\mathrm{pb},\mathrm{norm}}\rangle \mu { \mathcal V }({\boldsymbol{r}}){ \mathcal L }(\mu ){dA}.\end{eqnarray} \tag{ 5.24 }$$

Thus, we can expand the atmosphere lookup tables by adding another dimension, A_pb, and interpolate as before to obtain the correct amount of interstellar extinction added to the model (Jones et al. 2018, in preparation).

To demonstrate the effect of reddening on eclipsing binaries, consider a pair of main-sequence stars with temperatures 6500 and 5500 K, orbiting in a circular orbit with an i = 90° inclination. Say we observe such a binary in Johnson B and V passbands. We will correct for extinction in two ways: by using the standard approach and by using the wavelength-dependent approach. We will neglect all other effects (limb darkening, gravity brightening, reflection effect) so that we evaluate solely the impact of reddening.

The logic is as follows: for the given phase, we compute the effective spectrum of the binary by summing the Doppler-shifted spectra of the visible surfaces of both components. To this intrinsic spectrum we apply interstellar extinction as a function of wavelength. We then multiply this reddened spectrum by the instrumental response function (filter transmission × detector response) and integrate over the passband wavelength range to obtain the flux.

For the standard approach we use the same intrinsic spectrum without applying the wavelength-dependent extinction; instead, we divide the intrinsic spectrum by the extinction constant and integrate as before. We can then readily compare the two fluxes as a function of phase. There is a catch, though: as the Cardelli et al. (1989) formula depends on the wavelength, the extinction constant will also depend on the wavelength. So which wavelength do we choose? We could naively take the effective wavelength of the passband, i.e., the wavelength that splits the passband integral into halves. That is not what we need, though; we need a wavelength that splits the product ${ \mathcal P }(\lambda )A(\lambda ){I}_{\lambda }({T}_{\mathrm{eff}},\mathrm{log}\,g,\ldots )$ into halves. To compute this, we need to repeat the same wavelength integration process as before just to get a suitable effective wavelength. Table 5.1 summarizes effective wavelengths computed as a geometric mean of different quantities: passband, reddened passband, intrinsic spectrum, reddened spectrum, and the "true" composite reddened spectrum. Clearly, effective wavelength depends strongly on the effective temperature of the observed object and on the color excess $E(B-V)$ (Figure 5.10). Figure 5.11 depicts the discrepancy between the properly calculated light curve and the one obtained by subtracting a ${\lambda }_{\mathrm{eff}}$-calculated constant.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 5.1.  Summary of Different Approaches to Calculate the Wavelength to Be Used for the Dereddening Constant.

Approach	${\lambda }_{\mathrm{eff}}$	${\rm{\Delta }}{m}_{B}$	${\epsilon }_{{\rm{\Delta }}{m}_{B}}$
Rigorously calculated value	4470.2 Å	4.03	0.00
Passband transmission	4410.8 Å	4.09	0.06
Passband transmission + reddening law	4452.1 Å	4.05	0.02
Intrinsic spectrum	4436.3 Å	4.06	0.03
Effective (reddened) spectrum	4583.6 Å	3.90	−0.13

${\lambda }_{\mathrm{eff}}$ is determined by requiring that the area under the spectrum on both sides is equal. ${\rm{\Delta }}{m}_{B}$ is the value of extinction in the B passband, and ${\epsilon }_{{\rm{\Delta }}{m}_{B}}$ is the deviation from the rigorously calculated value. All values are calculated for $E(B-V)=1$ at quarter phase.

Even if the subtraction constant is properly calculated, the light curves will still exhibit measurable differences in both minima (Figure 5.11). This is because of the color index (effective temperature) dependence on phase, most notably during eclipses. For the case at hand, the difference in B magnitude is ∼1.75 mmag. For the case of a B0V–K4III binary, the same difference inflates to ∼0.2 mag (Prša & Zwitter 2005).

If light curves in three or more photometric passbands are available, it is possible to uniquely determine the color excess value $E(B-V)$ by comparing color indices in and out of eclipse. The reddening may thus be properly introduced as a parameter that can be adjusted during the fitting process.

Atmospheric extinction is an equivalent but a better-posed problem. As interstellar extinction depends on $E(B-V)$, atmospheric extinction depends on air mass, which (unlike $E(B-V)$) is a directly measurable quantity. The downside of atmospheric extinction compared to interstellar extinction is its short-scale time variability: due to the rapidly changing nature of atmospheric convection cells, the presence of cirrus clouds, and water vapor, air mass changes unpredictably during the course of observations.

To model atmospheric extinction, we use the equation triplet for Rayleigh-ozone-aerosol extinction sources given by Forbes et al. (1996) and summarized by Pakštiene & Solheim (2003):

$$\begin{eqnarray}&&{k}_{\mathrm{total}}={k}_{\mathrm{Rayleigh}}+{k}_{\mathrm{ozone}}+{k}_{\mathrm{aerosol}}.\end{eqnarray} \tag{ 5.25 }$$

Rayleigh scattering is given by the following expression:

$$\begin{eqnarray}&&{k}_{\mathrm{Rayleigh}}={k}_{\mathrm{Rayleigh}}^{\mathrm{STP}}\left(\frac{P}{{P}_{\mathrm{STP}}}\right)\left(\frac{{T}_{\mathrm{STP}}}{T}\right),\end{eqnarray} \tag{ 5.26 }$$

where P and T are the pressure and temperature during observations; P_STP = 101,300 Pa and ${T}_{\mathrm{STP}}=273.15\,{\rm{K}}$ are the standard temperature and pressure, respectively. Rayleigh scattering for these standard values is given by:

$$\begin{eqnarray}&&{k}_{\mathrm{Rayleigh}}^{\mathrm{STP}}=9.4977\cdot {10}^{-3}{e}^{-h/7.996}\frac{1}{{\lambda }^{4}}{\left(\frac{n(\lambda )-1}{n(\lambda =1\mu {\rm{m}})-1}\right)}^{2},\end{eqnarray} \tag{ 5.27 }$$

where λ is the wavelength in $\mu {\rm{m}}$. The formula is normalized to $\lambda =1\mu {\rm{m}}$, with n being the refractive index and the coefficient 9.4977 × 10⁻³ corresponding to the optical thickness of the atmosphere at $\lambda =1\,\mu {\rm{m}}$. This is adjusted for the height h (in km) of the observatory above sea level, with the density scale height for the lower troposphere assumed to be 7.996 km (Forbes et al. 1996). The normalized index of refraction term may be expressed in terms of wavelength by an empirical relation proposed by Penndorf (1957):

$$\begin{eqnarray}&&\left(\frac{n(\lambda )-1}{n(\lambda =1\mu {\rm{m}})-1}\right)=0.23465+\frac{107.6}{146-{\lambda }^{-2}}+\frac{0.93161}{41-{\lambda }^{-2}}.\end{eqnarray} \tag{ 5.28 }$$

Ozone extinction depends on the amount α_ozone of ozone in the vertical column of the atmosphere with cross section 1 cm² and the ozone absorption coefficient ${\kappa }_{\mathrm{ozone}}(\lambda )$:

$$\begin{eqnarray}&&{k}_{\mathrm{ozone}}=1.1\,{\alpha }_{\mathrm{ozone}}\,{\kappa }_{\mathrm{ozone}}(\lambda ),\end{eqnarray} \tag{ 5.29 }$$

where 1.1 corresponds to the ozone optical thickness at 320 nm. Values of ${\kappa }_{\mathrm{ozone}}$ are tabulated, e.g., in Cox (2000).

Aerosol scattering on particles with size of the same order of magnitude as the wavelength is the most complicated contribution to the overall atmospheric extinction. It is modeled by the power law:

$$\begin{eqnarray}&&{k}_{\mathrm{aerosol}}=A{\lambda }^{-b},\end{eqnarray} \tag{ 5.30 }$$

where A is the extinction at zenith for $\lambda =1\,\mu {\rm{m}}$ and b is the coefficient that depends on the size and distribution of aerosol particles. These two parameters are very difficult to assess; fortunately, the aerosol contribution is usually very small, even negligible (Pakštiene & Solheim 2003). We neglect its influence in our simulation.

In terms of binary stars, atmospheric extinction depends on the wavelength of the observed light and on the air mass of observations, which is of course the same for both binary components; we are thus left only with wavelength dependence at some inferred air mass. We can therefore compare the flux change imposed on the intrinsic spectrum by reddening and by atmospheric extinction (Figure 5.12). For weakly and moderately reddened eclipsing binaries atmospheric extinction dominates the blue parts of the spectrum owing to Rayleigh scattering, but for larger color excesses ($E(B-V)\gtrsim 0.5$) interstellar extinction is the dominant factor throughout the spectrum. Since air mass can be measured reliably from observations of comparison stars, the effect of atmospheric extinction can in principle be removed from the data.

Zoom In Zoom Out Reset image size