Hyper Illumination of Exoplanets: Analytical and Numerical Approaches

This work describes the illumination of exoplanets whose orbits are close enough to their host star that the finite angular size of their host star causes hyper illumination, in which more than 50% of the planet receives light. Such exoplanets include the hot Jupiters KELT-9 b (64.5% illuminated) and Kepler-91 b (69.6% illuminated). We describe the geometry of three primary illumination zones: the fully illuminated zone, penumbral zone, and unilluminated zone. The integrals required to determine the incident radiation as a function of position from the substellar point on the exoplanet are explained and derived, and the analytical solution is presented within the fully illuminated zone. We find that the illumination predicted by our model is greater at the substellar point than the typical plane-parallel ray model used would suggest. In addition, it is greater within the region of the penumbral zone extending into the antistellar side of the exoplanet. Finally, we compare our model to that used in starry, an open-source software package used to create albedo maps. It appears that starry may be overestimating the illumination of closely orbiting exoplanets because the foreshortening of the area element of the host star is not included in its calculation.


Introduction
Over 5000 exoplanets have been discovered and characterized using a variety of methods, including precise photometric data from instruments such as Kepler, which allow researchers to estimate planetary parameters such as albedo and temperature. 5To do so, scientists must first create mathematical models of the system in question.In this work, we concern ourselves with the modeling of the illumination of an exoplanet.The illumination of an exoplanet by its host star determines its atmospheric characteristics by contributing to the energy balance of the atmosphere, and the light reflected and absorbed by the planet.
For most exoplanetary systems, it is sufficient to assume that the planet is far enough away from its host star that its rays arrive as plane-parallel rays, as shown in Figure 1.In which case, the illumination of the exoplanet is greatest at the substellar point and decreases with the cosine of the zenith angle (Cowan et al. 2013, Equation (38)).Many methods of modeling the phase curve of exoplanets for the purposes of parameter estimation and albedo mapping make use of this assumption (Madhusudhan & Burrows 2012;Farr et al. 2018;Haggard & Cowan 2018;Berdyugina & Kuhn 2019;Kawahara 2020;Heng et al. 2021;Teinturier et al. 2022).The assumption has also been used in ocean glint studies by comparing a planet's apparent albedo at different orbital phases (Lustig-Yaeger et al. 2018).It is also implemented as part of 1D climate modeling code, such as PICASO (Batalha et al. 2019), and 3D climate modeling code such as ROCKE-3D (Way et al. 2017) and the Community Atmosphere Model (Neale et al. 2010;Wolf 2017).
Yet, for exoplanets that orbit very close to their host star the plane-parallel ray illumination assumption breaks down due to the finite angular size of the star and its area elements, and the presence of penumbral effects.Here, we define such exoplanets as being hyper illuminated.In Lillo-Box et al. (2014), the authors note that the hot Jupiter Kepler-91 b was approximately 70% illuminated due to its proximity to its host star.In addition, Guillot (2010) note that the use of the plane-parallel ray approximation would result in an overestimation of the dayside/nightside temperature difference for exoplanets experiencing more than 50% illumination.In Léger et al. (2011), the authors compare the modeled temperature distribution of CoRot-7 b using the plane-parallel ray approximation to a model that accounts for the finite angular size of its host star in their Figure 4, noting that the temperature effects within the penumbral region are significant.To correctly characterize such exoplanets, the spherical and finite angular size of an exoplanet's host star must be taken into account, as demonstrated in Knuth et al. (2017) and Carter (2018).
The reflection effect between binary stars serves as a starting point.In Sen (1948), the author describes how to approach the problem if one of the stars is treated as a point source rather than assuming plane-parallel rays.Later, Kopal (1959) proposed an analytical solution for spherical stars of finite size, which requires an approximation out to fourth order in the ratio of the stellar radii to their separation of centers.In addition, a full numerical treatment of the reflection effect between binary stars is described in Wilson (1990), including multiple reflection effects in which the reflected light by one star causes heating in its companion star.The approximation used in Kopal (1959) is planned to be implemented in the Transit and Light Curve Modeller (Csizmadia et al. 2023).In Kopal (1959), the author considers the reflected light, and hence illumination, of the larger of the two stars.Carter (2018) noted that this analysis of the reflected intensity within regions of the star that are partially illuminated does not extend to that of exoplanets because exoplanets are always smaller than their host star; therefore, changes to the limits of integration must be implemented.
Finally, a numerical integration of an exoplanet's illumination due to a star of finite size is implemented within starry (Luger et al. 2022) to allow for albedo mapping of transiting exoplanets via the analysis of reflected light.However, the method used in starry neglects the foreshortening of the stellar area elements included in Kopal (1959) and Carter (2018).In this work, Section 2 describes the physical theory required to analyze hyper-illuminated exoplanets, including the extent of the illumination in Section 2.1, the illumination due to a point source in Section 2.2.1, and the extension to a spherical source in Section 2.2.2.The analytical solution within the fully illuminated zone is described in Section 2.3, and numerical approaches in Section 2.4.We highlight our results in Section 3, using three exoplanets as examples utilizing theoretical illumination maps.Finally, we conclude by exploring the significance of our results in Section 4. A list of symbols is given in Table 4 in the Appendix.

