Novel instellation calculations for close-in exoplanets

Context. The distribution of instellation at the top of a planet’s atmosphere or surface is usually calculated using the inverse-square law of radiation. This is based on the assumption that the host star is far enough to be considered a point-sized Lambertian source. This assumption, which works well for the solar system planets and most exoplanets, must be revised for close-in exoplanets. Aims. Our objective is to derive accurate instellation patterns for close-in exoplanets, for which the e ﬀ ects of the spherical shape of the star must be taken into account. Methods. We ﬁrst derived an analytical formula of the insolation as a function of substellar longitude, taking the star and the planet as 3-D bodies, and incorporating the limb darkening e ﬀ ects of the star. We then developed a numerical model to compute the distribution of the insolation on close-in planets as a function of substellar longitude for a wide range of stellar and planetary properties. Results. We observe signiﬁcant deviations in instellation values and distribution on close-in exoplanets, due in similar proportions to the physical size of the star and stellar limb-darkening e ﬀ ects. We ﬁnd that the insolation at the susbtellar point is always higher – by as much as 21 % for known exoplanets – than the standard calculation. We ﬁnd that the substellar longitude of the terminator can signiﬁcantly extend on the nightside, from 90 ◦ up to 110 ◦ for known exoplanets. This deviation from standard calculations is attributed to the incorporation of the physical size of the star. The size of the star is a relative quantity based upon how close we are to the star. Due to the proximity of close-in exoplanets to their host stars, the angular size of the star as seen from these exoplanets can be as high as 36 degrees, that is, 72 times larger than the Sun as seen from earth. Hence, it becomes imperative to consider the star as a sphere to compute the instellation. We present and provide the revised instellation and zenith angle patterns for close-in rocky planets (LHS 3844b, Kepler-10 b, TRAPPIST-1b, etc.) and gaseous planets (KELT 9b, WASP 18b, WASP 121b, etc.) of interest for in-depth characterization. Our python code InstellCa for calculating insolation on planets is made available to the community for further studies.


Introduction
Instellation is one of the major factors that govern the climate of an exoplanet. Precise instellation values and their distribution above the atmosphere is therefore essential for climate modelling. The inverse-square (IS) law of instellation suffices for a vast majority of the extra-solar planet population. It turns out that the IS law is based on a few facile assumptions. With our approach, we attempt to reform these assumptions. To understand these, we assume the star as a point-sized object at the origin of our coordinate system demonstrated as a bright yellow spot in Figure 1. The distance of the star to the substellar point in this assumption is elementary and given by a − R p where a is the semi-major axis that is, (R p + R + d) viz., the sum of planetaryand-stellar radius and the distance between the surfaces. We can then write the IS law as: where L b is the bolometric luminosity. This equation has been used for previous studies of close-in exoplanets to construct atmospheric thermal structures (Stevenson et al. 2014). Since the star has a physical size, we see the distance of the substellar point from the surface of the star is not a constant but instead variable at every point on the star. Additionally, due to a large distance, the incident radiation in IS law is assumed to have a planar wave-front. Due to this, the angle subtended by the radiation is dependent on only one variable that is, the latitude (θ p ). Also, the IS law doesn't account for the limb darkening effects of the star. For exoplanets with semi-major axes less than ten times their hoststar radius, these assumptions break down. When we derive an equation taking close-in exoplanets into account, it is imperative to reform these assumptions. The paper is structured as follows.
In section 2, we first derive an equation for instellation using a 3-D geometrical formulation and limb darkening effect for a surface element on the star and then integrate it for the whole star. The integration of this equation by analytic means happens to be complicated and cumbersome. We then use numerical integration methods to solve this integral. In the results section, we have explored the effects of our model for various close-in Hot Jupiters, rocky and eccentric planets. Subsequently, we also discuss the cases where our model significantly differs from the IS law and the reason behind this difference. Finally, we discuss a few ways our model can be tested with observational data and further scope for improvement.

