Bolometric Treatment of Irradiation Effects: General Discussion and Application to Binary Stars

We are primarily focusing on radiating bodies with a simple geometry, especially close to spherical, where we can describe the redistribution processes and identify the surface parts associated with flux incidence and emission. For strongly deformed bodies, such as contact binaries, we currently lack sufficient insight to propose a realistic yet tractable model because, as of yet, the thermodynamical properties are not well understood. Substantial effort to better understand radiative transfer in strongly deformed bodies is underway (Kochoska et al. 2018).

In the stationary or quasi-stationary state, we assume that the flux is strictly conserved at all times, meaning that the net incident flux is also emitted at the same time. The absorbed part of the irradiance, (1 − ρ)F_in, is redistributed over the entire surface, and it increases the intrinsic exitance by δF₀. We can describe flux redistribution by defining a redistribution operator $\hat{{ \mathcal D }}$,

$$\begin{eqnarray}&&\delta {F}_{0}=\hat{{ \mathcal D }}\left[(1-\rho ){F}_{\mathrm{in}}\right]=\hat{{ \mathcal D }}({\mathbb{1}}-\hat{{\rm{\Pi }}}){F}_{\mathrm{in}}.\end{eqnarray} \tag{ 17 }$$

This redistribution operator could in principle be time-dependent in the quasi-stationary case, but for the strictly stationary state redistribution, the operator does not depend on time. The increased intrinsic exitance ${F}_{0}^{{\prime} }$ can be written as

$$\begin{eqnarray}&&{F}_{0}^{{\prime} }={F}_{0}+\delta {F}_{0}\quad {F}_{\mathrm{out}}={F}_{0}^{{\prime} }+\hat{{\rm{\Pi }}}{F}_{\mathrm{in}},\end{eqnarray} \tag{ 18 }$$

where we work under the assumption that the intrinsic exitance F₀ is not affected by irradiation. By construction, the redistribution operator $\hat{{ \mathcal D }}$ maps positive-valued functions defined over the surface to positive-valued functions and conserves their integrals over the surface.

Generally, we can write the redistribution operator $\hat{{ \mathcal D }}$ at time t as

$$\begin{eqnarray*}&&\hat{{ \mathcal D }}f({\boldsymbol{r}})={\int }_{{ \mathcal M }}K({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} ;t)f({\boldsymbol{r}}^{\prime} )\,{dA}({\boldsymbol{r}}^{\prime} ),\end{eqnarray*}$$

where $K:{ \mathcal M }\times { \mathcal M }\to {{\mathbb{R}}}_{+}$ is a positive kernel with the following normalization property:

$$\begin{eqnarray*}&&{\int }_{{ \mathcal M }}K({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} ;t)\,{dA}({\boldsymbol{r}})=1\quad \forall {\boldsymbol{r}}\in { \mathcal M }.\end{eqnarray*}$$

This imposes the conservation of total flux at a given moment in time over the surface:

$$\begin{eqnarray}&&{\int }_{{ \mathcal M }}[\hat{{ \mathcal D }}f]({\boldsymbol{r}})\,{dA}({\boldsymbol{r}})={\int }_{{ \mathcal M }}f({\boldsymbol{r}})\,{dA}({\boldsymbol{r}}).\end{eqnarray} \tag{ 19 }$$

More generally, we can formulate the kernel by using an auxiliary function $G:{ \mathcal M }\times { \mathcal M }\to {{\mathbb{R}}}_{+}$:

$$\begin{eqnarray*}&&K({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} ;t)=\displaystyle \frac{G({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} ;t)}{{\int }_{{ \mathcal M }}G({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} ;t){dA}({\boldsymbol{r}})}.\end{eqnarray*}$$

When the flux is uniformly distributed over the whole surface, G_uniform ≡ 1 and the redistribution operator $\hat{{{ \mathcal D }}_{\mathrm{uni}}}$ is very simple:

$$\begin{eqnarray*}&&{\hat{{ \mathcal D }}}_{\mathrm{uni}}f({\boldsymbol{r}})=\displaystyle \frac{1}{A}{\int }_{{ \mathcal M }}f({\boldsymbol{r}}^{\prime} ){dA}({\boldsymbol{r}}^{\prime} ).\end{eqnarray*}$$

We model a local redistribution by introducing a distance measure (e.g., a geodesic) on the surface, $d:{ \mathcal M }\times { \mathcal M }\to {{\mathbb{R}}}_{+}$, and a weight function, $g:{{\mathbb{R}}}_{+}\to {{\mathbb{R}}}_{+}$, that determines the ratio of the flux that is transported from the irradiated element to any other element on the surface at distance d. Then, the kernel describing the redistribution of incident flux is written as

$$\begin{eqnarray}&&{G}_{\mathrm{loc}}({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} ;t)=g\left(d({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} )\right).\end{eqnarray} \tag{ 20 }$$

The redistribution operator associated with this kernel is denoted by ${\hat{{ \mathcal D }}}_{\mathrm{loc}}$. The weight function g depends on the energy transport in the atmosphere, which we do not (readily) know. However, we can make reasonable assumptions that are likely to hold. Namely, we take g to be monotonically decreasing and diminishing to zero for arguments larger than a given threshold value l, for example, $g(x)=\exp (-x/l)$ or g(x) = 1−x/l, where l is proportional to the optical depth in the atmosphere. In the limit l → 0, local redistribution reflects all incoming flux: ${\hat{{ \mathcal D }}}_{\mathrm{loc}}={\mathbb{1}}$ at l = 0. Therefore, local redistribution with small l (w.r.t. to the size of object) will have an effect similar to increasing reflection. For spherical bodies of radius R, it is much more convenient to use the ratio l/R as a parameter determining the threshold value.

In the case of rotating stars, the axis of rotation breaks the isotropic symmetry, so excess flux tends to be reradiated at latitudes similar to where it was received. This gives rise to flux conservation in the latitudinal direction. It is meaningful to define a distance on the surface, ${d}_{\perp }:{ \mathcal M }\times { \mathcal M }\to {{\mathbb{R}}}_{+}$, along the rotation axis $\hat{{\boldsymbol{s}}}$. The kernel can then be written as

$$\begin{eqnarray}&&{G}_{\mathrm{lat}}({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} ;t)=g\left({d}_{\perp }({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} ;\hat{{\boldsymbol{s}}})\right).\end{eqnarray} \tag{ 21 }$$

and the corresponding redistribution operator is labeled as ${\hat{{ \mathcal D }}}_{\mathrm{lat}}$.

Let us make a small digression here and note that, in the case of a phase delay between heating and reradiation, it is possible to incorporate the lag into the presented framework by using coordinates shifted horizontally w.r.t. the rotation axis $\hat{{\boldsymbol{s}}}$. The kernel for that case would be written as

$$\begin{eqnarray*}&&{G}_{\mathrm{shift}}({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} ;t)=G({{\boldsymbol{R}}}_{\phi \hat{{\boldsymbol{s}}}}{\boldsymbol{r}},{\boldsymbol{r}}^{\prime} ;t),\end{eqnarray*}$$

where ϕ is the angle by which the location of irradiated element is shifted relative to the location of the radiating element, and ${{\boldsymbol{R}}}_{{\boldsymbol{\omega }}}$ is the rotation matrix about the axis of rotation ω. We expect that the angle ϕ is positive and proportional to the angular velocity of the star; in general it can also depend on the position and time, as long as the quasi-stationarity of redistribution assumed here is not violated. Note that the horizontal shift does not affect the overall energy balance.

We can describe individual redistribution processes by the corresponding operator ${\hat{{ \mathcal D }}}_{i}$ and form the overall redistribution operator $\hat{{ \mathcal D }}$ as a weighted sum of individual ${\hat{{ \mathcal D }}}_{i}$:

$$\begin{eqnarray}&&\hat{{ \mathcal D }}=\displaystyle \sum _{i}{w}_{i}{\hat{{ \mathcal D }}}_{i},\quad \displaystyle \sum _{i}{w}_{i}=1,\end{eqnarray} \tag{ 22 }$$

where w_i are positive real numbers. We can view Equation (22) as the decomposition of the redistribution operator $\hat{{ \mathcal D }}$ into generators of a certain type of redistribution, where the weights quantify the amount of energy redistributed by the corresponding process.

The radiosity for Wilson's model can be obtained by substituting ${F}_{0}^{{\prime} }$ from Equation (18) into Equations (14):

$$\begin{eqnarray}&&{F}_{\mathrm{out}}={F}_{0}+\left[\hat{{ \mathcal D }}({\mathbb{1}}-\hat{{\rm{\Pi }}})+\hat{{\rm{\Pi }}}\right]{\hat{{ \mathcal L }}}_{\mathrm{LD}}{F}_{\mathrm{out}}.\end{eqnarray} \tag{ 23 }$$

In turn, the updated intrinsic exitance ${F}_{0}^{{\prime} }$ is given by the radiosity:

$$\begin{eqnarray}&&{F}_{0}^{{\prime} }={F}_{0}+\hat{{ \mathcal D }}({\mathbb{1}}-\hat{{\rm{\Pi }}}){\hat{{ \mathcal L }}}_{\mathrm{LD}}{F}_{\mathrm{out}}.\end{eqnarray} \tag{ 24 }$$

Similarly, the irradiation for the Lambertian reflection model can be written as