Theory
In the following, we outline the geometry and calculations required to determine the light incident, i.e., the illumination, on a point due to a source in the context of points on a planet illuminated by its host star.Here, the planet and star are treated as spherical.When referenced, the exoplanetary radius (R p ) is taken to be the distance from its center to its surface, where we treat the top of the atmosphere as the effective surface for exoplanets with atmospheres.In the following, we restrict ourselves with the illumination due to the finite angular size of the host star and do not consider the refraction of light by the atmosphere of the planet, which shall be left to future work.
We begin by considering the extent of the illumination of the exoplanet in Section 2.1, where we present a table of planets which are hyper illuminated.Next, the method of calculating the illumination will be described in Section 2.2.We conclude with a presentation of the analytical and numerical solutions in Sections 2.3 and 2.4.

Illumination Extent
In this work, we are concerned with the illumination of one sphere (an exoplanet) by a spherical illumination source (a star).First, we consider in broad terms how much each part of the exoplanet is illuminated by the star using Figure 2. Points within the inner tangents between the star and exoplanet receive light from the entire apparent disk of the host star, and we shall refer to this region as the fully illuminated zone.Points beyond the outer tangents are in the unilluminated zone and receive no light from the host star.Points between the two are in the penumbral zone, and receive light from part of the host star.
The penumbral zone is further split into two subzones.The division occurs at η pen in Figure 2, defined by points whose horizon intersect with the center of the host star.Points in the starside penumbral zone receive light from more than half of the apparent disk of the host star (PZ1), the others from less than half the apparent disk (PZ2).The division between the two is important when defining integration limits, as shown in Carter (2018, Sections 6.2 and6.3) andPerera &Knuth (2023).
The angles shown in Figure 2 may be determined in terms of the parameters of the system using trigonometry and are as follows (Luger et al. 2022, Equation (4) ;Kopal 1959;Carter 2018, Equation (5.3.89) where R s is the radius of the host star, R p is the radius of the exoplanet, and r s is the separation of centers of the two.
With the aid of the angles defined in Equation (1), one can calculate the fractional surface area of each zone as , where A zone is the surface area of a given zone.To determine the area illuminated (σ lit ), it is easiest to subtract the unilluminated area of the spherical cap defined by η un , i.e., A R 2 1 sin , from the surface area of the planet.Doing so reveals that the fractional area illuminated is given by R R r If the star is treated like a point source, then we may assume that R s /r s = 1, and we find that the fractional area illuminated is R r 0.5 1 , which is equivalent to the fractional area of the spherical cap delimited by η pen , i.e., with apex angle 90°− η pen measured from the substellar point.If we further assume that the exoplanet is small compared to the star-planet separation, we naturally find that the fractional area illuminated is one-half and return to the plane-parallel ray approximation used in works such as Sobolev (1975).
Finally, we present in Table 1 a series of confirmed exoplanets, with their parameters and excess illumination.The excess illumination is the fractional surface area over onehalf illuminated, or σ excess = σ lit − 0.5.The area of excess illumination is calculated using A R 4 p excess excess 2 . Planets appearing in Table 1 include hot Jupiters, super-Earths, sub-Neptunes, and the terrestrial planet Trappist-1 d, illustrating the range of planet types impacted by hyper illumination.

Illumination Calculation
In the following two sections, we describe the method used to calculate the illumination of each point on a planet due to a point source (Section 2.2.1), then alter it to apply to an extended illumination source (Section 2.2.2).

Illumination by a Point Source
We begin by considering the geometry between a point illumination source located at the center of the host star with Cartesian coordinate location r x x y y z z s s s s ˆˆ= + + from the center of the exoplanet located at O p , and a point P located on the surface of the exoplanet at p xx yy zz ˆˆ= + + in Figure 3.The point-source illumination will be used for comparison, and to build up the formulas used to determine that of an extended source.
For a point-source illuminator like that shown in Figure 3, the intensity (I pt ) of light received at p is given by for points on the exoplanet for which the point source is visible.In Equation (3), πS is the flux per unit normal, which may be   Note.Listed are the default parameters for a selection of confirmed exoplanets obtained from the NASA Exoplanet Catalog (2023 July 13), along with their excess illumination (σ excess ) and the surface area of the excess illumination (A excess ).We have used the semimajor axis of the orbit as r s .For reference, the total area of India, 3.29 × 10 6 km 2 (CIA 2023), is comparable to the excess illumination of Trappist-1 d.
written as πS = L s /ρ 2 for a point source of luminosity L s at distance ρ from p (Kopal 1959).From Figure 3, we see that ρ = r s − p, where where θ is the latitude of p, f is the longitude of p −δ is the polar angle of r s , and −θ s is the azimuthal angle of r s .The foregoing coordinate system is set such that the illumination maps in Section 3 are oriented so that +y-axis is up the page, the +x-axis is to the right, and the +z-axis is out of the page; therefore, θ is measured from the x-z plane with positive values up the page, and f from the +z-axis in the x-z plane with positive values to the right.
-, such that the dot product in Equation (3) may be expressed as From the foregoing, we find that where y x y sin cos and b is the semiminor axis of the ellipse defining the orientation of the hemisphere facing the illumination source.Combining Equations (3)-( 5) and normalizing to the stellar luminosity (L s ) reveals that Here, we are making use of the max() function to set the intensity to zero for points on the exoplanet that cannot see the illumination source.
For star-planet separations much larger than the planetary radius, Equation (5) reduces to cos cos i x J » , and ρ ≈ r s in Equation ( 7).In such cases, ϑ i is equivalent to the zenith angle (g¢ in Figure 4), which is the spherical angle between p and the substellar point.The intensity of light is then approximately that given by the plane-parallel ray approximation: in units of the stellar luminosity (Sobolev 1975).

Illumination by an Extended Source
To apply the foregoing to extended sources, i.e., those of finite angular size, we will follow Kopal (1959) and Carter (2018, Section 5.3).We neglect ellipsoidal variations and restrict our exploration to spherical illumination sources, as depicted in Figure 4. Here, a single point on the host star's surface is represented by Q, and our task is to determine the total intensity of light received at point P on the exoplanet from all visible Q's.A point Q is visible from P so long as cos x and cos q¢ are greater than zero.
For each Q visible at P, the intensity of light received is given by where I Q is the intensity of light emitted by a point Q, p ˆ• r¢ accounts for the projection of that intensity in the direction of P, and dΩ Q,P is the differential solid angle of Q as viewed from P. Each point Q represents the center of some area element of the star dσ Q ; therefore, we make use of I Q dΩ Q,P , rather than simply S as was done in Equation (3).In this way, we account for the way in which the differential solid angles of the area elements of a sphere contribute to the intensity, as opposed to a single point's contribution.
Next, the differential solid angle is given by the projected differential surface area divided by the separation between the area element of Q and P squared (Boas 2006, Equation (10.20)).Using Figure 4, see that the projected differential surface area is given by dr dA and the separation is r¢, so that The projection of the differential surface area accounts for the foreshortening of regions with differential surface area dA Q as viewed at P. Finally, using Figure 4, we see that Inserting the foregoing into Equation (9) for extended illumination sources.The quantity max 0, cos cos x ¢ ensures that only Q visible at P contribute to the illumination of P. In principle, Equation (11) may be used to account for limb darkening by allowing I Q to vary with cos q¢ (Mandel & Agol 2002).In this work, we will largely restrict ourselves to cases of zero limb darkening and leave its inclusion to future work.
It is important to note that for large star-planet separations, we may assume that dr cos , cos cos x g » ¢ and r s r¢ » .In which case, we return to the plane-parallel ray approximation as described by Equation (8) after integration over the surface area of the star and normalization (Carter 2018, Section 5.2).

Analytical Solutions
We turn now to the analytical integration of Equation (11).Doing so requires one to determine cos x, cos q¢, and r¢ in terms of the spherical coordinates of Q as described in Kopal (1959), Carter (2018, Section 5.3), and Perera & Knuth (2023).It is convenient to use a coordinate system centered at P with azimuthal angle ω Q and polar angle θ Q , each measured from r.In such a system dA R d d sin , where R s is the radius of the host star.The shift in our coordinate system to one centered at P, with zˆr = in Figure 4, gives the following line of cosines for dr Q ˆ: dr Our task is to determine cos x and cos q¢ in terms of coordinates on the host star (R s , θ Q , and ω Q ), and constants for the point P (ρ, α, and g¢).
First, it is useful to determine the separation of Q and P, r¢ in Figure 4, using the law of cosines for !QO s P: R R 2 cos .13 Further, we find that the new definition of cos x in Figure 4 requires that In addition, it may be shown that Making use of Equations ( 14)-( 15), we find that We will now insert Equation (16) into Equation (11) and add integration limits, allowing us to remove the max() function by ensuring that the limits are such that cos cos 0 x q¢ > .It is convenient to write this integral as a summation of two integrals as follows: , ; , cos , 17 18b ò ò w w q q r q r q r r q q q w = -- The limits of integration depend on the location of the point P because its location determines the nature of the intersection of the horizon with the stellar disk.
As shown in Kopal (1959) and Carter (2018, Section 5.3.2), the integration must be performed over a spherical cap for points within the fully illuminated zone because there is no intersection of the horizon with the stellar disk.The limits on ω Q are then simply ω 0 = 0 and ω f = 2π, and we will find it convenient to rewrite the integral in terms of r¢ instead of θ Q .The integration within this zone is analytical and will be described in Section 2.3.1.
For points within the penumbral zone, we must account for the horizon's intersecting with the apparent disk of the host star.This may be described using the geometric depth of the eclipse (λ), as shown in Figure 5.In our coordinate system, the law of sines and some trigonometric identities can be used to show that where the positive sign is used for points in PZ1 and the negative sign for points in PZ2 (see Figure 2).The division between the zones occurs at λ = 0, which is delimited by a spherical cap on the planet whose apex angle is 90°− η pen .
The division into two penumbral zones introduces unique challenges to the integration that are neglected in Kopal (1959), as discussed in (Carter 2018, Section 5.3.3-5.3.4).The integration of Equation ( 11) is not analytical within the penumbral zones due to the dependence of the integration limits on λ 1,2 ; therefore, expansion approximations must be used (Carter 2018).An approximation to fourth order of R p /r s , R s /r s , and products of the two will be left to future work.

Fully Illuminated Zone
We turn now to the intensity of incident light for points within the fully illuminated zone, which are the points within a spherical cap centered at the substellar point with apex angle 90°− η full (see Figure 2).
For such points, the integral A[0, 2π; θ 0 , θ f ] (Equation ( 18a)) evaluates to zero because the integration over of cos Q w is zero.In addition, the evaluation of Equation (18b) over ω Q becomes ò p q q p r q r q r r q q q = ´-- because there is no dependence of I Q on ω Q .We may simplify B[0, 2π; θ 0 , θ f ] by solving Equation (13) for cos Q q and differentiating in terms of r¢ to find that The minimum value of r¢ occurs for Q on ρ in Figure 4 and is given by R s min r r ¢ = -.Its maximum value occurs for Q such that r¢ is tangent to the surface of the star, or max r ¢ = R s 2 2 r -.The integration of Equation (18b) then becomes for points in the fully illuminated zone.
Thus far, we have not restricted ourselves to stars of uniform intensity and could account for limb darkening by allowing I Q to vary with cos q¢, and hence on r¢.We will now restrict ourselves to uniform stars with zero limb darkening, but note that Equation ( 22) is analytical for linear limb darkening (Kopal 1959).Letting I Q = I 0 = constant, we find that the final evaluation of Equation (18b) in the fully illuminated zone is Application of trigonometry reveals that cos 2 a r may be written as Finally, the intensity of light as a function of position on the exoplanet in units of the stellar luminosity ) within the fully illuminated zone is given by In Section 3, we present instellation maps within the fully illuminated zone produced using Equation (26).