Method
The IS law is well defined for a point-sized source. A spherical body can be assumed to be made up of many such pointsize sources. To calculate the net irradiance or instellation due to a spherical source(star) on a planet, we need to integrate the radiation emanating from all points on this star. Previous studies by Wilson (1990), Budaj (2011), Chen & Rhein (1969), and Horvat et al. (2019 have used a fairly similar approach to calculate the bolometric radiance for modeling the light curves of close binary stars. They also incorporate the reflection of radiation between the binary components to estimate the radiosity. Our model doesn't include the reflection effects, and we do not intend to calculate the radiosity and radiance. Instead, we calculate the irradiance due to the star on the planet. We work strictly with bolometric quantities. Also, our model uses non-linear limb darkening law as opposed to the linear and quadratic laws used in previous models. The work by Wood (1973) is useful to understand the essence of our model. Wood (1973) uses an approach very similar to ours. However, in his analysis, he aims to establish an analytical relation using quasi-empirical methods without integration. With our approach, we attempt to include all the intricacies of the problem and take advantage of the advanced computational power of numerical methods of integration to derive our results. We get the net instellation upon integrating all values for the stellar surface visible to the planet. Assuming the star and planet to be spherical, we choose an arbitrary point on the star and calculate the irradiance on an arbitrary point on the planet as illustrated in Fig. 1.
2.1. Description of the star-planet system in 3-D geometry We can express any point in Cartesian coordinates on the stellar sphere as: where θ and φ define the latitude and longitude of the host star, respectively. Likewise, we can express any point in Cartesian coordinates on the planetary sphere as: x p = R p cosθ p sinφ p y p = R p + R + d + R p cosθ p cosφ p = a + R p cosθ p cosφ p z p = R p sinθ p where θ p and φ p define the latitude and longitude of the planet, respectively. Here, d is the shortest distance between the surfaces of the planet and the star and (R p + R + d) is the semi-major axis(a). The star is assumed to have a total bolometric luminosity L b and effective temperature T . For a surface element dS on The standard equation for irradiance or the power per unit area incident on a surface with an arbitrary area vector, for a source with power dL, is expressed as: (McCluney 2014) We calculate these instellation values along the substellar longitude on the planet; this makes the calculation simple since the planet is considered as a 2-D body. We also need to incorporate the subtle geometrical factors and limb darkening effects for such close systems. We discuss these effects in the following subsection.

Apparent surface
The integration limits of the instellation on a close-in planet depend on the star-planet distance. The equation determining the angle on the stellar surface visible to the planet is obtained using trigonometric arguments as: Therefore, we have:

Stellar surface above local horizon
For higher latitudes on close-in exoplanets, a part of the stellar surface is hidden beneath the local horizon. Our model takes this into account by integrating the surface above the local horizon. An incident ray grazing through the local horizon would be perpendicular to the normal vector. Therefore, for obtaining the angle for the local horizon, we put cosA=0 in equation 14.
Since this effect happens along the latitude, it would only affect the polar angle(θ ). When we solve specifically for θ , we get a quadratic equation with the following roots, This equation determines the limit of integration for higher latitudes. Therefore, equation 3 now becomes: The equation 3 in this case would be: ∀π − θ p > θ 1 the above equation would be used since for higher latitudes we have some part of the star hidden by the horizon.
In addition to this, we know that the planet also has a definite size, and this causes the lower limit to increase from zero since the pole of the planet is above the equatorial plane of the star. Therefore, ∀π − θ p > θ 3 equation 4 becomes, The angle θ 3 can be calculated using basic trigonometry as, However, this factor only comes into play when the planet is much closer to its host star. In fact, so close that the distance between the stellar and planetary surface is comparable to the planetary radius.

Limb Darkening
Hitherto we have assumed the star to be a perfect black-body. However, the stellar photosphere is prone to absorption and scattering of the emitted light that adds a level of intricacy to the former assumption. A star appears to be more bright around the centre and less bright near the limbs. In other words, there exists a temperature gradient across the stellar photosphere. As a result, each surface element on the star has a different temperature. Therefore, it becomes imperative to include the effect of this gradient. We use the non-linear limb darkening model proposed by Claret (2000) due to its accuracy and applicability to the whole HR diagram. The equation for limb darkening in the Claret (2000) model is given by: Further we define, Here I is the bolometric intensity, a 1 , a 2 , a 3 , a 4 are bolometric limb darkening coefficients, and µ is given by, Where, ψ is the angle between the line of sight and the normal to the surface of the star, as shown in Figure 1. This angle can be calculated using equation 13 by replacing the normal to the planetary surface by normal to the stellar surface. For all practical purposes, we can ignore the radius of the planet to make the calculation easy.

Final equation for instellation distribution
If we include the stellar limb darkening effects and other geometrical effects, then our equation culminates to: Now, the luminosity is also a function of ψ. The luminosity of a surface element on the star depends on the optical depth at that surface and hence the temperature of the surface.
Here, T o is the temperature at an optical depth of 1. The analytical solution for this integral is intricate to obtain due to the complex nature of the limits. It is possible to use standard solutions of elliptic integrals to solve the integral analytically. However, it cannot be done without introducing over-simplistic assumptions, and it is cumbersome to solve. Therefore, we use numerical methods to solve the integral.
We calculate the limb darkening coefficients according to a given temperature, surface gravity and metallicity using the Claret(2000) database. We use the limb darkening coefficients computed using the PHOENIX model for stellar temperatures from 2000 to 9800 K, log of surface gravity ranging from 3.5 to 5, and assuming solar metallicity. We perform a multi-linear interpolation to obtain the required coefficients.

Zenith angle calculation
We calculate the mean zenith angle weighted over the area and irradiance at a given latitude(θ) on the planet:

Instellation plots
We present the results for a Kepler-10 b, LHS 3844 b, TRAPPIST-1 b and CoRoT-7 b in figure 3. Additional plots are presented in the appendix.

Instellation at the poles
An important distinction between the two approaches comes at the poles of the exoplanet. The IS law predicts zero instellation since it follows a cosine relation of instellation with latitude. Therefore, for a latitude of 90 degrees, the instellation becomes zero. Our model takes into account the precise angle made by all surface elements on the star with the point on the planet under consideration and hence gives precise, non-zero instellation values. The contour plot-1 in figure 4 shows how the ratio between our model and the standard case varies with varying stellar radii and semi-major axes. We see higher values of the ratio at the upper-left extreme of the plot since the IS law heavily underestimates the instellation at the poles.

Maximum latitude on the night side
Another interesting aspect of this model is the distribution of instellation on the planetary hemisphere not facing the star. Unlike the IS model, the instellation here extends towards the night-side of the exoplanet. Contour plot-2 shows the maximum latitude up to which instellation is non-zero. If we consider the latitude to go beyond its maximum value of 90 degrees along the sub-stellar longitude, then we see that for some exoplanets the instellation becomes zero at nearly 110-degree latitude or 20 degrees from the pole to the night-side of the exoplanet. We define this angle as the nocturnal threshold or the angle beyond which we see the night side. For far-off planets, the angle remains 90 degrees, and this is indicative of the fact that the IS model suffices for such exoplanets.

Discussion
Our model makes the theoretical basis for estimating the amount of irradiation in the proximity of the stellar surface. The results of our model can be easily tested for the Sun. The Parker solar probe which is bound to reach its perihelion at nearly 0.04 AU from the Sun in December 2025 would be useful to obtain observational data at a close distance to Sun. The Poynting flux measured by the FIELDS instrument on the solar probe would be useful, as an observational benchmark, to verify the results. We also suspect that our results can have a significant impact on GCM models for close-in exoplanets. Our model predicts significant instellation at higher latitudes on the night-side hemisphere of the exoplanets. We surmise that this might have some significance in the phase offsets observed for Hot-Jupiters like WASP-12b. An earlier study by Schwartz et al. (2017) aims to understand these phase offsets, and it would be interesting to see how their model might be affected if they use our revised instellation patterns. Also, our results might have some significance for explaining the brightness offset in phase curves of WASP-18b Arcangeli et al. (2019).

Conclusion
Our model currently derives the limb darkening coefficients using the database provided by Claret (2000). However, there have been more recent studies in this area, for instance, Claret & Bloemen (2011) that have not calculated bolometric coefficients since observations are wavelength specific. We could improve our model if we obtain revised bolometric limb darkening coefficients from these recent studies.
To better understand the behaviour of these patterns, we need an analytical estimation of instellation as a function of the stellar radius, semi-major axis, limb darkening coefficients and latitude. This would be useful to get insights to understand the most extreme cases that might exist but have not been observed. There is a scope for solving the integral without the limb darkening effects for the sub-stellar point. However, the inclusion of limb darkening adds greater complexities to the problem. Due to this, understanding the role of limb darkening on the instellation patterns becomes intricate since the final integral becomes too complicated to be solved using standard elliptic integral solutions. We need a comprehensive mathematical solution for the integral to interpret the results better.

Assumptions in adopting the limb darkening model
In the following section, we discuss the assumptions we took to incorporate the limb darkening effects in our model. The effective temperature is the observational parameter that can be obtained from exoplanet catalogues. However, our model needs the central temperature of the star. Therefore, we need a relation between the effective temperature and central temperature T 0 . For obtaining this, we assume that the bolometric luminosity for the sun is not affected by the limb darkening considerations. In other words, the power emitted by the star should be the same regardless of what limb darkening model we choose.

S
σT Here, the LHS is the total bolometric luminosity of the star, and we use numerical methods to solve the integral on the RHS. There are also some caveats when we compare the metallicity definitions in the Claret model with the observational metallicity data. The number of α elements are known for all the stars of the Hypatia catalog and also for the sun (Asplund et al. 2009). The α In other words, the instellation values are not very sensitive to small changes in metallicity. Fig.3 shows how the additional factor varies for the stars of the Hypatia catalog (Hinkel et al. 2014)

Orbital Dynamics
We have included the orbital elements into our model to include the variations due to different positions in an orbit. For highly eccentric planets, there is a significant variation in periastron and apastron instellation. For such planets, we provide the annual mean of instellation over one complete orbit.

True anomaly
The true anomaly is essential to calculate the distance of the planet from its host star at any given instant in its orbit.

Eccentricity
The eccentricity of a planetary orbit is related to the true anomaly and radius as: where ν is the true anomaly and a is the semi-major axis. If we incorporate this additional factor in equation 10, then we can obtain the instellation on the planet according to its position in the orbit.

Obliquity of rotational axis
The declination is related to the obliquity and true anomaly as, The declination will determine the extent of the shift along the latitude axis in the plot of instellation vs. latitude. This is merely a change in the reference angle. For the periastron(ν = 0), the declination would be equal to the obliquity and peak of instellation values will be at a latitude equivalent to the obliquity instead of 0-degree latitude.

Additional Results
For Earth, this model is in accordance with the IS model and correctly predicts the value of the solar constant to be 1354.04 W/m 2 which is well in agreement with the experimental solar constant measurements (Fröhlich 2000).  The InstellCa code for computing instellation values is written in Python, and it uses numerical integration techniques from the scientific python(SciPy) library. It reads the required input from the observational data available at the Extrasolar planet encyclopedia or the NASA exoplanet archive. As an output, it generates a curve of Instellation vs. sub-stellar longitude in addition to the curve for the standard point-size approximation of the star to enable the user to compare the difference between the two approaches.
-The code needs multiple parameters to run. It reads these parameters from catalogue files available at either NASA Exoplanet archive or Exoplanet.eu. The user needs to download the current Exoplanet data from these databases in a csv format. For NASA archive, the data should be downloaded for all columns. -The InstellCa package consists of the Python code, user guide, exoplanet observational data and limb darkening data. The downloaded csv file should be saved in the exoplanet catalog sub-directory of InstellCa package. -The code can be run for multiple exoplanets depending on the user input.
-The code asks the user for a choice between manual entry of parameters or to read from exoplanet.eu or NASA archive.
-Subsequently, the code asks the user to input the name of the exoplanet. The user should carefully mention the name exactly as mentioned in the catalog file with attention to case sensitivity. -The code needs three crucial parameters to run: semi-major axis, stellar radius, stellar effective temperature. If either of these is missing in the catalog the code stops and displays the warning message "Crucial parameter missing". In this case, the user should manually input the missing parameters. -If the planetary radius, stellar mass or metallicity is missing in the catalog data, the code specifically asks for these. If the eccentricity is missing, the code assumes a circular orbit. -In case the user needs to enter a different value for one or two parameters but does not want to input all the other values, the code gives an option to change one or two values. For changing one of the values the input is as follows: -Stellar radius:1 planetary radius:2 effective temperature:3 semi-major axis:4 For changing two values the user needs to feed the following input: -Stellar radius, planetary radius:12 -Stellar radius, effective temperature:13 -Stellar radius, semi-major axis:14 planetary radius, effective temperature:23 planetary radius, semi-major axis:24 effective temperature, semi-major axis:34 -The code also asks the user for choosing the type of limb darkening law to used for the instellation calculation. If the user chooses the Claret(2000) non-linear law, the code automatically computes these coefficients using pre-existing data. In case the pre-existing data is insufficient, the user is asked to enter all the coefficients. For linear and quadratic laws, the user needs to input the bolometric coefficients. -For highly eccentric orbits(e>0.3) the code computes an annual mean of instellation over one complete orbit.
-As an output, the code generates a plot of Instellation vs. sub-stellar longitude and saves it according to a user-defined filename.