$$\begin{eqnarray}&&{F}_{\mathrm{in}}={\hat{{ \mathcal L }}}_{\mathrm{LD}}{F}_{0}+\left[{\hat{{ \mathcal L }}}_{\mathrm{LD}}\hat{{ \mathcal D }}({\mathbb{1}}-\hat{{\rm{\Pi }}})+{\hat{{ \mathcal L }}}_{{\rm{L}}}\hat{{\rm{\Pi }}}\right]{F}_{\mathrm{in}}.\end{eqnarray} \tag{ 25 }$$

Equations (23) and (25) are the main theoretical results of the paper and represent a unification of a specific reflection scheme and irradiation redistribution under one irradiation framework. The solution of these equations determines the intrinsic exitance ${F}_{0}^{{\prime} }$ and radiosity F_out:

$$\begin{eqnarray}\begin{array}{rcl}{F}_{0}^{{\prime} } & = & {F}_{0}+\hat{{ \mathcal D }}({\mathbb{1}}-\hat{{\rm{\Pi }}}){F}_{\mathrm{in}},\\ {F}_{\mathrm{out}} & = & {F}_{0}^{{\prime} }+\hat{{\rm{\Pi }}}{F}_{\mathrm{in}}.\end{array}\end{eqnarray} \tag{ 26 }$$

The solutions of the irradiation models, i.e., the updated exitance ${F}_{0}^{{\prime} }$ and radiosity F_out, determine the bolometric intensity of the radiating bodies. For Wilson's model, the limb-darkened intensity from Equation (8) yields

$$\begin{eqnarray}&&I(\hat{{\boldsymbol{e}}},{\boldsymbol{r}})=\displaystyle \frac{{F}_{\mathrm{out}}({\boldsymbol{r}})}{{D}_{0}({\boldsymbol{r}})}D(\hat{{\boldsymbol{e}}},{\boldsymbol{r}}),\end{eqnarray} \tag{ 27 }$$

while for the Lambertian reflection model we get

$$\begin{eqnarray}&&I(\hat{{\boldsymbol{e}}},{\boldsymbol{r}})=\displaystyle \frac{{F}_{0}^{{\prime} }({\boldsymbol{r}})}{{D}_{0}({\boldsymbol{r}})}D(\hat{{\boldsymbol{e}}},{\boldsymbol{r}})+\displaystyle \frac{1}{\pi }\left({F}_{\mathrm{out}}({\boldsymbol{r}})-{F}_{0}^{{\prime} }({\boldsymbol{r}})\right).\end{eqnarray} \tag{ 28 }$$

The presented Lambertian irradiation (reflection and redistribution) model is an exact bolometric description of this process under three assumptions: (1) Lambertian reflection is wavelength independent, (2) redistribution does not affect the limb-darkening of the surface, and (3) intrinsic exitance is not affected by the irradiation.

To quantify the impact of Lambertian correction, we can expand the updated intrinsic emission ${F}_{0}^{{\prime} }$ and radiosity F_out. For Wilson's model, we get

$$\begin{eqnarray}\begin{array}{rcl}{F}_{\mathrm{out}} & = & {F}_{0}+\left[\hat{{\rm{\Pi }}}+\hat{{ \mathcal D }}({\mathbb{1}}-\hat{{\rm{\Pi }}})\right]{\hat{{ \mathcal L }}}_{\mathrm{LD}}{F}_{0}\\ & & +\ \underline{\left[\hat{{\rm{\Pi }}}+\hat{{ \mathcal D }}({\mathbb{1}}-\hat{{\rm{\Pi }}})\right]{\hat{{ \mathcal L }}}_{\mathrm{LD}}\left[\hat{{\rm{\Pi }}}+\hat{{ \mathcal D }}({\mathbb{1}}-\hat{{\rm{\Pi }}})\right]{\hat{{ \mathcal L }}}_{\mathrm{LD}}{F}_{0}}+\ldots ,\end{array}\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}\begin{array}{rcl}{F}_{0}^{{\prime} } & = & {F}_{0}+\hat{{ \mathcal D }}({\mathbb{1}}-\hat{{\rm{\Pi }}}){\hat{{ \mathcal L }}}_{\mathrm{LD}}{F}_{0}\\ & & +\ \underline{\hat{{ \mathcal D }}({\mathbb{1}}-\hat{{\rm{\Pi }}}){\hat{{ \mathcal L }}}_{\mathrm{LD}}\left[\hat{{\rm{\Pi }}}+\hat{{ \mathcal D }}({\mathbb{1}}-\hat{{\rm{\Pi }}})\right]{\hat{{ \mathcal L }}}_{\mathrm{LD}}{F}_{0}}+\ldots ,\end{array}\end{eqnarray} \tag{ 30 }$$

whereas in the redistribution models based on the Lambertian reflection, the expansions are written as

$$\begin{eqnarray}&&\begin{array}{l}{F}_{\mathrm{out}}={F}_{0}+\left[\hat{{\rm{\Pi }}}+\hat{{ \mathcal D }}({\mathbb{1}}-\hat{{\rm{\Pi }}})\right]{\hat{{ \mathcal L }}}_{\mathrm{LD}}{F}_{0}\\ \quad +\ \underline{\left[\hat{{\rm{\Pi }}}+\hat{{ \mathcal D }}({\mathbb{1}}-\hat{{\rm{\Pi }}})\right]\left[{\hat{{ \mathcal L }}}_{\mathrm{LD}}\hat{{ \mathcal D }}({\mathbb{1}}-\hat{{\rm{\Pi }}})+{\hat{{ \mathcal L }}}_{0}\hat{{\rm{\Pi }}}\right]{\hat{{ \mathcal L }}}_{\mathrm{LD}}{F}_{0}}+\ldots ,\end{array}\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}\begin{array}{rcl}{F}_{0}^{{\prime} } & = & {F}_{0}+\hat{{ \mathcal D }}({\mathbb{1}}-\hat{{\rm{\Pi }}}){\hat{{ \mathcal L }}}_{\mathrm{LD}}{F}_{0}\\ & & +\ \underline{\hat{{ \mathcal D }}({\mathbb{1}}-\hat{{\rm{\Pi }}})\left[{\hat{{ \mathcal L }}}_{\mathrm{LD}}\hat{{ \mathcal D }}({\mathbb{1}}-\hat{{\rm{\Pi }}})+{{ \mathcal L }}_{0}\hat{{\rm{\Pi }}}\right]{\hat{{ \mathcal L }}}_{\mathrm{LD}}{F}_{0}}+\ldots .\end{array}\end{eqnarray} \tag{ 32 }$$

By comparing the expressions for F_out (${F}_{0}^{{\prime} }$) in different reflection approaches we see that they differ in the underlined second-order terms. The difference is typically small and likely not measurable at the current level of precision. However, it is conceptually important as it corresponds to a different physical description of the surface boundary energy balance.

Energy is conserved when the difference between the total emitted flux L_out and the total incident flux L_in equals the total intrinsic flux L₀:

$$\begin{eqnarray}&&{L}_{\mathrm{out}}-{L}_{\mathrm{in}}={L}_{0}\quad {L}_{* }={\displaystyle \int }_{{ \mathcal M }}{F}_{* }({\boldsymbol{r}}){dA}({\boldsymbol{r}}),\end{eqnarray} \tag{ 33 }$$

where * = out, in, 0.

If we want to account for the processes that are not included in the energy balance (such as scattering), then the total flux per Equation (19) is not conserved. The losses can occur at different levels:

• (a)  
  the decrease in the non-reflected part of the incident light at the surface of the irradiated star, described by the scalar function ξ(r) ∈ [0, 1]:

  $$\begin{eqnarray}&&\delta {F}_{0}=\hat{{ \mathcal D }}(1-\hat{{\rm{\Pi }}})\hat{{\rm{\Xi }}}{F}_{\mathrm{in}},\end{eqnarray} \tag{ 34 }$$

  where we introduce the auxiliary operator $\hat{{\rm{\Xi }}}:f\mapsto \xi f$ to describe the losses;
• (b)  
  the decrease of energy in the interior of the irradiated star, described by a "lossy" redistribution operator $\hat{{ \mathcal D }}^{\prime}$:

  $$\begin{eqnarray}&&\delta {F}_{0}=\hat{{ \mathcal D }}^{\prime} ({\mathbb{1}}-\hat{{\rm{\Pi }}}){F}_{\mathrm{in}},\end{eqnarray} \tag{ 35 }$$

  where the redistribution operator has the following property for an arbitrary positive function F:

  $$\begin{eqnarray}&&{\int }_{{ \mathcal M }}[\hat{{ \mathcal D }}^{\prime} F]({\boldsymbol{r}}){dA}({\boldsymbol{r}})\leqslant {\int }_{{ \mathcal M }}F({\boldsymbol{r}}){dA}({\boldsymbol{r}})\,;\quad \mathrm{or}\end{eqnarray} \tag{ 36 }$$
• (c)  
  the decrease in the emergent light at the surface of the radiating star, described analogously to case (a):

  $$\begin{eqnarray}&&\delta {F}_{0}=\hat{{\rm{\Xi }}}\hat{{ \mathcal D }}({\mathbb{1}}-\hat{{\rm{\Pi }}}){F}_{\mathrm{in}}.\end{eqnarray} \tag{ 37 }$$

  which is just a specific case of the previous with $\hat{{ \mathcal D }}^{\prime} =\hat{{\rm{\Xi }}}\hat{{ \mathcal D }}$, from the modeling point of view.