Numerical Solutions
To determine the incident light outside of the fully illuminated and unilluminated zones, one must use either analytical approximations (Kopal 1959;Carter 2018, Sections 5.3.3 and 5.3.4Here, we outline numerical integration approaches to the problem.We begin by describing a numerical approach to evaluating Equation (11) then continue with the approach used in starry (Luger et al. 2022), and finally compare the two.

Evaluating I ext,cos cos
x q¢ To evaluate Equation (11) numerically, we distribute N source points uniformly over the effective disk of the host star, as is done in starry (Luger et al. 2022, Section 4.1), and let dA Q = 1/N.Then we distribute points over the exoplanet in a grid.For each pair of points P and Q we calculate dI ext,cos cos x q¢ from Equation (11), and determine the total illumination, I ext,cos cos x q¢ , for each P by summing over all Qʼs.Example instellation maps are presented in Section 3.

starry Solution
Here, we will consider only the solution presented in Luger et al. (2022) for surfaces with zero roughness and uniform albedo because we are only concerned with the illumination model used in starry, not reflected light or occultations.
First, to determine the illumination at a point on the exoplanet, starry uses (Luger et al. 2022, Equations (A11) and (A12)) I r 1 max 0, cos , 27 where ϑ i is the angle between the normal of p and the location of the point source, or r s in Figure 3.6 Note that Equation (27) differs from Equation (7), but the two are approximately equivalent for cases in which R p = r s .
To evaluate the illumination by an extended, spherical source, starry averages over the contribution of a set of N source points uniformly distributed over the apparent disk of the host star (Luger et al. 2022, Section 4.1).In doing so, starry effectively sums over N point sources using Equation (27) with the replacement of r s with |r Q | = |r s + dr Q | (see Figure 4) and averages over the visible points by multiplying by dA Q = 1/N. 7f one were to neglect foreshortening provided by assuming cos 1 q¢ » in Equation (11), i.e., dI I dA max 0, cos , 28 x one would arrive at a similar solution to that given in starry.
In cases for which R p = r s , r Q r¢ » , and cos cos i x J » , Equation (28) reduces to Equation (27) with the replacement of r s = r Q .In Section 3, we include a comparison of I ext,cos x to I ext,cos cos x q¢ and I ext,starry .

Results
Here, we explore the intensity calculated by each of the methods outlined in Section 2. In the following, we always take r R z s s = , so that the substellar point is located at p R z p = .To account for planetary spin and obliquity, one need only apply the appropriate rotations for the exoplanet's rotation about the substellar point.
In Section 3.1, we begin with intensity maps of latitude versus longitude normalized to the (constant) stellar intensity (I star ) in parts per million (ppm).Following the intensity maps, we will present subtractions between the methods as a function of longitude along the intensity equator, which is a great circle perpendicular to the plane of the page in Figure 4 and passing through r s (latitude = 0°).Next, we consider the intensity at the substellar point for each model as the star-planet separation is varied to ensure that they converge for large star-planet separations.We conclude in Section 3.2 by exploring the effects of varying the number of source points, N, on the computations of I ext,starry and I ext,cos cos x q¢ , by comparing them to I full and their computation times.

Intensity Maps and Subtractions
In the following, we will be using three exoplanets for illustrative purposes and whose parameters are described in Table 2. KELT-9 b and Kepler-91 b are both hot Jupiters with illumination in excess of 60% of their surface area.In addition, we consider a terrestrial planet, Trappist-1 d (Gillon et al. 2017;Turbet et al. 2020;Agol et al. 2021), whose illumination is much closer to 50%.Throughout the remainder of this section, we use the mean values for planetary and stellar parameters, and calculate the intensity values assuming a star-planet separation equal to the semimajor axis of the planetary orbit unless otherwise noted.
To begin, the intensity maps shown in Figures 6-8 are labeled as follows.First, those using numerical methods are labeled I ext,starry (Luger et al. 2022), I ext,cos cos x q¢ (Equation ( 11)), and I ext,cos x (Equation ( 28)), and were generated using 25 source points distributed uniformly over the effective stellar disk.Three analytical methods are shown for comparison.The I full maps show the results produced by Equation (26) within the fully illuminated zone only, the I pt,starry method uses the point-source approximation given in Equation ( 27), and the I plane method is generated using the plane-parallel ray approximation given in Equation (8).
Shown in Figure 6 are intensity maps for KELT-9 b generated using the six methods described above.We note that the intensity predicted by I ext,starry exceeds that of the analytical model everywhere within the fully illuminated zone.The maximum difference occurs at the substellar point and is 118.8263ppm of I star .In addition, the I ext,cos cos x q¢ model Figure 6.Intensity maps for KELT-9 b using six methods to calculate the illumination.In the figure, I pt,starry and I ext,starry use starry's point-source and extendedsource models, and I ext,cos x uses Equation (28) as described in Section 2.4.2.I ext,cos cos x q¢ uses the methods described in Section 2.2 by Equations ( 9), (10), and (11), and I full is the analytical solution to Equation (11) within the fully illuminated zone given by Equation (26).For the numerical models (I ext,starry , I ext,cos cos x q¢ , and I ext,cos x ), 25 source points were used.Note that I ext,starry and I ext,cos x suggest greater illumination compared to the other methods, with I ext,starry overestimating the intensity at the substellar point by 118.8263 ppm of I star compared to I full .See text for descriptions of each method listed in the maps.The black dashed line delimits the fully illuminated zone.The white dotted lines are contours lines, where 100% is the maximum intensity across all six maps.matches that of the analytical model, I full , within the fully illuminated zone well.The I ext,cos cos x q¢ model also predicts a more rapid decrease in intensity moving from points near the substellar point toward the unilluminated zone than that predicted by I ext,starry .Finally, we note that the I ext,cos x method, which neglects foreshortening as described in Section 2.4.2, is a better match to I ext,starry , indicating the foreshortening to be a major contributor to the difference between starry's method and that presented in this work (see Equation (11)).
The intensity maps for Kepler-91 b are shown in Figure 7.The maximum difference between I starry and I full is much less than that of KELT-9 b at only 21.2204 ppm of I star .The maximum difference occurs at the substellar point, and I starry predicts greater illumination over a larger area of the planet than both I ext,cos cos x q¢ and I full , as was seen for KELT-9 b.Both the fully illuminated and unilluminated zones of Kepler-91 b are smaller than those of KELT-9 b due to the differences in the radii of the exoplanet, star, and orbit.
The last set of intensity maps in Figure 8 shows the illumination for Trappist-1 d.This exoplanet's orbit is much larger than that of KELT-9 b and Kepler-91 b when scaled to its host star's size.As a result, we find that the differences between the six solutions are less pronounced, with a maximum of only 0.0242 ppm of I star .Note that the decreased brightness of the I ext,cos cos x q¢ method compared to the I full method is likely due to the very small magnitude of illumination, so that numerical errors are proportionally greater.Here, we find that starry performs well, potentially because Trappist-1 is more point-like, i.e., cos cos i x J » and cos 1 q¢ » , and starry makes use of the point-source approximation that I ext,cos cos x q¢ does not.We explore these differences later in Section 3.2., where we discuss how the models perform as the number of source points are varied.
Figure 9 shows subtractions between the six methods explored in Figure 6 for KELT-9 b along the intensity equator.Note that the largest differences consistently occur at the substellar point.Further, differences are greatest between those methods that do or do not include foreshortening, e.g., I ext,starry − I ext,cos cos x q¢ .As expected, the differences shown by the I ext,cos cos x q¢ − I full line are much smaller, indicating a good match between the numerical and analytical solutions within the fully illuminated zone.For exoplanets with extremely close-in orbits, like KELT-9 b, the plot of I ext,cos cos x q¢ − I plane in Figure 9 shows that neglecting the finite angular size of the host star will result in underestimating the illumination at the substellar point and within PZ1.The subtraction of I full − I plane is similar to that of I ext,cos cos x q¢ − I plane within the fully illuminated zone, again indicating a good match between the numerical and analytical solutions and supporting the assumption that the numerical solution within the penumbral zone would be consistent with an analytical one.
In Figure 10, we present plots of the substellar intensity estimated using the six methods using the parameters for KELT-9 b given in Table 2, except varying the value of the semimajor axis.As expected, the difference between the six methods decreases with increasing star-planet separation because the star's angular size approaches that of a point source.It is only for very close-in exoplanets that we expect hyper illumination to have a strong influence on the incident light, and hence the reflected light, of an exoplanet.

Varying N
Next, we explore the effect of varying the number of source points on the host star (N) used in the numerical calculation of I ext,starry and I ext,cos cos x q¢ .We begin by exploring the instellation along the intensity equator.We will then explore the instellation at the substellar point as a function of N.
In Figure 11, we plot the outputs of I ext,starry and I ext,cos cos x q¢ along the intensity equator, and I full within the fully illuminated zone.The values of N are chosen such that the number of points can be uniformly distributed over the effective stellar disk.
From Figure 11, we observe that the greatest instellation and the greatest deviation from I full occurs at the substellar point (longitude = 0°).In addition, we see that both I ext,starry and I ext,cos cos x q¢ match I full within the fully illuminated zone well for N = 1 and for N > 4. Yet, I ext,starry does not approach I full as well as I ext,cos cos x q¢ for KELT-9 b or Kepler-91 b, for which the foreshortening of area elements on the host star's limbs will be greater.It is worth recalling that I ext,cos cos x q¢ provides estimates of the instellation beyond the fully illuminated zone and the hemisphere facing the host star for N > 1, which motivates its use over I full and I plane which do not apply over the whole exoplanet.
To better understand the effect N has on the precision of the numerical integration, we next explore the absolute differences The black-and cyan-dotted vertical lines delimit the unilluminated and fully illuminated zones, respectively, where the penumbral zone is located between the two lines.See legend for subtractions.The inset is a zoomed-in plot within the illuminated region to highlight the smaller differences.
Figure 10.Intensity at the substellar point for KELT-9 b analogs using six models of illumination.In all cases, the models converge for large star-planet separations.The green dashed line marks the scaled star-planet separation given by a/R s using the parameters given in Table 2. See text for legend labeling conventions.
potentially because the illumination source is more point-like and the approximations used in Equation (27) better apply.We note that the differences for Trappist-1 d are likely not detectable with current instruments.
Finally, we wish to explore the balance between computation time and the precision granted by increasing the number of source points used in the numerical calculation.We do so by comparing the difference between I ext,starry , I ext,cos cos x q¢ , and I full in Table 3 with computation time.Similar to Figure 12, we see that increasing the number of source points results in the outputs of I ext,starry and I ext,cos cos x q¢ approaching a nonzero difference to I full .Given the information presented in Table 3, we recommend that 25 source points be used to calculate the illumination as granting the ideal balance between computation time and precision using this paper's numerical integration, I ext,cos cos x q¢ .

Discussion and Conclusions
In this work, we derived the equations required to determine the illumination of an exoplanet due to a spherical host star as a function of position on the planet.Importantly, we include the foreshortening of the surface area of the host star in determining an expression for the solid angle of its surface element in Equation (10).We show that neglecting this component results in an overestimation of the illumination, especially at the substellar point (Figures 6-9).
To our knowledge, all global climate modeling and atmospheric models assume that the incident radiation from the host star reaches the top of the atmosphere as plane-parallel rays.This results in the instellation being maximized at the substellar point, decreasing to zero with cos g¢, and the exoplanet experiences only 50% illumination.This work shows that the plane-parallel ray assumption underestimates the illumination at the substellar point, and neglects radiation extending beyond the starside of the exoplanet for small star-planet separations, such as that of KELT-9 b.In future work, we will provide the models for incident light for a given grid resolution or type to the global climate or atmospheric models to determine the effects on model outputs.
In addition, future work will explore the potential impact of the modified instellation near the poles of terrestrial exoplanets orbiting M-dwarf stars, like Trappist-1 d, using a global climate model such as ROCKE-3D.Way et al. (2021) show that the extent of polar ice can influence global mean surface temperatures by several degrees due to ice-albedo feedback.The changes to the instellation near the poles due to hyper illumination may affect polar ice in a low-obliquity world in a way that has yet to be explored.
To improve the estimation of albedos and temperature distributions of exoplanets determined using light-curve data, one must account for the reflected light of the exoplanet.In the case of hyper-illuminated exoplanets, light will not only be reflected from the starside, but also from the part of the antistarside of the exoplanet within PZ2, as shown in Figure 2. To determine the normalized reflected luminosity of the exoplanet as a function of phase angle (ò), one would integrate Equation (17) multiplied by the reflection intensity distribution ( , , ( ) g g f ¢ ), and the cosine of the angle of foreshortening (cos g) over the visible area of the exoplanet: , cos .29 The evaluation of Equation (29) must be performed in a piecewise fashion due to the piecewise nature of the integration , respectively.Odd values of N are plotted with unfilled markers, and even values with filled markers.Note that all of the averages are nonzero, but they approach zero for larger star-planet separations.Both models demonstrate a "wiggle" of the absolute error.I ext,starry performs better for large star-planet separations for which the host star is more point-like.Notes.Table listing the difference in parts per million (ppm) of the stellar intensity between the numerical and analytical models of the substellar point illumination.KELT-9 b is used as the test case.Times for starry include initializing the maps and calculating the intensity.All sets were run on an Intel ® Core™ i3-4160T CPU @ 3.10 GHz and are the average time of 10 calculations.Given our findings, we recommend using 25 source points.
limits of Equation (18), which varies with each illumination zone.
The evaluation of Equation (29) within the fully illuminated zone is not analytical, and instead must be calculated numerically, or by using an approximation by expanding Equation (26) using Legendre polynomials (Kopal 1959;Carter 2018).Future work will include calculating Equation (29) within each zone and incorporating the results into modeling software such as EXONEST (Knuth et al. 2017;Carter 2018;Perera & Knuth 2023).Doing so will improve estimates of planetary parameters such as temperature distributions, which can be degenerate with reflected light (Placek et al. 2014;Perera & Knuth 2023), scattering albedo, and planetary radius.2020), theano (The Theano Development Team et al. 2016), symPy (SymPy Development Team 2023).