Case (a) could be used to describe losses due to scattering from the surface that are not taken into account by the reflection; case (b) could mimic the absorption of energy by processes inside of the star that violate flux conservation in a semi-stationary regime, e.g., altering the internal dynamics of the envelope in convective stars (Ruciński 1969); and case (c) could model obstructions for the re-emission of redistributed irradiation (e.g., some irregularities on the surface of the stars in the form of spots, convective cells, etc.). How to translate the mentioned processes into the redistribution model is beyond the scope of this paper.

If the loss and reflection coefficients ξ and ρ are constant over the surface, we can express flux losses L_loss = L₀ − (L_out − L_in) in cases (a) and (b) as

$$\begin{eqnarray*}&&{L}_{\mathrm{loss}}={L}_{\mathrm{in}}(1-\rho )(1-\xi ),\end{eqnarray*}$$

where L_in cannot be written in a simple form as it depends on the radiosity operators. Setting ξ = 1 eliminates the losses. Flux conservation can be described by the three fractions, all with respect to the incident flux: the part of the incident flux reflected from the surface, R_refl; the part of the flux absorbed and then redistributed across the surface, R_redistr; and the part of the flux that is lost at the surface, R_lost:

$$\begin{eqnarray}&&{R}_{\mathrm{refl}}+{R}_{\mathrm{redistr}}+{R}_{\mathrm{lost}}=1,\end{eqnarray} \tag{ 38 }$$

where R_refl ≡ ρ, R_redistr ≡ ξ (1 − ρ), and R_lost ≡ (1 − ξ) (1 − ρ).

The local effective temperature of a surface element is defined by the radiosity F_out of blackbody emission:

$$\begin{eqnarray}&&{T}_{\mathrm{eff}}({\boldsymbol{r}})=\displaystyle \frac{1}{\sigma }{F}_{\mathrm{out}}{\left({\boldsymbol{r}}\right)}^{\displaystyle \frac{1}{4}},\end{eqnarray} \tag{ 39 }$$

where σ is the Stefan–Boltzmann constant. In the absence of reflection, the radiosity of a star equals the intrinsic exitance F₀ and the corresponding intrinsic effective temperature is denoted by T_eff,0, which is distributed across the surface according to the adopted gravity darkening model (i.e., von Zeipel 1924 for radiative photospheres). In the presence of reflection and redistribution, the intrinsic exitance gets updated to ${F}_{0}^{{\prime} }$ and the radiosity equals F_out. The effective temperature associated with the updated intrinsic exitance is

$$\begin{eqnarray}&&{T}_{\mathrm{eff},0}^{{\prime} }({\boldsymbol{r}})={T}_{\mathrm{eff},0}({\boldsymbol{r}}){\left(\displaystyle \frac{{F}_{0}^{{\prime} }({\boldsymbol{r}})}{{F}_{0}({\boldsymbol{r}})}\right)}^{\displaystyle \frac{1}{4}},\end{eqnarray} \tag{ 40 }$$

and the effective local temperature associated with radiosity is equal to

$$\begin{eqnarray}\begin{array}{rcl}{T}_{\mathrm{eff}}({\boldsymbol{r}}) & = & {T}_{\mathrm{eff},0}^{{\prime} }({\boldsymbol{r}}){\left(\displaystyle \frac{{F}_{\mathrm{out}}({\boldsymbol{r}})}{{F}_{0}^{{\prime} }({\boldsymbol{r}})}\right)}^{\displaystyle \frac{1}{4}}\\ & = & \ {T}_{\mathrm{eff},0}^{{\prime} }({\boldsymbol{r}}){\left(1+\rho ({\boldsymbol{r}})\displaystyle \frac{{F}_{\mathrm{in}}({\boldsymbol{r}})}{{F}_{0}^{{\prime} }({\boldsymbol{r}})}\right)}^{\displaystyle \frac{1}{4}}.\end{array}\end{eqnarray} \tag{ 41 }$$

The results of the irradiation framework presented here are bolometric quantities, which are wavelength independent. Because of that, we can only synthesize wavelength-independent (bolometric) observations. That said, the treatment described above lends itself readily to the wavelength-dependent re-emission approximation analogous to that of Wilson (1990). We outline this procedure for a given wavelength-dependent plane-parallel stellar atmospheric model, distribution of local effective temperature T_eff,0, and other properties of the atmosphere across the isolated star:

• 1.  
  The atmospheric model determines the spectral intensity on the surface of the star given by

  $$\begin{eqnarray}&&{I}_{\mathrm{atm}}(\lambda ,\cos \theta ,T,\,\ldots ),\end{eqnarray} \tag{ 42 }$$

  where λ is the wavelength, θ is the angle from the normal to the plane, T is the local effective temperature, and "..." marks all other parameters that determine properties of the atmosphere. From the spectral intensity we calculate the intrinsic exitance F₀ at each point r of the surface,

  $$\begin{eqnarray}&&\,\begin{array}{c}{F}_{0}({\boldsymbol{r}})=2\pi {\int }_{0}^{\infty }d\lambda {\int }_{0}^{1}d\mu \\ \,\mu {I}_{\mathrm{atm}}(\lambda ,\mu ,{T}_{\mathrm{eff},0}({\boldsymbol{r}}),\,\ldots ).\end{array}\end{eqnarray} \tag{ 43 }$$
• 2.  
  From the intrinsic exitance F₀, albedo ρ, and parameters determining redistribution, we calculate here the presented irradiation framework using the updated exitance ${F}_{0}^{{\prime} }$ and radiosity F_out.
• 3.  
  The radiosity F_out determines the local effective temperature T_eff (Equation (41)). Following Wilson (1990), we may use the effective temperature as the new local temperature in the spectral intensity,

  $$\begin{eqnarray}&&{I}_{\mathrm{atm}}(\lambda ,\cos \theta ,{T}_{\mathrm{eff}}({\boldsymbol{r}}),\,\ldots ),\end{eqnarray} \tag{ 44 }$$

  in order to calculate non-bolometric observables.

The outlined procedure is effectively a reinterpretation of bolometric results in the spectral sense and can only be seen as a rough approximation for the truly wavelength-dependent irradiation framework that would involve ray-tracing the light coming from each surface element to the observer, a complicated and computationally extensive scheme beyond the scope of this paper. The approximate procedure outlined above is currently the standard way of dealing with this technical issue. In order to be consistent throughout the paper, we focus purely on bolometric processes and bolometric observations, but using the described procedure one can also model passband-dependent observations.

We start the discretization by partitioning the surface ${ \mathcal M }$ into distinct subsets:

$$\begin{eqnarray}&&{ \mathcal M }=\bigcup _{i}{{ \mathcal S }}_{i}\quad {{ \mathcal S }}_{i}\bigcap {{ \mathcal S }}_{j}=0\quad i\ne j.\end{eqnarray} \tag{ 45 }$$

with their area equal to

$$\begin{eqnarray}&&{A}_{i}={\int }_{{{ \mathcal S }}_{i}}{dA}({\boldsymbol{r}}).\end{eqnarray} \tag{ 46 }$$

We approximate functions defined on the surface as functions constant over any surface element ${{ \mathcal S }}_{i}$, called the piecewise constant function over ${ \mathcal M }$. A piecewise constant approximation $\tilde{f}$ of an integrable function f defined on the surface ${ \mathcal M }$ is given by

$$\begin{eqnarray}\tilde{f}({\boldsymbol{r}})=\displaystyle \sum _{i}{f}_{i}{\chi }_{i}({\boldsymbol{r}})\quad {\chi }_{i}({\boldsymbol{r}})=\left\{\begin{array}{lll}1 & : & {\boldsymbol{r}}\in {{ \mathcal S }}_{i}\\ 0 & : & \mathrm{otherwise}\end{array}\right.,\end{eqnarray} \tag{ 47 }$$

with the expansion coefficients f_i expressed as

$$\begin{eqnarray}&&{f}_{i}=\displaystyle \frac{1}{{A}_{i}}{\int }_{{{ \mathcal S }}_{i}}f({\boldsymbol{r}}){dA}({\boldsymbol{r}}),\end{eqnarray} \tag{ 48 }$$

where χ_i(r) is a characteristic function of ${{ \mathcal S }}_{i}$ on ${ \mathcal M }$. The set {χ_i(r)} is a functional basis of piecewise constant functions. We treat the vector f = (f_i) as a discretized version of the function f.

Here, the considered operators are linear and therefore can be written in the form

$$\begin{eqnarray}&&\hat{{ \mathcal O }}f({\boldsymbol{r}})={\int }_{{ \mathcal M }}H({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} )f({\boldsymbol{r}}^{\prime} ){dA}({\boldsymbol{r}}^{\prime} ),\end{eqnarray} \tag{ 49 }$$

where H(r, r') is a kernel function that depends on the operator we are considering. We have a piecewise constant function $f={\sum }_{i}{f}_{i}{\chi }_{i}$ with expansion coefficients f_i. We approximate its image $\hat{{ \mathcal O }}f$ with a piecewise constant function ${\sum }_{i}{f}_{i}^{{\prime} }{\chi }_{i}$ with the expansion coefficients given by

$$\begin{eqnarray}&&{f}_{i}^{{\prime} }=\displaystyle \sum _{j}{O}_{i,j}{f}_{j},\end{eqnarray} \tag{ 50 }$$

where the matrix elements O_i,j are expressed as

$$\begin{eqnarray}&&{O}_{i,j}=\displaystyle \frac{1}{{A}_{i}}{\int }_{{{ \mathcal S }}_{i}}{\int }_{{{ \mathcal S }}_{j}}H({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} ){dA}({\boldsymbol{r}}^{\prime} ){dA}({\boldsymbol{r}}).\end{eqnarray} \tag{ 51 }$$

The matrix O = [O_i,j] is the discretized version of the operator $\hat{{ \mathcal O }}$. In the case of the radiosity operator, O_i,j are the generalizations of the view factors (Modest 2013). For the discretized version of the redistribution operator $\hat{{ \mathcal D }}$, denoted by the matrix D = [D_i,j], the flux conservation takes the form

$$\begin{eqnarray}&&\displaystyle \sum _{i}{A}_{i}{D}_{i,j}={A}_{j}.\end{eqnarray} \tag{ 52 }$$

The irradiation framework is implemented by extending the open-source package PHOEBE, where the working surface ${ \mathcal M }$ is a mesh of triangles that approximates the true shape of astrophysical bodies. It supports different geometrical bodies, e.g., aligned and misaligned Roche-shaped stars and isolated rotating stars. The triangular discretization of Roche-shaped bodies was already used in, e.g., Hendry & Mochnacki (2000) and Pribulla (2012), just to name a few. We consider two discretization schemes for operators, per-triangle discretization and per-vertex discretization, as described in Prša et al. (2016b). In both cases, we simplify the expressions for the matrix elements (Equation (51)) to

$$\begin{eqnarray}&&{O}_{i,j}\approx {A}_{j}H({{\boldsymbol{r}}}_{i},{{\boldsymbol{r}}}_{j}),\end{eqnarray} \tag{ 53 }$$

where r_i and r_j are the surface element locations and A_i are the corresponding areas.

In the decomposition of the redistribution operators, we need the distances across the surfaces. The calculation of distances on the triangular meshes is computationally very expensive; see, e.g., Martínez et al. (2005). For astrophysical bodies close to spherical, we can frequently approximate the distances by those on the sphere, which is computationally tractable. Well-detached stars in a binary configuration certainly fall into this category. Consider our object of interest packed inside a sphere of radius R and center c, satisfying

$$\begin{eqnarray}&&\displaystyle \sum _{i}{\left(\parallel {{\boldsymbol{r}}}_{i}-{\boldsymbol{c}}{\parallel }^{2}-{R}^{2}\right)}^{2}=\min .\end{eqnarray} \tag{ 54 }$$

Next, define an operator to obtain the radial unit vector at a point on the sphere w.r.t. the center c,

$$\begin{eqnarray}&&\hat{{\mathbb{P}}}{\boldsymbol{r}}=\displaystyle \frac{{\boldsymbol{r}}-{\boldsymbol{c}}}{\parallel {\boldsymbol{r}}-{\boldsymbol{c}}\parallel }.\end{eqnarray} \tag{ 55 }$$

The geodesic distance between the two points on the mesh, (r₁, r₂), is approximated by the distance between the corresponding points on the sphere:

$$\begin{eqnarray}&&d({{\boldsymbol{r}}}_{1},{{\boldsymbol{r}}}_{2})=R\,\arccos ((\hat{{\mathbb{P}}}{{\boldsymbol{r}}}_{1})\cdot (\hat{{\mathbb{P}}}{{\boldsymbol{r}}}_{2})).\end{eqnarray} \tag{ 56 }$$

Furthermore, if the local curvature of the mesh is small, i.e., the mesh points are dense and the distances between neighboring points considered for local redistribution do not differ much from the geodesic, we can approximate the distances with the Euclidean form:

$$\begin{eqnarray}&&d({{\boldsymbol{r}}}_{1},{{\boldsymbol{r}}}_{2})\approx | | {{\boldsymbol{r}}}_{1}-{{\boldsymbol{r}}}_{2}| | .\end{eqnarray} \tag{ 57 }$$

The deviation from the surface, measured along the axis ${\boldsymbol{\phi }}$, ($\parallel {\boldsymbol{\phi }}\parallel =1$), is then

$$\begin{eqnarray}\begin{array}{rcl}{d}_{\perp }({{\boldsymbol{r}}}_{1},{{\boldsymbol{r}}}_{2}) & = & R\,\arccos \left(\sqrt{(1-{u}_{1}^{2})(1-{u}_{2}^{2})}+{u}_{1}{u}_{2}\right),\\ {u}_{i} & = & {\boldsymbol{\phi }}\cdot \hat{{\mathbb{P}}}{{\boldsymbol{r}}}_{i}.\end{array}\end{eqnarray} \tag{ 58 }$$

The local and latitudinal redistribution models, based on d and d_⊥, are schematically depicted in Figure 2. The flux incident on the surface element centered on the vertex depicted in black is redistributed over the surface elements depicted in yellow, associated with the vertices depicted in red. The fraction of the incident flux that is redistributed over the elements depends on the weight function discussed before.

Zoom In Zoom Out Reset image size

For most practical cases, it suffices to assume discrete redistribution models (i.e., local, latitudinal and global) with constant reflection coefficients. The irradiation parameters ρ_refl, ρ_loc, ρ_lat, and ρ_uni, which are associated with the reflection and the local, latitudinal, and uniform redistributions, respectively, are used to construct operator weights (see Equation (22)):

$$\begin{eqnarray}\begin{array}{rcl}{w}_{\mathrm{loc}} & = & \displaystyle \frac{{\rho }_{\mathrm{loc}}}{1-{\rho }_{\mathrm{refl}}},\quad {w}_{\mathrm{lat}}=\displaystyle \frac{{\rho }_{\mathrm{lat}}}{1-{\rho }_{\mathrm{refl}}},\\ {w}_{\mathrm{uni}} & = & \displaystyle \frac{{\rho }_{\mathrm{uni}}}{1-{\rho }_{\mathrm{refl}}},\end{array}\end{eqnarray} \tag{ 59 }$$

which are non-negative and sum up to 1. In turn, they determine the redistribution matrix for a given object,

$$\begin{eqnarray}&&{\boldsymbol{D}}={w}_{\mathrm{loc}}{{\boldsymbol{D}}}_{\mathrm{loc}}+{w}_{\mathrm{lat}}{{\boldsymbol{D}}}_{\mathrm{lat}}+{w}_{\mathrm{uni}}{{\boldsymbol{D}}}_{\mathrm{uni}}.\end{eqnarray} \tag{ 60 }$$

The principal advantage of these irradiation parameters, i.e., ρ_refl, ρ_loc, ρ_lat, and ρ_uni, is that they add up to 1 for each body separately (assuming no losses), and each individual parameter represents the fraction of the total incoming flux redistributed by the given irradiation process.

Notice that the uniform redistribution matrix for an ith body D_uni,i can be expressed as a projection:

$$\begin{eqnarray}&&{{\boldsymbol{D}}}_{\mathrm{uni},i}=\displaystyle \frac{1}{{A}_{i}}{\left[1,\ldots ,1\right]}^{T}[{A}_{i,1},{A}_{i,2},\,\ldots ,\,{A}_{i,{N}_{i}}],\end{eqnarray} \tag{ 61 }$$

where A_i,j is the area of the jth surface element (j in [1, N_i]) on the ith body and ${A}_{i}={\sum }_{j}{A}_{i,j}$ is the total area of the body. Using this property, we can significantly speed up calculations related to uniform redistribution. This is taken into account in the implementation of the irradiation framework in the extension of PHOEBE 2.1.

We follow the presented discretization procedure and approximate all operators by matrices and all functions defined on the surface by vectors with their entries representing average quantities on surface elements. The discretized limb-darkened ${\hat{{ \mathcal L }}}_{\mathrm{LD}}$ and Lambertian ${\hat{{ \mathcal L }}}_{{\rm{L}}}$ radiosity operators are represented by the matrices L_LD and L_LD, respectively, and the redistribution operator $\hat{{ \mathcal D }}$ is described by the matrix D, and the reflection operator $\hat{{\rm{\Pi }}}$ is approximated by Π. The radiosity F_out, intrinsic exitance F₀, updated exitance ${F}_{0}^{{\prime} }$, and irradiation F_in are approximated by the vectors F_out, F₀, ${{\boldsymbol{F}}}_{0}^{{\prime} }$, and F_in, respectively.

The essential reflection–redistribution equations are given by Equation (23) for Wilson's reflection and by Equation (25) for the Lambertian reflection model. Following the discretization rules, these can be written in matrix form as

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{\mathrm{out}}={{\boldsymbol{G}}}_{{\rm{W}}}+{{\boldsymbol{Q}}}_{{\rm{W}}}{{\boldsymbol{F}}}_{\mathrm{out}}\quad \ \mathrm{for}\ \mathrm{Wilson}^{\prime} {\rm{s}}\ \mathrm{model},\end{eqnarray} \tag{ 62 }$$

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{\mathrm{in}}={{\boldsymbol{G}}}_{{\rm{L}}}+{{\boldsymbol{Q}}}_{{\rm{L}}}{{\boldsymbol{F}}}_{\mathrm{in}}\quad \mathrm{for}\ \mathrm{Lambertian}\ \mathrm{model},\end{eqnarray} \tag{ 63 }$$

where we, for compactness, introduce the auxiliary vectors G_W and G_L given by

$$\begin{eqnarray}&&{{\boldsymbol{G}}}_{{\rm{L}}}={{\boldsymbol{F}}}_{0}\quad {{\boldsymbol{G}}}_{{\rm{W}}}={{\boldsymbol{L}}}_{\mathrm{LD}}{{\boldsymbol{F}}}_{0},\end{eqnarray} \tag{ 64 }$$

and the matrices Q_s written as

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{Q}}}_{{\rm{W}}} & = & \left[{\boldsymbol{D}}({\mathbb{1}}-{\boldsymbol{\Pi }})+{\boldsymbol{\Pi }}\right]{{\boldsymbol{L}}}_{\mathrm{LD}}\\ {{\boldsymbol{Q}}}_{{\rm{L}}} & = & {{\boldsymbol{L}}}_{\mathrm{LD}}{\boldsymbol{D}}({\mathbb{1}}-{\boldsymbol{\Pi }})+{{\boldsymbol{L}}}_{{\rm{L}}}{\boldsymbol{\Pi }}.\end{array}\end{eqnarray} \tag{ 65 }$$