Appendix Table of Symbols
In Table 4 we provide a list of symbols used in this work that includes each symbol, a description of the symbol, and a reference to the symbol's use in the text.

Figure 1 .
Figure 1.Example of the plane-parallel ray model of illumination of an exoplanet in which the yellow arrows represent the rays of light incident on the exoplanet at right.For such exoplanets, the instellation is proportional to the cosine of the zenith angle (g¢) because all angles of incidence (ξ) are equal to the zenith angle.The illumination extends to cover only the hemisphere of the exoplanet facing its host star.

Figure 2 .
Figure 2. Shown are the illumination zones of an exoplanet (right) due to illumination by a spherical star (left) of radius R s illuminating a spherical exoplanet of radius R p , whose centers are separated by distance r s .The angles shown delimit the fully illuminated zone (η full , in yellow), the penumbral zones (η pen ), and the unilluminated zone (η un , in black).The penumbral zone nearer the host star (orange) is denoted PZ1 in the text, and that farther from the host star (maroon) is denoted PZ2.The inner tangents are associated with the dashed lines, the outer tangents with the dotted-dashed lines, and the penumbral zones with the dotted lines.Figure adapted from Carter (2022).
Figure 2. Shown are the illumination zones of an exoplanet (right) due to illumination by a spherical star (left) of radius R s illuminating a spherical exoplanet of radius R p , whose centers are separated by distance r s .The angles shown delimit the fully illuminated zone (η full , in yellow), the penumbral zones (η pen ), and the unilluminated zone (η un , in black).The penumbral zone nearer the host star (orange) is denoted PZ1 in the text, and that farther from the host star (maroon) is denoted PZ2.The inner tangents are associated with the dashed lines, the outer tangents with the dotted-dashed lines, and the penumbral zones with the dotted lines.Figure adapted from Carter (2022).