The vectors G_W and G_L represent the intrinsic exitance and limb-darkened irradiated exitance, respectively, whereas the matrices Q_s represent the products of the discretized radiosity operator, reflection coefficients, and the redistribution operator associated with a given reflection–redistribution scheme described by Equations (23) and (25).

By design, the matrices Q_s scale linearly with the discretized radiosity operators and, consequently, its norm is ≤1. Because of that, Equations (62)–(63) are convergent and can be solved iteratively:

$$\begin{eqnarray}&&{{\boldsymbol{F}}}^{(k+1)}={{\boldsymbol{G}}}_{s}+{{\boldsymbol{Q}}}_{s}{{\boldsymbol{F}}}^{(k)}\quad \mathrm{for}\,k=0,1,\,\ldots ,\end{eqnarray} \tag{ 66 }$$

with s = W, L labeling the reflection–redistribution scheme and using the initial condition ${{\boldsymbol{F}}}^{(0)}={{\boldsymbol{G}}}_{s}$. The vector F represents F_out and F_in in Wilson's and the Lambertian reflection model, respectively. Accurate irradiation calculations can be very time consuming. It involves constructing all matrices and solving a large sparse system of linear equations, as described in this section. Therefore, it is useful to have an approximate model to determine the magnitude of irradiation effects in order to decide whether they need to be taken into account given the required precision. We provide a detailed discussion in Appendix B.

We estimate the time complexity of the different phases of the irradiation framework as implemented in the extension of PHOEBE, i.e., triangulating bodies, obtaining properties of the mesh (area of elements, total area, total volume, ...), calculating radiosity operator matrices describing reflection, calculating redistribution operator matrices, and finally, solving the linear system.

Let us assume that we have m convex bodies and the surface of the ith body is partitioned into N_i elements, which in our case are triangles. The generation of a triangular mesh covering the surface of bodies and the calculation of the properties of the mesh are a standard part of PHOEBE 2.1 and are performed in ${\sum }_{i}{ \mathcal O }({N}_{i})$ operations, where ${ \mathcal O }$ signifies the limiting behavior of a function; see, e.g., Širca & Horvat (2018). The radiosity operator matrices L_LD and L₀ are sparse matrices in all practical cases, because the visibility between surface elements is frequently obstructed. Their construction in general takes ${ \mathcal O }({({\sum }_{i}{N}_{i})}^{2})$ operations, but for convex bodies, the number of operations reduces to ${ \mathcal O }({\sum }_{i\gt j}\,{N}_{i}{N}_{j})$. The redistribution matrix D has a block diagonal form, where each block corresponds to a separate body. A block associated with the ith body has a dimension N_i × N_i. In order to calculate all of the blocks, we need ${\sum }_{i}{ \mathcal O }({N}_{i}^{2})$ operations. By setting N_i = N, the computational costs of constructing the redistribution and radiosity operator matrices are equal to ${ \mathcal O }({{mN}}^{2})$ and ${ \mathcal O }(m(m-1){N}^{2})$, respectively. Notice that the latter grows quadratically with the number of bodies, whereas the former grows only linearly. The matrices Q_s (s = W, L) (Equation (65)) are of dimension M × M, where $M={\sum }_{i}{N}_{i}$ is the number of surface elements, and are typically still sparse with ${N}_{\mathrm{nonzero}}={ \mathcal O }({\sum }_{i\gt j}\,{N}_{i}{N}_{j})\,+{\sum }_{i}{ \mathcal O }({N}_{i}^{2})$ non-zero elements. The system of equations determined by the matrices Q_s are solved iteratively. Assuming we need N_it iterations, the solution can be found in N_itN_nonzero operations.

The times needed for the different phases of the irradiation framework as a function of the number of surface elements N = N_i for a simple binary system of two identical Roche-shaped stars are depicted in the Figure 3. Notice that the times needed to generate individual matrices and obtain a solution are of the same order of magnitude and are by far the most costly part of the irradiation framework. These times have a clear quadratic dependence on N in comparison to times to generate the mesh and calculate its properties, which scale linearly with N.

Zoom In Zoom Out Reset image size

Consider a system of two identical spheres with a constant exitance F₀. The consequence of mutual irradiation is the updated exitance ${F}_{0}^{{\prime} }$ and radiosity F of the two spheres. We compute it by using Lambertian reflection and a specific redistribution model with a linear weight function g(x) = 1 − x/l without considering losses. The results are depicted in Figure 4 as a density plot of ${F}_{0}^{{\prime} }$ and F across the surfaces of the spheres. Common to all redistribution models is that the increase of the reflection coefficient decreases the incident flux redistributed over the surface and, as a consequence, a decrease in ${F}_{0}^{{\prime} }$; the increase of the area over which the incoming flux is redistributed makes ${F}_{0}^{{\prime} }$ more uniform across the surface and, on average, decrease in size.

Zoom In Zoom Out Reset image size

In order to highlight the effects of both processes, reflection and redistribution, we choose a reflection coefficient ρ = 0.3, which is large enough to have notable reflection and small enough to enable significant redistribution. Figure 4 shows that the uniform redistribution produces a uniform ${F}_{0}^{{\prime} }$ and represents a bias value for radiosity, as indicated by Equation (18); the local redistribution generally gives the largest ${F}_{0}^{{\prime} }$ and F, and both distributions have a similar shape; and lastly, the latitudinal redistribution increases exitance at latitudes that are most strongly illuminated.

A mean-field approximation of two-sphere case discussed here with a focus on average radiosity and irradiance is presented in Appendix B and can be used to check the order of magnitude of the reflection–redistribution effects.

We apply the introduced irradiation models by calculating LCs for a binary star with redistribution switched both on and off. We consider simplified reflection in which reflection coefficients are constant across a body. For presentation purposes, the bodies are spherical, but the theory is valid for any geometrically defined surface. In addition, we discuss how well we can discriminate between the different redistribution effects based on the LCs.

Consider a model LC with certain redistribution parameters ${\boldsymbol{x}}\in {{\mathbb{R}}}^{p}$ calculated at N time stamps: ${\boldsymbol{C}}({\boldsymbol{x}})=[{C}_{i}({\boldsymbol{x}}){]}_{i=1}^{N}\in {{\mathbb{R}}}^{N}$. In the vicinity of the parameters x (i.e., for perturbed parameters ${\boldsymbol{x}}^{\prime} ={\boldsymbol{x}}+\delta {\boldsymbol{x}}\in {{\mathbb{R}}}^{p}$), we can approximate this LC vector by its Taylor expansion:

$$\begin{eqnarray*}&&{\boldsymbol{C}}({\boldsymbol{x}}^{\prime} )={\boldsymbol{C}}({\boldsymbol{x}})+{\boldsymbol{C}}^{\prime} ({\boldsymbol{x}})\delta {\boldsymbol{x}}+{ \mathcal O }(\parallel \delta {\boldsymbol{x}}{\parallel }^{2}).\end{eqnarray*}$$

The vector norms in use here are denoted by $\parallel \cdot \parallel$ and for a vector x is defined as $\parallel {\boldsymbol{x}}\parallel =\sqrt{{{\boldsymbol{x}}}^{T}{\boldsymbol{x}}}$. The discrepancy between the LC vector at the parameters x and at the perturbed parameters x' is measured by the norm of the difference of the LC vectors, written as

$$\begin{eqnarray}\begin{array}{rcl}\parallel {\boldsymbol{C}}({\boldsymbol{x}}^{\prime} )-{\boldsymbol{C}}({\boldsymbol{x}}){\parallel }^{2} & = & \delta {{\boldsymbol{x}}}^{T}{\boldsymbol{C}}{{\prime} }^{T}({\boldsymbol{x}}){\boldsymbol{C}}^{\prime} ({\boldsymbol{x}})\delta {\boldsymbol{x}}\\ & & +\ { \mathcal O }(\parallel \delta {\boldsymbol{x}}{\parallel }^{3}).\end{array}\end{eqnarray} \tag{ 67 }$$

The discrepancy can be quantified by performing a singular value decomposition (SVD) of the LC vector derivative,

$$\begin{eqnarray}&&{\boldsymbol{C}}^{\prime} ({\boldsymbol{x}})={\boldsymbol{U}}({\boldsymbol{x}}){\boldsymbol{\Sigma }}({\boldsymbol{x}}){{\boldsymbol{V}}}^{T}({\boldsymbol{x}}),\end{eqnarray} \tag{ 68 }$$