Figure 3 .
Figure 3. Illustration of a point P on an exoplanet with origin O p located at p xx yy zz ˆˆ= + + being illuminated by a point-source host star at r x x y y z z s s s s ˆˆ= + + with associated directions, vectors, and coordinates.

Figure 4 .
Figure 4. Illustration of a point P on the planet being illuminated by a point Q on the host star with associated directions and lengths.The point Q is located at distance r q from O p and r¢ from P. The star's center is located at O s a distance r s from the center of the planet O p .Note that cos i J is now given by r p Q ˆ• ˆ, and cos x by p ˆ• r¢, similar to Figure 3, but with the shift of r s to r Q = r s + dr Q and ρ to r d Q r r ¢ = + .Adapted from Carter (2018).

Figure 5 .
Figure 5. Shown is the geometry required to determine the geometric depth of the eclipse (λ), which defines the polar angle on the host star at which the horizon at P intersects the apparent stellar disk.A: the horizon intersects the host star such that more than half of the apparent disk is visible.This situation applies for points within PZ1 for which 2  a p , δ = π − |λ| − |β|, and 2 | | b p a = -.B: the case for P for which the horizon blocks more than half of the apparent stellar disk, i.e., for points within PZ2.In this case 2 a p > , δ = π − |λ| − |β|, and 2 | | b a p = -.Adapted from Carter (2018).
Parameters of the test case planets used in this work.The illumination was calculated using Equation (2) assuming the mean values for the planetary and stellar radii, and semimajor axis given above.References.(a) Gaudi et al. (2017), (b) Esteves et al. (2015), (c) Berger et al. (2018), and (d) Agol et al. (2021).

Figure 7 .
Figure 7. Intensity maps for Kepler-91 b using six models of illumination.See text and Figure 6 for labeling conventions.As in Figure 6, I ext,starry overestimates the illumination compared to I full .

Figure 8 .Figure 9 .
Figure 8. Intensity maps for Trappist-1 d using six models of illumination.See text and Figure 6 for labeling conventions.Note that the differences between the models are less pronounced than seen in Figures 6 and 7 due to the much smaller value of R p /r s .

Figure 11 .
Figure 11.Shown are intensity calculations using I ext,starry and I ext,cos cos x q¢ along the intensity equator for three exoplanets and varying the numbers of source points N. Included for reference is the analytically determined values within the fully illuminated zone I full .See legend for line styles.From top to bottom, the exoplanets are KELT-9 b, Kepler-91 b, and Trappist-1 d.Within the fully illuminated zone, we see that I ext,cos cos x q¢ performs better than I ext,starry for the two hot Jupiter exoplanets.Both models perform well for Trappist-1 d, with I ext,starry matching better.

Figure 12 .
Figure 12.The absolute differences |I ext,starry − I full | and I I ext,cos cos full | | x q¢ for three exoplanets are shown at the substellar point.From top to bottom, the exoplanets are KELT-9 b, Kepler-91 b, and Trappist-1 d.The dashed line and dotted-dashed line show the average values of |I ext,starry − I full | and I I ext,cos cos full | | x q¢ Figures3 and 4 Figure 2

Table 3
Numerical Method Comparison