where ${\boldsymbol{U}}({\boldsymbol{x}})\in {{\mathbb{R}}}^{N\times N}$ and ${\boldsymbol{V}}({\boldsymbol{x}})\in {{\mathbb{R}}}^{p\times p}$ are orthogonal matrices and ${\boldsymbol{\Sigma }}({\boldsymbol{x}})\in {{\mathbb{R}}}^{N\times p}$ is a diagonal matrix with singular values on the diagonal: the maximum and the minimum value on the diagonal are σ_max and σ_min, respectively. For details on SVD, see, e.g., Širca & Horvat (2018). The discrepancy of the LC vector in the limit of small perturbations is then bound by the singular values

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{\min }({\boldsymbol{x}})\parallel \delta {\boldsymbol{x}}\parallel & \leqslant & \parallel {\boldsymbol{C}}({\boldsymbol{x}}+\delta {\boldsymbol{x}})-{\boldsymbol{C}}({\boldsymbol{x}})\parallel \\ & \leqslant & \ {\sigma }_{\max }({\boldsymbol{x}})\parallel \delta {\boldsymbol{x}}\parallel .\end{array}\end{eqnarray} \tag{ 69 }$$

When σ_min is zero, there is a linear combination of irradiation effects that do not produce changes in the LCs, resulting in a degenerate case. Note that σ_min and σ_max have a dimension and, consequently, they scale linearly with the amplitude of the LC.

Typically, the difference in LCs C(x + δx) − C(x) can be measured up to a certain noise level. Let us denote the discrepancy between the measured and computed LCs as N and treat it as noise.⁶ The changes in the irradiation parameters δx are detectable if

$$\begin{eqnarray}&&\parallel \delta {\boldsymbol{x}}\parallel \geqslant \displaystyle \frac{\parallel {\boldsymbol{N}}\parallel }{{\sigma }_{\max }({\boldsymbol{x}})}\equiv {\epsilon }_{\mathrm{sufficient}},\end{eqnarray} \tag{ 70 }$$

and all changes in the irradiation parameters are measurable if

$$\begin{eqnarray}&&\parallel \delta {\boldsymbol{x}}\parallel \geqslant \displaystyle \frac{\parallel {\boldsymbol{N}}\parallel }{{\sigma }_{\min }({\boldsymbol{x}})}\equiv {\epsilon }_{\mathrm{total}}.\end{eqnarray} \tag{ 71 }$$

To quantify how well we can discriminate between effects, we need to consider how strongly the LC varies due to changes in irradiation parameters about x. We quantify the variation by the ratio between the largest and the smallest responses in the LC vector equal to the condition number (Širca & Horvat 2018)

$$\begin{eqnarray}&&\kappa ({\boldsymbol{x}})=\displaystyle \frac{{\epsilon }_{\mathrm{total}}}{{\epsilon }_{\mathrm{sufficient}}}=\displaystyle \frac{{\sigma }_{\max }({\boldsymbol{x}})}{{\sigma }_{\min }({\boldsymbol{x}})}\geqslant 1.\end{eqnarray} \tag{ 72 }$$

In order to clearly separate the effect of different parameters on the LC, κ needs to be as large as possible and $\parallel \delta {\boldsymbol{x}}\parallel \geqslant {\epsilon }_{\mathrm{total}}$. Note that, in the degenerate case, the conditional number κ is infinite and so is _total.

We are studying irradiation effects in a toy binary system motivated by NN Serpentis, an eclipsing binary system composed of a white dwarf (primary star—P) and red dwarf (secondary star—S) with an orbital period of 0.13 days (Qian et al. 2009; Parsons et al. 2010), although without its recently discovered circumbinary disk (Hardy et al. 2016). The parameters of the toy system are listed in Table 1, where only the red dwarf is subjected to irradiation redistribution. The lobes of the stars in the toy system are described in Roche geometry. The redistribution of irradiation is applied to the system as described in Section 5, with irradiation in this two-body system described by eight parameters: ρ_refl,b, ρ_loc,b, ρ_lat,b, and ρ_uni,b for bodies b = S, P.

Table 1.  Parameters for an NN Serpentis-like System Based on Data in Parsons et al. (2010)

Parameter	White Dwarf	Red Dwarf
Atmosphere	Blackbody	Blackbody
Exponent in gravity brightening, g	1	0.32
Polar radius, R (R⊙)	0.0211	0.147
Effective temperature, Teff (K)	57000	3500
Masses, M (M⊙)	0.535	0.111
Fraction of reflection, ρrefl	1	0.6
Synchronicity parameter, Fsync	1	1
Fillout factor, fRa	0.0472	0.773
Limb darkening (LD):
Model	Logarithmic	Logarithmic
Coefficient xLD	0.5	0.5
Coefficient yLD	0.5	0.5
Orbit:
Period, P (day)	0.1300801714
Eccentricity,	0
Systemic velocity, γ (km s−1)	0
Inclination, ι (degree)	89.6
Mass ratio, M2/M1	0.207
Semimajor axis, a (R⊙)	0.934

Note.

^aThe fillout factor follows the definition in Kallrath & Milone (2009, Ch. 3.1.6).

Download table as:  ASCIITypeset image

We choose the white dwarf's albedo to be ρ_refl,P = 1 and the red dwarf's albedo to be ρ_refl,S = 0.6. This means for the former that ρ_loc,P = ρ_lat,P = ρ_uni,P = 0 and for the latter that it can be subjected to irradiation redistribution effects. The synthetic bolometric LCs and radial velocity curves are calculated in the limiting cases of redistribution in the red dwarf: in the absence of redistribution (ρ_loc,S = ρ_lat,S = ρ_uni,S = 0), where we are confronted with losses; in the presence of only local redistribution (ρ_loc,S = 0.4, ρ_lat,S = ρ_uni,S = 0); in the presence of only latitudinal redistribution (ρ_lat,S = 0.4, ρ_loc,S = ρ_uni,S = 0); and in the presence of only uniform redistribution (ρ_uni,S = 0.4, ρ_loc,S = ρ_lat,S = 0).

We assume that the albedos (i.e., the fraction of the flux that is reflected directly from the surface) of the primary and secondary stars are ρ_P = 1 and ρ_S = 0.6, respectively, and that the radii (l; see Equations (20)–(21)) of the latitudinal and local redistribution processes are given by l/R = 20° ≐ 0.35. The bolometric LCs and radial velocity curves of such a system are presented in Figure 5. Note that the calculation of the approximate passband-dependent models would be possible using Wilson's approach of spectral re-processing of bolometric results, as outlined in Section 4.3, which we elect to forego on account of clarity and consistency.

Zoom In Zoom Out Reset image size

The redistribution increases the emitted flux. The strongest increase in comparison to that without redistribution, around 2%, is noticeable with local redistribution, as we can seen from Figure 5(a). This is because it effectively increases the reflectivity of the object, resulting in an increased radiosity/fluxes around the secondary eclipse, when the reflection is at its maximum. The latitudinal and uniform redistributions have a similar effect on the LCs with an approximate 8% increase from the case without redistribution. Because the uniform redistribution spreads the incoming flux over a wider area (the whole body) than the latitudinal redistribution, the flux of the former case is necessarily smaller than the flux of the latter case.

The redistribution does not affect radial velocities as strongly as fluxes; see Figure 5(c). Interestingly, we obtain a similar radial velocity curves for pairs of latitudinal and uniform redistributions, and local and no redistributions. This is a consequence of the local redistribution affecting only a small area on the surface in comparison to the uniform and latitudinal redistributions.

Differences between the Wilson and Lambertian reflection schemes in LCs and radial velocities are of the order of magnitude 10⁻⁷ and 10⁻⁸ and presented in Figures 5(b) and (c), respectively. Currently, these differences are not measurable in practice. The fluxes obtained using Wilson's approach are larger than those using Lambertian reflection, which seems to be generally true in a binary configuration. This follows from the fact that the limb-darkened (Wilson) diffusion of light amplifies the intensities in the direction nearly normal to the surface in comparison to Lambertian diffusion, where D/D₀ > 1/π for μ ≈ 1 for almost all limb-darkening coefficients. The intensities in the direction nearly orthogonal to the surface are the most important in the transfer of energy between the bodies, yielding an increase in reflected fluxes and consequently an increase in radiosity. The differences between the radial velocity curves obtained using either the Wilson or Lambertian reflection and some redistribution type are pairwise similar in shape for the local and without redistributions, and the uniform and latitudinal redistributions. The differences of the latter pair are generally smaller, because the redistribution to a wider area has the effect of blurring out local differences in radiosity across the lobe.

Following the analysis presented in Section 6.2, we calculate σ_max and σ_min by varying the redistribution parameters of the second star for the cases discussed in Figure 5 and find that σ_max ≈ 0.48, σ_min = 0.00285, and the resulting conditional number κ ≈ 169 (Equation (72)) This suggests that these cases are far from degenerate, and according to Equation (69), the effect of the redistribution in the LC can vary by up to two orders in magnitude, at the some fixed magnitude of the redistribution parameters.

All operators introduced in this paper, i.e., the reflection $\hat{{\rm{\Pi }}}$, redistribution $\hat{{ \mathcal D }}$ and radiosity ${\hat{L}}_{l}$ (l = L, LD) operators, preserve the positivity of the functions, and all functions describing the irradiation process, i.e., irradiance F_in, radiosity F_out, and exitances F₀ and ${F}_{0}^{{\prime} }$, are defined on the surfaces of bodies and are non-negative.

In the mean-field approximation approach, we decompose a function F defined on a surface ${ \mathcal S }$ of area A into its surface average $\langle F\rangle$, defined as

$$\begin{eqnarray*}&&\langle F\rangle =\displaystyle \frac{1}{A}\int F({\boldsymbol{r}}){dA}({\boldsymbol{r}}),\end{eqnarray*}$$

and the deviation δF from it, writing $F=\langle F\rangle +\delta F$. Then, an action of some operator $\hat{{ \mathcal O }}$ on the function F reads

$$\begin{eqnarray*}&&\hat{{ \mathcal O }}F=\hat{{ \mathcal O }}\langle F\rangle +\hat{{ \mathcal O }}\delta F.\end{eqnarray*}$$

Here, for the operators and functions used, the first term on the right is the dominant contribution. In the mean-field approximation, we drop the term depending on fluctuations and approximate the surface average of operator image of an function as

$$\begin{eqnarray*}&&\langle \hat{{ \mathcal O }}F\rangle \approx \langle F\rangle \langle \hat{{ \mathcal O }}\rangle \,\quad \mathrm{and}\quad \langle \hat{{ \mathcal O }}\rangle \equiv \langle \hat{{ \mathcal O }}1\rangle .\end{eqnarray*}$$

We call $\langle \hat{{ \mathcal O }}\rangle$ the average operator and is defined as the surface average of the image of a constant function equal to 1.

Let us discuss a system of two convex bodies, labeled A and B, with Lambertian reflection from the surfaces and introduce the intrinsic exitances F_0,b, updated intrinsic exitances ${F}_{0,b}^{{\prime} }$, and exiting F_out,b and entering F_in,b radiosities defined for both bodies b = A, B. The updated intrinsic exitances ${F}_{0,b}^{{\prime} }$ and exiting radiosities F_out,b are expressed as

$$\begin{eqnarray}\begin{array}{rcl}{F}_{0,b}^{{\prime} } & = & {F}_{0,b}+{\hat{{ \mathcal D }}}_{b}(1-{\hat{{\rm{\Pi }}}}_{b}){F}_{\mathrm{in},b}\quad \mathrm{and}\\ {F}_{\mathrm{out},b} & = & {F}_{0,b}^{{\prime} }+{\hat{{\rm{\Pi }}}}_{b}{F}_{\mathrm{in},b},\end{array}\end{eqnarray} \tag{ 73 }$$

where ${\hat{{ \mathcal D }}}_{b}$ and ${\hat{{\rm{\Pi }}}}_{b}$ are the redistribution operator and reflection operators associated with the body b, respectively. Consequently, we can separate the irradiation Equation (25) into a system of two equations, each dealing with the irradiation of the considered star:

$$\begin{eqnarray}\begin{array}{rcl}{F}_{\mathrm{in},{\rm{B}}} & = & {\hat{{ \mathcal L }}}_{\mathrm{LD},{\rm{A}}\to {\rm{B}}}{F}_{0,{\rm{A}}}+\left[{\hat{{ \mathcal L }}}_{\mathrm{LD},{\rm{A}}\to {\rm{B}}}{\hat{{ \mathcal D }}}_{{\rm{A}}}({\mathbb{1}}-{\hat{{\rm{\Pi }}}}_{{\rm{A}}})\right.\\ & & \left.+\ {\hat{{ \mathcal L }}}_{{\rm{L}},{\rm{A}}\to {\rm{B}}}{\hat{{\rm{\Pi }}}}_{{\rm{A}}}\right]{F}_{\mathrm{in},{\rm{A}}},\end{array}\end{eqnarray} \tag{ 74 }$$

$$\begin{eqnarray}\begin{array}{rcl}{F}_{\mathrm{in},{\rm{A}}} & = & {\hat{{ \mathcal L }}}_{\mathrm{LD},{\rm{B}}\to {\rm{A}}}{F}_{0,{\rm{B}}}+\left[{\hat{{ \mathcal L }}}_{\mathrm{LD},{\rm{B}}\to {\rm{A}}}{\hat{{ \mathcal D }}}_{{\rm{B}}}({\mathbb{1}}-{\hat{{\rm{\Pi }}}}_{{\rm{B}}})]\right.\\ & & \left.+\ {\hat{{ \mathcal L }}}_{{\rm{L}},{\rm{B}}\to {\rm{A}}}{\hat{{\rm{\Pi }}}}_{{\rm{B}}}\right]{F}_{\mathrm{in},{\rm{B}}}.\end{array}\end{eqnarray} \tag{ 75 }$$

The radiosity operators ${\hat{{ \mathcal L }}}_{\mathrm{LD},b\to b^{\prime} }$ and ${\hat{{ \mathcal L }}}_{{\rm{L}},b\to b^{\prime} }$ describe the limb-darkened and Lambertian diffusion of light from the surface of body b onto the surface of body b'. For the sake of simplicity, we assume a constant reflection on both bodies (b = A, B):

$$\begin{eqnarray*}&&{\hat{{\rm{\Pi }}}}_{b}={\rho }_{b}{\mathbb{1}},\quad {\rho }_{b}=\mathrm{const}.\end{eqnarray*}$$

The surface of the bodies can be divided into an illuminated part, called day side, and a non-illuminated part, called night side. The irradiation Equations (74)–(75) describe processes only on the day side, and consequently, F_in,b is non-zero only on that side of body b = A, B. Additionally, we may notice that

$$\begin{eqnarray}&&{F}_{\mathrm{out},b}({\boldsymbol{r}})={F}_{0,b}^{{\prime} }({\boldsymbol{r}})+{\hat{{\rm{\Pi }}}}_{b}{F}_{\mathrm{in},b}({\boldsymbol{r}})\quad {\boldsymbol{r}}\in \mathrm{dayside},\end{eqnarray} \tag{ 76 }$$

$$\begin{eqnarray}&&{F}_{\mathrm{out},b}({\boldsymbol{r}})={F}_{0,b}^{{\prime} }({\boldsymbol{r}})\quad {\boldsymbol{r}}\in \mathrm{nightside}.\end{eqnarray} \tag{ 77 }$$

Next, we introduce the averages over the whole surface $\langle \cdot \rangle$, over the day side $\langle \cdot {\rangle }_{\mathrm{day}}$, and over the night side $\langle \cdot {\rangle }_{\mathrm{night}}$ of a considered body. The areas of the entire surface and of the day side of body b are denoted by A_b and A_b,day, respectively. Taking into account Equation (77), we can conclude for body b that

$$\begin{eqnarray}&&\langle {F}_{\mathrm{out},b}\rangle =p\langle {F}_{\mathrm{out},b}{\rangle }_{\mathrm{day}}+(1-p)\langle {F}_{0,b}^{{\prime} }{\rangle }_{\mathrm{night}},\end{eqnarray} \tag{ 78 }$$

where we use the ratio of the areas p = A_b,day/A_b.

We start the approximation of the irradiation model by decomposing all involved quantities, i.e., F_0,b, ${F}_{0,b}^{{\prime} }$ and F_out,b, into their day- and night-side counterparts. Then, we take the surface average of the irradiation Equations (74)–(75) and the quantities over both sides of the bodies separately. The irradiation equations are only defined on the day sides, and so the night side averages yield zero. We approximate the actions of the operators using a mean-field approach, whereby we replace the functions defined across the day and night sides with their corresponding average; for details, see Appendix B.1. Following this idea, we approximate the averages of the operators acting on a function F defined over the day side of body b according to the next rules:

$$\begin{eqnarray}&&\langle {\hat{{ \mathcal L }}}_{l,b\to b^{\prime} }F{\rangle }_{\mathrm{day}}\approx {L}_{l,b\to b^{\prime} }\langle F{\rangle }_{\mathrm{day}},\end{eqnarray} \tag{ 79 }$$

$$\begin{eqnarray}&&\langle {\hat{{ \mathcal D }}}_{b}F{\rangle }_{\mathrm{day}}\approx {\eta }_{b}\langle F{\rangle }_{\mathrm{day}},\end{eqnarray} \tag{ 80 }$$

$$\begin{eqnarray}&&\langle {\hat{{ \mathcal D }}}_{b}F{\rangle }_{\mathrm{night}}\approx (1-{\eta }_{b})\langle F{\rangle }_{\mathrm{day}},\end{eqnarray} \tag{ 81 }$$

where l = LD, L is labeling different surface behaviors. Here we introduce the model constants L_l,b→b' and η_b that describe the effective action of the radiosity and redistribution operators on the surface-averaged incoming radiosity $\langle {F}_{\mathrm{in},{\rm{b}}}{\rangle }_{\mathrm{day}}$. More precisely, L_l,b→b' represents the average ratio of emitted energy transferred from points on body b to body b', and η_b quantifies the average ratio of absorbed energy that is re-emitted on the same side.

The averages of the quantities describing the irradiation process over the day and night sides are then given by

$$\begin{eqnarray}&&\langle {F}_{0,b}^{{\prime} }{\rangle }_{\mathrm{day}}=\langle {F}_{0,b}{\rangle }_{\mathrm{day}}+(1-{\rho }_{b}){\eta }_{b}\langle {F}_{\mathrm{in},b}{\rangle }_{\mathrm{day}},\end{eqnarray} \tag{ 82 }$$

$$\begin{eqnarray}&&\langle {F}_{0,b}^{{\prime} }{\rangle }_{\mathrm{night}}=\langle {F}_{0,b}{\rangle }_{\mathrm{night}}+(1-{\rho }_{b})(1-{\eta }_{b})\langle {F}_{\mathrm{in},b}{\rangle }_{\mathrm{day}},\end{eqnarray} \tag{ 83 }$$

$$\begin{eqnarray}&&\langle {F}_{\mathrm{out},b}{\rangle }_{\mathrm{day}}=\langle {F}_{0,b}^{{\prime} }{\rangle }_{\mathrm{day}}+{\rho }_{b}\langle {F}_{\mathrm{in},b}{\rangle }_{\mathrm{day}},\end{eqnarray} \tag{ 84 }$$

$$\begin{eqnarray}&&\langle {F}_{\mathrm{out},b}{\rangle }_{\mathrm{night}}=\langle {F}_{0,b}^{{\prime} }{\rangle }_{\mathrm{night}}.\end{eqnarray} \tag{ 85 }$$

By introducing additional auxiliary coefficients,

$$\begin{eqnarray}&&{T}_{b\to b^{\prime} }=(1-{\rho }_{b}){\eta }_{b}{L}_{\mathrm{LD},b\to b^{\prime} }+{\rho }_{b}{L}_{{\rm{L}},b\to b^{\prime} },\end{eqnarray} \tag{ 86 }$$

$$\begin{eqnarray}&&{G}_{b^{\prime} }={L}_{\mathrm{LD},b\to b^{\prime} }\langle {F}_{0,b}{\rangle }_{\mathrm{day}},\end{eqnarray} \tag{ 87 }$$

the average of the irradiation Equations (74)–(75) can be rewritten into a simple system of two scalar equations involving variables $\langle {F}_{\mathrm{in},b}{\rangle }_{\mathrm{day}}$ for body b = A, B:

$$\begin{eqnarray}&&\langle {F}_{\mathrm{in},{\rm{A}}}{\rangle }_{\mathrm{day}}={G}_{{\rm{A}}}+{T}_{{\rm{B}}\to {\rm{A}}}\langle {F}_{\mathrm{in},{\rm{B}}}{\rangle }_{\mathrm{day}},\end{eqnarray} \tag{ 88 }$$

$$\begin{eqnarray}&&\langle {F}_{\mathrm{in},{\rm{B}}}{\rangle }_{\mathrm{day}}={G}_{{\rm{B}}}+{T}_{{\rm{B}}\to {\rm{A}}}\langle {F}_{\mathrm{in},{\rm{B}}}{\rangle }_{\mathrm{day}}.\end{eqnarray} \tag{ 89 }$$

The solution of this system is equal to

$$\begin{eqnarray}&&\left[\begin{array}{c}\langle {F}_{\mathrm{in},{\rm{A}}}{\rangle }_{\mathrm{day}}\\ \langle {F}_{\mathrm{in},{\rm{B}}}{\rangle }_{\mathrm{day}}\end{array}\right]=\displaystyle \frac{1}{1-{T}_{{\rm{B}}\to {\rm{A}}}{T}_{{\rm{A}}\to {\rm{B}}}}\left[\begin{array}{c}{G}_{{\rm{A}}}+{T}_{{\rm{B}}\to {\rm{A}}}{G}_{{\rm{B}}}\\ {T}_{{\rm{A}}\to {\rm{B}}}{G}_{{\rm{A}}}+{G}_{{\rm{B}}}\end{array}\right].\end{eqnarray} \tag{ 90 }$$

Let us now assume that the bodies are perfect spheres of radii r_A and r_B, and their centers separated by distance d, as depicted in Figure 6. We set the intrinsic exitance F_0,b to be constant over the surface of each star separately. We are interested in the limit r_b ≪ d in which we can approximate model constants as

$$\begin{eqnarray}&&{L}_{l,b\to b^{\prime} }\approx \displaystyle \frac{{r}_{b}^{2}}{2{d}^{2}},\end{eqnarray} \tag{ 91 }$$

$$\begin{eqnarray}{\eta }_{b}\approx \left\{\begin{array}{lll}\displaystyle \frac{1}{2} & : & \mathrm{global}\ \mathrm{or}\ \mathrm{latitudinal}\ \mathrm{redistribution}\\ 1 & : & \mathrm{local}\end{array}\right..\end{eqnarray} \tag{ 92 }$$

We find numerically that the coefficient L_l,b→b' is independent of the type of energy diffusion, labeled by l. An explicit formula for the coefficient is given in Appendix B.3. In the following, we compare the average exiting radiosity $\langle {F}_{\mathrm{out},b}\rangle$ and average updated exitance $\langle {F}_{0,b}^{\prime} \rangle$ as functions of the distance d between the bodies obtained in a one-dimensional model, given by

$$\begin{eqnarray}&&\langle {F}_{\mathrm{out},b}\rangle ={F}_{0,b}+\displaystyle \frac{1}{2}\langle {F}_{\mathrm{in},b}{\rangle }_{\mathrm{day}},\end{eqnarray} \tag{ 93 }$$

$$\begin{eqnarray}&&\langle {F}_{0,b}^{\prime} \rangle ={F}_{0,b}+\displaystyle \frac{1}{2}(1-{\rho }_{b})\langle {F}_{\mathrm{in},b}{\rangle }_{\mathrm{day}},\end{eqnarray} \tag{ 94 }$$

and that obtained from numerical calculations by discretizing the surface into triangles, as described in Section 5. In the considered limit, we may take the area of the day and night sides to be identical; for details, see Appendix B.3. The results are depicted in Figure 7 and show a good qualitative agreement between the two approaches, especially in the limit of larger separations between bodies.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Let us consider a binary system composed of two spheres, labeled A and B, with radii r_A and r_B, and their centers separated by a distance d, as depicted in Figure 6. The area of the illuminated side of the sphere b is given by

$$\begin{eqnarray}&&{A}_{b,\mathrm{day}}=2\pi {r}_{b}^{2}\left(1-\displaystyle \frac{{r}_{b}^{{\prime} }-{r}_{b}}{d}\right).\end{eqnarray} \tag{ 95 }$$

We assume that the limb-darkening law is constant across the surface. The coefficients L_l,b→b' with l = L, LD introduced in Appendix B describing the effect of the radiosity operator on the average incoming radiosity can be identified as the averages of radiosity operator, based on the previous subsection:

$$\begin{eqnarray}&&{\left\langle {\widehat{{ \mathcal L }}}_{l,b\to b^{\prime} }\right\rangle }_{\mathrm{day}}\equiv {L}_{l,b\to b^{\prime} }.\end{eqnarray} \tag{ 96 }$$

In the considered case, this coefficient can be expressed as an integral over all possible pairs of points on the two spheres with their line of sight unobstructed:

$$\begin{eqnarray}\begin{array}{rcl}{L}_{l,b\to b^{\prime} } & = & \displaystyle \frac{{\left({r}_{b}{r}_{b^{\prime} }\right)}^{2}}{{A}_{b^{\prime} ,\mathrm{day}}}\displaystyle \int d{\rm{\Omega }}({\hat{{\boldsymbol{n}}}}_{b})\\ & & \times \ \displaystyle \int d{\rm{\Omega }}({\hat{{\boldsymbol{n}}}}_{b^{\prime} })U({\rm{\Delta }}{\boldsymbol{r}}\cdot {\hat{{\boldsymbol{n}}}}_{b})U({\rm{\Delta }}{\boldsymbol{r}}\cdot {\hat{{\boldsymbol{n}}}}_{b^{\prime} })\\ & & \cdot \ \displaystyle \frac{({\rm{\Delta }}{\boldsymbol{r}}\cdot {\hat{{\boldsymbol{n}}}}_{b})({\rm{\Delta }}{\boldsymbol{r}}\cdot {\hat{{\boldsymbol{n}}}}_{b^{\prime} })}{\parallel {\rm{\Delta }}{\boldsymbol{r}}{\parallel }^{4}}\displaystyle \frac{{D}_{l}(\widehat{{\rm{\Delta }}{\boldsymbol{r}}}\cdot {\hat{{\boldsymbol{n}}}}_{b})}{{D}_{l,0}},\end{array}\end{eqnarray} \tag{ 97 }$$

where integrations are carried out over full solid angles, and we should remind ourselves that D_l is the limb-darkening function and D_l,0 is its integral over a hemisphere:

$$\begin{eqnarray}&&{D}_{l,0}=2\pi {\int }_{0}^{1}d\mu \,{D}_{l}(\mu )\mu .\end{eqnarray} \tag{ 98 }$$

Here we use the unit-step function U(x) = {1: x ≥ 0; 0: otherwise}, and the vector connecting the pairs of points is equal to

$$\begin{eqnarray}&&{\rm{\Delta }}{\boldsymbol{r}}={r}_{2}{\hat{{\boldsymbol{n}}}}_{2}+d\hat{{\boldsymbol{k}}}-{r}_{1}{\hat{{\boldsymbol{n}}}}_{1}.\end{eqnarray} \tag{ 99 }$$

We find numerically that the leading order of the average operator in the limit $d\to \infty$ behaves as

$$\begin{eqnarray*}&&{L}_{l,b\to b^{\prime} }\sim \displaystyle \frac{{r}_{b}^{2}}{2{d}^{2}},\end{eqnarray*}$$

and this behavior seems to be independent of the chosen limb-darkening law labeled by the index l.