Constraining Mornings and Evenings on Distant Worlds: A new Semianalytical Approach and Prospects with Transmission Spectroscopy

Our semianalytic approach to the problem has been implemented in the catwoman library (Jones & Espinoza 2020), which is fully documented ⁴ and available on Github. ⁵ In practice, catwoman's code base is that of batman (Kreidberg 2015), and as such, the library inherits most of the high-level functionalities of this latter library. A catwoman light curve thus receives as input the time-of-transit center t₀, the period P of the orbit, the scaled semimajor axis a/R_*, the inclination i of the orbit with respect to the plane of the sky, the eccentricity e and argument of periastron ω of the orbit, and a set of limb-darkening coefficients for any of the laws already available in batman. On top of these, catwoman takes as input the radii of each of the stacked semicircles, R_p,1 and R_p,2, and the angle φ between the axis that connects them and the vector that follows the direction of motion in the orbit (see Figure 1).

The motivation behind allowing one to define the angle φ in the light-curve generation comes from the possibility of being able to detect the sky-projected spin–orbit misalignment of the planet, which is something the eclipse mapping technique for both light curves (Williams et al. 2006; Rauscher et al. 2007) and radial velocities (Nikolov & Sainsbury-Martinez 2015) is able to do in principle. As will be shown in Section 3.3, detecting the effect of asymmetric light curves due to morning/evening terminator structural and/or compositional inhomogeneities almost guarantees the possibility of putting constraints on this angle, and ignorance of its value does not have a great impact on the detectability of the effect. One important point to consider on this parameter is that this defines the instantaneous angle between the axis that joins the semicircles and the direction of the orbital motion (see Figure 1; orbital motion indicated with a dashed arrow). Because orbits as projected in the plane of the sky are curved in general, this means the axis that joins the semicircles rotates when compared against a straight line projected in this plane. This effect has been implemented within catwoman as well (see Appendix); we validate this implementation against a numerical implementation in the next subsection.

In order to validate the semianalytical approach presented here and implemented in the catwoman library, we built a numerical model that is also able to generate asymmetric light curves due to terminator inhomogeneities but through a completely independent and straightforward (albeit "brute-force") approach. While by construction, catwoman is able to reach any desired precision level (as that is a parameter that can be modified and is tested for convergence before running light-curve model evaluations), our objective with this alternative approach is to validate and illustrate that catwoman is indeed able to generate asymmetric light curves with accuracies of at least 10 ppm, which are the noise limits ultraprecise spectrophotometers like the upcoming JWST will be able to reach (Greene et al. 2016). Our implementation of this numerical scheme is also available in Github ⁶ and is detailed below.

Our approach to this numerical version of the light curve generation of asymmetrical transits is very similar to that of von Paris et al. (2016) and consists of simply discretizing the plane of the sky into n_p × n_p "pixels" centered around the target star. Pixels within the star are filled with values between zero and 1 according to a given intensity profile I(μ)/I(1), while positions that include either the planet or the sky are filled with zeros. The precision on the light curves generated by this scheme can thus be optimized by simply increasing n_p . In practice, this is implemented by populating a matrix of dimensions (n_p , n_p ), on which we first fill all pixels within a distance of n_p /2 from the center of this matrix (i.e., (n_p /2, n_p /2)) with intensities given by the defined intensity profile (a quadratic law in the case of our implementation); all other pixels are filled with zeros. With this, we sum all of the pixel values to compute our out-of-transit flux. Our algorithm then, using as input the coordinates of the center of the planet with respect to a reference frame centered on the star (X, Y) at each time step and the input angle φ, computes the slope of the orbital motion s = dY/dX by simple differences at each time step i, i.e., s_i = (Y_i+1 − Y_i )/(X_i+1 − X_i ). This is then used to compute the instantaneous rotation of the axis that joins the stacked semicircles with respect to the orthogonal system that defines the (X, Y) positions as $\arctan ({s}_{i})+\varphi$. This axis is then used to separate the areas covered by both semicircles, pixels inside of which are set to zero.

We use this simple numerical scheme to validate the semianalytical framework developed in this work by computing a set of cases including a challenging one in which the planetary orbit is significantly curved. This latter case allows us, in turn, to verify that our method outlined in Appendix is correctly accounting for the rotation of the axis that joins the semicircles with respect to the orthogonal system that defines the (X, Y) positions of the planet. To generate this, we simulate an exoplanet with a period of 3 days, time-of-transit center t₀ = 0, scaled semimajor axis a/R_* = 1.5, and zero eccentricity for three inclinations: i = 55°, 72°, and 89°. For the star, we define a quadratic limb-darkening profile with u₁ = 0.3 and u₂ = 0.2. As for the physical properties of the planet, we assume it to have asymmetric terminator regions with R_p,1 = 0.1, R_p,2 = 0.09, and φ = −45°. Planetary positions (X, Y) were obtained using catwoman (which uses the exact same method as batman to calculate them) for 100 equally spaced time stamps between −0.5 and 0.5 days. We performed numerical simulations with n_p = 2500, 5000, 10,000, 20,000, and 40,000 (i.e., doubling the number of pixels on each side of our matrix) and found that the maximum flux changes roughly halved as well between each of those runs. These changes reached 4 ppm between n_p = 20,000 and 40,000, which we consider as our maximum error on the fluxes of our numerical scheme when selecting this latter number for n_p . Simulations using both our numerical (with n_p = 40,000) and semianalytical (through the catwoman library, with a maximum error set to 1 ppm ⁷ ) schemes are presented in Figure 2. As can be observed, the differences between both are very small; they reach peak differences of less than 7 ppm, most of which are explained by the errors defined by our numerical scheme.

Zoom In Zoom Out Reset image size

For all practical purposes, these limits give us confidence that our semianalytical framework works as expected for precisions that are better than current and near-future instruments such as JWST, which is expected to reach about 10 ppm light-curve precision (Greene et al. 2016). We note that the speed increase of the catwoman library in comparison to the numerical implementation is huge; catwoman takes a couple hundred milliseconds to generate a light curve in a 2.9 GHz Intel Core i9 processor. The numerical implementation takes tens of seconds to generate the same model, although the latter is a pure Python code, whereas the catwoman library is a mixture of Python and C. In general, in experiments made with this processor, catwoman takes about twice the time batman takes to generate a light curve. This is consistent with the fact that the catwoman code base is inherited from the batman one and goes to show that the analytical part of catwoman is as fast as batman's—only we perform it twice, one for each of the stacked semicircles.

In order to perform the simulations for JWST-like observations, we needed to calculate a typical cadence for observations to be taken by the observatory for time-series exposures. In this work, we focus on observations aiming to constrain the effect using NIRISS/SOSS, as this instrument allows us to obtain spectra all the way down to 0.6 μm through the combination of data from orders 1 (1–3 μm) and 2 (0.6–1 μm). Given that the largest limb asymmetries seem to be in the transition between optical and near-IR wavelengths (see, e.g., Kempton et al. 2017; Powell et al. 2019), we believe this will usually be the instrument of choice to characterize the effect for bright targets (with the alternative being, of course, NIRSpec/PRISM for fainter targets). Considering that the reset time for this instrument is relatively short (a couple of seconds), it sufficed for our work to know the typical integration time of a JWST observation with SOSS. For a solar-type, V = 11 star, according to the JWST Exposure Time Calculator ⁸ (Pontoppidan et al. 2016), "saturation" is attained at about 20–60 groups per integration for NIRISS/SOSS, which implies maximum integration times between 40 and 80 s per data point in the time series. We arbitrarily decided to use 20 s for the cadence of our simulations in order to simulate observations trying to target half-saturation values, which has been a typical strategy for HST observations. ⁹

In this case, we tried three different simulations in order to illustrate the impact of different assumptions in this one-transit, 20 s cadence case: (a) one in which everything but the radii of both sides of the exoplanet R_p,1 and R_p,2 is known, (b) one in which everything but the radii and the limb-darkening coefficients of the star is known, and, finally, (c) one in which everything but the radii, the limb-darkening coefficients of the star, and the time-of-transit center is known. The limb-darkening law that was assumed to generate and fit the light curves was a quadratic law with (u₁, u₂) = (0.3, 0.2), which are representative values for solar-type stars. In our juliet fits, we assumed the uninformative priors for two-parameter limb-darkening laws proposed by Kipping (2013) for the cases in which limb darkening was assumed to be unknown; for the time-of-transit center, a uniform prior with a width of 4.8 hr around the real predicted time-of-transit center was imposed. The results for these simulations are presented in Figure 3.

Zoom In Zoom Out Reset image size

As can be observed from the simulations, the asymmetric transit light curves should be detectable (i.e., ${\rm{\Delta }}\mathrm{ln}Z\gtrsim 2$) for morning-to-evening depth differences above around 25 ppm for a wide range of precisions, at least in white light (i.e., adding all of the flux over the entire wavelength range of a given instrument), where JWST observations should achieve tens of ppm precisions per point in the transit light curves, and as long as the ephemerides are well known and constrained. If they are not, however, as can be observed in the right panel of Figure 3(a), the actual detection of the effect becomes extremely challenging because, as was already noted by von Paris et al. (2016), changes in the time-of-transit center in a symmetric model can account for the asymmetry in the light curve. The changes in this timing are very small; only a couple of seconds of shifts in the time-of-transit center suffice to mimic the asymmetry in the transit light curves (see Section 4.3 for details). This implies that to detect this effect in white light, very precise timing is are needed in order to claim a detection.

It is important to note that although from the above results, the detection of the effect directly in the white-light light curves even with JWST-like precision seems relatively challenging to do with only one transit in the absence of precise timing constraints, the observatory has the advantage that it can perform spectrophotometry; thus, the effect can be detected through the wavelength dependence of the radii at each side of the terminator region, as was highlighted by Powell et al. (2019; see also Section 4). In particular, NIRISS/SOSS can produce extremely precise (tens of ppm) white-light transit light curves in orders 1 and 2, which can be used to claim a detection of the effect using these white-light transit light curves alone at much higher significance levels (and thus be sensitive to much lower morning-to-evening depth differences) than the ones shown here. Light curves like these, in addition, should provide very precise morning and evening depths. Figures 3(b) and (c) show an example light-curve fit on a 50 ppm-precision light curve, where the injected morning-to-evening depth difference was 270 ppm. As can be seen in Figure 3(b), both evening and morning depths are highly correlated but nonetheless provide precise constraints on each of about 16 ppm in this case, giving a retrieved morning-to-evening depth difference of Δδ = 257 ± 32 ppm, fully consistent with the input value. ¹⁰ Most of the information to constrain those depths come from ingress and egress, as is evident in the residuals (blue for catwoman, red for batman) of Figure 3(c). We provide a deeper understanding of the relation between light-curve precision and morning and evening depth precisions in Section 4.

Although the TESS mission has a significantly smaller aperture than JWST, the cadence and types of observations the mission does are excellent for the detection of asymmetries in transit light curves. The mission not only attains an exquisite precision, but it is also able to observe several transits of the same exoplanet, mitigating the problem we observed with only one transit in JWST-like observations like the ones simulated in the previous subsection. We performed the same simulations that we did for JWST but with a TESS-like cadence of 2 minutes, where we only consider observations on a 27 day period (i.e., one TESS sector). Interestingly, the three cases that we tried in the previous subsection (all known, all but limb darkening known, and all but limb darkening and the time-of-transit center known) all resulted in practically identical results; we show the one corresponding to the case in which all of the parameters are assumed to be known but the radius, limb-darkening coefficients, and time-of-transit center in Figure 4.

Zoom In Zoom Out Reset image size

As can be observed, the results are very similar to the ones of JWST. This is a combination of the fact that there is about a 27-fold increase in the number of transits, which helps with the sixfold increase on the cadence of the observations as compared to the JWST ones. The fact that there are more transits, in addition, helps with the problem JWST will face related to the ephemerides, where in our analysis, of course, there is an implicit assumption regarding no possible deviations from strict periodicity in the transit times. We reiterate, however, that our simple simulations in Section 3.1 did not consider the huge advantage that JWST has over TESS regarding the ability to measure wavelength-dependent transits, which should in turn break the degeneracy with the ephemerides as already suggested by Powell et al. (2019). We delve deeper into the benefits of wavelength-dependent JWST transit observations for detecting the limbs of exoplanets in Section 4.

Although few targets attain the precisions at which one might statistically distinguish between an asymmetric and a symmetric model directly from the transit light curves with TESS in one sector, the fact that many targets are observed by more than one sector makes this effect within reach of what TESS is currently able to detect. Targets in the JWST continuous viewing zone are particularly appealing to try to detect this effect. For example, WASP-62b (V = 10.3), for which the median per-point precision was 880 ppm during the prime mission in a 2 minute cadence, has been observed to date in more than 20 sectors, providing a combined per-point precision per transit of about at least $880/\sqrt{20}\approx 200$ ppm—a promising precision level to constrain the effect of asymmetric limbs if we assume that a simulation like the one in Figure 4 also applies to an exoplanet like WASP-62b. The recent HST study by Alam et al. (2021) also makes this target particularly appealing to detect this effect, as the atmospheric retrievals performed on its spectrum point toward it having a colder temperature (≈800 K) than what is expected from its equilibrium temperature (≈1400 K). This is one of the key hints that there might indeed be differences between the limbs of this exoplanet, as suggested by the work of MacDonald et al.(2020).

As a final test of the detectability of the effect, we explore whether our ignorance of the angle φ can impact it; we take our JWST-like simulation as a proxy for studying this, given the similarity in the shape of the detectability maps of JWST and TESS presented in Figures 3 and 4. To explore this, we use a transit light curve whose parameters are defined by the same ones as in the previous experiments, but in this case, we set φ = 45° as the ground truth and a uniform prior between −90° and 90° for the parameter in our fits. Our results for a JWST-like simulation (using the same cadence as in Section 3.1) are shown in Figure 5, where, in addition to the detectability map, we also show a portion of the posterior distribution of a simulation with the same properties as the one shown in Figure 3 for the case in which φ was fixed, i.e., a light-curve precision of 50 ppm and a morning-to-evening depth difference of 270 ppm.

Zoom In Zoom Out Reset image size

As can be seen, the detectability region (i.e., the medium and dark blue regions) of the plot has shifted by a small amount with respect to the one presented in Figure 3, implying that a slightly better light-curve precision is needed in order to detect the effect if the angle φ is not known a priori. In addition, given our large prior on φ, we actually detect two posterior modes (of which we show only one in Figure 5): one allowing for positive and one allowing for negative values of φ, each swapping the posterior distributions between the evening and morning limb. This was expected by construction; a model with φ and a given value of R_p,1 = a and R_p,2 = b is the same as a model with −φ with R_p,1 = b and R_p,2 = a. If we take the mode with positive values for φ, we note that in this case, the evening and morning depths themselves are much more uncertain. This gives rise to a retrieved morning-to-evening depth difference of ${250}_{-113}^{+374}$ ppm—again consistent but much less precise than the constraint assuming a known value of φ presented in Section 3.1. The constraint on φ itself is also not very precise; for this particular simulation, we obtain $\varphi ={32}_{-20}^{+33}$°, fully consistent with the input value but not very constraining to understand the true underlying projected planetary spin angle.

Before moving into the next section, it is important to reiterate that the precisions and detectability limits shown here were obtained for a very conservative, worst-case scenario system with a very small ingress/egress duration. As will be shown in Section 4.1, the odds of detecting the effect on systems that have better prospects for it (i.e., systems with longer ingress/egress durations) are in reality much higher. The lower limits we set here thus seem promising for the detection of the effect with current and near-future instrumentation.

While in Section 3, we presented lower limits on the statistical detection of the effect on transit light curves based on Bayesian evidence, an important aspect of interpreting transit light-curve fits with the semianalytic model presented in this work will involve constraining the actual measured transit depths of each side of the planet. This will be useful not only to extract transit spectra of the different limbs when using wavelength-dependent light curves such as the ones to be obtained by JWST but also to compute the maximum possible transit depth differences allowed by the data when analyzing broadband data such as that from missions like TESS.

An important detail to consider when extracting transit depths from asymmetric limbs is the fact that the observable quantities that are directly constrained by the data are the areas of each of the semicircles through the transit depths each of them produce. In symmetric models, where the limbs are assumed to be equal, the transit depth is simply $\delta =({R}_{p}^{2}/{R}_{* }^{2})$, the projected area of the planet over the projected area of the star. In the asymmetric case, however, the projected area of the semicircles is the quantities of interest, the transit depth of each limb being given by ${\delta }_{i}=(1/2)({R}_{p,i}^{2}/{R}_{* }^{2})$. This is what effectively defines the transit spectrum of each limb and is, in turn, what should be used to compare against theoretical transmission spectroscopy models.

Figure 6 shows how the precision of the transit depths of each limb depends on the light-curve precision, as well as the precision of the entire area of the planet, defined by the depth δ₁+δ₂, for the case of the exoplanet simulated in Section 3 with an ingress/egress duration of τ = 4.6 minutes (solid lines). As can be seen, the precision on the transit depth of the entire planet is always much smaller than the corresponding for both semicircles, but the relationship between the two is not simple, as the transit depths of the semicircles are highly correlated with each other. Indeed, the transit depth of the entire planet is constrained by the entire light curve, whereas the transit depths of each side of the planet (sampled by the semicircles in our model) are mostly constrained by ingress and egress. This implies, in turn, that this latter precision would, of course, increase on systems with larger ingress/egress durations, which in many cases might be the optimal ones to target in order to maximize the chances to unveil this effect.

Zoom In Zoom Out Reset image size

The simulations presented in Section 3 for the τ = 4.6 minute ingress/egress durations were performed to show lower limits on the detection of the effect, and even in those cases, the odds were very favorable given the current (e.g., TESS) and future (e.g. JWST) light-curve precisions and cadences. Hot Jupiters typically have longer ingress/egress durations, and some of the already characterized ones by missions like, e.g., HST, show good prospects for the detection of the effect as well. As an example, we repeat the JWST simulations in Section 3 for HAT-P-41b, which we select as is one of the most thoroughly characterized ultrahot Jupiters in transmission—all the way from near-UV, optical, and up to near-infrared wavelengths (Lewis et al. 2020; Wakeford et al. 2020; Sheppard et al. 2021). We tune the physical and orbital parameters of the system to the ones used in Wakeford et al. (2020), which imply a 23.9 minute ingress/egress duration. The cadence (53 s) and number of data points (500) for our simulations are set to the ones optimized by PandExo ¹¹ (Batalha et al. 2017) using NIRISS/SOSS as the instrument of choice such that SUBSTRIP256 does not saturate (which would be the setup of choice in order to obtain simultaneous spectroscopy in the near-infrared and optical through orders 1 and 2 for this target). ¹² For consistency, we set the limb-darkening coefficients to the average ones on order 1 of NIRISS/SOSS (but we note that these do not impact the overall precision and detectability, as was already shown in Section 3). The resulting precisions of this experiment are presented in Figure 6 as dashed lines. As can be observed, the precision change on the transit depths of each of the limbs is significant and ranges from a 60% to 70% improvement in it. The precision change in the transit depth of the entire planet, however, is much smaller (and driven mainly by the difference in the absolute transit depths and transit durations), which acts as a baseline in showing quantitatively how the prospects of detecting asymmetric limb differences are very sensitive to the ingress/egress duration.

One might argue that in Figure 6, one could observe two transits of the exoplanet with 4.6 minutes ingress/egress duration in order to match the signal-to-noise ratio of the exoplanet with the 23.9 minute ingress/egress duration. Although this would be true if the observatory only targeted the transit event, in practice, there are observational overheads (like, e.g., pre- and posttransit baselines and overall observatory overheads beyond clock time on target) that have to be included in that reasoning. For instance, the time recommended in the JWST documentation ¹³ one should spend in a target during a transit is given by the dwell Equation, which reads

$$\begin{eqnarray}&&{T}_{{\rm{dwell}}}=0.75+{\rm{MAX}}(2,{T}_{14})+{T}_{14}+1\ {\rm{hr}},\end{eqnarray} \tag{ 3 }$$

where T₁₄ is the transit duration in hours. The exoplanet with a 4.6 minute ingress/egress duration has a 1 hr total transit duration, which gives T_dwell = 4.75. The exoplanet with a 23.9 minute ingress/egress duration has a 3.6 hr transit duration, T_dwell = 8.95. Two transits of the 1 hr transit duration target would imply a requested time of 9.5 hr, which is at least half an hour more expensive than the 3.6 hr transit duration target—all this without considering extra observatory overheads. The conclusion, thus, is that the efficiency of the time and targets to be requested to detect the effect have to be studied on a case-by-case basis.

In addition to the above, it is also important to note that the posterior distributions of the limbs show typically nonnegligible correlations (see, e.g., Figures 3 and 5). This is important to consider because this covariance carries extra information that could become important when performing inference on the limbs (through, e.g., atmospheric retrievals). In our experiments, we observed that the shape of the posterior distribution (measured through the correlation, i.e., the covariance divided by the standard deviations of the marginal posterior distribution of each limb) remained fairly constant across the parameter space covered in our work; it simply shrinks as better precisions are achieved. We also found that the posterior distribution of the limbs is very well approximated in the case of a known angle ϕ by a 2D Gaussian distribution, which is characterized not only by the mean and standard deviation of the marginal posterior distributions but also by the covariance between the limbs. We provide a practical example of how to use this information in a retrieval framework in the next subsection.

In order to showcase the importance of directly using transit light curves to constrain the limbs of exoplanets and demonstrate that this indeed might be the most efficient avenue to constrain the transit spectrum of the limbs, we use HAT-P-41b as a case study, using both real HST data and simulated JWST observations.

We first analyze the HST data presented in Lewis et al. (2020), Wakeford et al. (2020), and Sheppard et al. (2021). In particular, we take the HST/WFC3 (1.1–1.6 μm), STIS G430L (0.3–0.5 μm), and G750L (0.5–0.9 μm) data from Sheppard et al. (2021) and the WFC3 UVIS G280 (0.2–0.8 μm) data from Wakeford et al. (2020). In order to perform the atmospheric retrievals, we use the framework developed in MacDonald et al. (2020), namely, one in which the transmission spectrum is fit with a model using two different atmospheric structures, each of which represents one of the limbs. To this end, we use and slightly modify the "chemically consistent" CHIMERA ¹⁴ atmospheric retrieval framework of Line et al. (2013), such that two transmission spectra are modeled at each iteration of the algorithm. Each of those has a different temperature and cloud structure, as well as different C/O ratios; the temperature at the top of the atmosphere is forced to always be colder for one limb than the other in order to simulate energy redistribution processes that might be happening between the limbs. Chemical equilibrium is assumed for each limb. Once those are computed, each model transit depth is multiplied by 1/2 in order to compute the spectra of each of the limbs, δ_i . These are then added together to form the "combined" transit depth that we compare against the HST transit spectrum. We decide to leave the metallicity, 10 bar radius, and overall temperature–pressure profile shape (other than the temperature at the top at the atmosphere) as common parameters between the limbs, as we assume the HST data would not be sensitive to those parameters. We also only use the data at wavelengths longer than 0.35 μm, as the publicly available opacities within CHIMERA only go down to this particular wavelength range.The retrieval is performed using nested sampling with the the pymultinest library (Buchner et al. 2014), which makes use of the MultiNest (Feroz et al. 2009) algorithm. The full set of priors used for our CHIMERA atmospheric retrievals in what follows are presented in Table 2.

Table 2. Priors and Parameters Used for Our Section 4.2 CHIMERA Atmospheric Retrievals

Parameter	Prior	Comment
Parameters Individual to Each Limb
Tirr (K)	U(1000, 2500)	Stellar input at
		top of atmosphere
log10C/O	U(−2, 0.3)	Log C/O ratio
log10 Kzz (cm2 s−1)	U(5, 11)	Log-eddy diffusion
		coefficient
fsed	U(0.5, 5)	Sedimentation efficiency
${\mathrm{log}}_{10}{P}_{{\rm{base}}}$ (bars)	U(−6, 1.5)	Cloud base pressure
${\mathrm{log}}_{10}{f}_{{\rm{cond}}}$	U(−15, −2)	Cloud condensate mixing
		ratio at cloud base
Parameters Common to Both Limbs
[M/H]	U(−2, 3)	Atmospheric metallicity
fR	U(0.5, 1.5)	Multiplicative factor to
		10 bar "fiducial" radius
${\mathrm{log}}_{10}{\kappa }_{{\rm{IR}}}$ (cm2 g−1)	U(−3, 0)	T/P profile IR opacity
${\mathrm{log}}_{10}{\gamma }_{1}$	U(−3, 0)	Log visual-to-IR opacity

Note. These are composed of two limbs, one forced to be colder than the other. Here U(a, b) stands for a uniform distribution between a and b. For detailed discussions of each parameter, see Line et al. (2013) for the general framework and Mai & Line (2019) for the implementation of the cloud model used in this work, which is that of Ackerman & Marley (2001).

Download table as:  ASCIITypeset image

The best-fit retrieval to the HST data using this two-limb retrieval framework is presented in the top panels of Figure 7; the corresponding retrieved temperature–pressure profile is also presented in the left panel of Figure 8. Posterior credibility intervals (C.I.s) for each parameter fit in the retrieval are presented in the second column of Table 3. As can be seen in Figure 7, while the retrieved transit spectrum follows the data fairly well, the individual retrieved limb spectra show a wide range of possible solutions, suggesting that they are fairly unconstrained by the data. Indeed, the overall posterior values for all parameters are largely indistinguishable between the limbs. For example, the temperatures at the top of the atmosphere are ${1210}_{-132}^{+227}$ and ${1522}_{-241}^{+519}$ K for the "cold" and "hot" limbs, respectively, which are completely consistent with each other. The rest of the parameters, while fairly unconstrained by the data, all point toward a cloudy nature in HAT-P-41b's atmosphere; our solution suggests a relatively high-altitude cloud-deck solution whose base is at about 10 mbars—a picture in good agreement with the forward modeling performed in Lewis et al. (2020). While the retrievals in that work preferred a near-UV opacity source shaping the spectrum instead of clouds, most of these retrievals were done by fitting the atomic and molecular abundances directly, without chemical constraints between them. The only chemical equilibrium retrieval performed in it needed a large temperature (>2500 K) to be able to build up the necessary H⁻ opacity to explain the data better. The chemical equilibrium retrieval performed here, while including H⁻ opacity, is much more constrained in temperature, as we did not allow the possibility for the temperature at the limbs to exceed 2500 K. Nonetheless, while we see a possible second mode in the temperature posterior in Figure 7, where the hot limb can reach temperatures of about 2300 K, we see no buildup of posterior density above this range. While our retrieval could thus be interpreted as an alternative solution to the observed HST data, we do not discuss this at length here, as the objective of this experiment was to only obtain a possible solution for a two-limb model given the HST data, in order to showcase the importance of constraining limb spectra with upcoming high-precision facilities, as we will illustrate and quantify below.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3. Posterior C.I.s of the Two-limb Retrieval Made on Real HST Data and Simulated JWST/NIRISS Data of HAT-P-41b

Parameter	Posterior C.I.	Posterior C.I.	Posterior C.I.
	(Transit Spectrum, HST)	(Transit Spectrum, JWST a )	(Limb Spectra, JWST a )
Parameters for the "Cold" Limb
Tirr (K)	${1210}_{-132}^{+227}$	${1197}_{-115}^{+130}$	${1237}_{-70}^{+59}$
${\mathrm{log}}_{10}C/O$	$-{0.85}_{-0.65}^{+0.56}$	$-{0.98}_{-0.57}^{+0.42}$	$-{0.99}_{-0.56}^{+0.42}$
${\mathrm{log}}_{10}{K}_{{zz}}$ (cm2 s−1)	${7.9}_{-1.6}^{+1.7}$	${8.0}_{-1.9}^{+1.9}$	${7.7}_{-1.7}^{+1.9}$
fsed	${3.4}_{-1.7}^{+1.6}$	${3.4}_{-1.7}^{+1.6}$	${3.3}_{-1.7}^{+1.8}$
${\mathrm{log}}_{10}{P}_{{\rm{base}}}$ (bars)	$-{2.1}_{-1.9}^{+1.9}$	$-{1.9}_{-2.5}^{+2.1}$	$-{1.8}_{-2.6}^{+2.1}$
${\mathrm{log}}_{10}{f}_{{\rm{cond}}}$	$-{7.0}_{-4.5}^{+3.1}$	$-{9.3}_{-3.5}^{+3.6}$	$-{9.2}_{-3.5}^{+3.9}$
Parameters for the "Hot" Limb
Tirr (K)	${1522}_{-241}^{+519}$	${1592}_{-108}^{+94}$	${1567}_{-76}^{+66}$
${\mathrm{log}}_{10}C/O$	$-{0.68}_{-0.78}^{+0.50}$	$-{1.16}_{-0.52}^{+0.49}$	$-{1.23}_{-0.48}^{+0.51}$
${\mathrm{log}}_{10}{K}_{{zz}}$ (cm2 s−1)	${8.0}_{-1.7}^{+1.7}$	${7.8}_{-1.8}^{+2.0}$	${7.8}_{-1.8}^{+2.1}$
fsed	${3.1}_{-1.5}^{+1.6}$	${3.3}_{-1.6}^{+1.7}$	${3.3}_{-1.7}^{+1.7}$
${\mathrm{log}}_{10}{P}_{{\rm{base}}}$ (bars)	$-{2.5}_{-1.9}^{+2.0}$	$-{1.9}_{-2.5}^{+2.1}$	$-{1.9}_{-2.5}^{+2.1}$
${\mathrm{log}}_{10}{f}_{{\rm{cond}}}$	$-{8.0}_{-4.0}^{+3.8}$	$-{9.1}_{-3.6}^{+3.7}$	$-{9}_{-3.7}^{+3.8}$
Parameters Common to Both Limbs
[M/H]	${2.06}_{-0.37}^{+0.24}$	${2.02}_{-0.041}^{+0.036}$	${2.022}_{-0.039}^{+0.037}$
fR	${1.0276}_{-0.0087}^{+0.0072}$	${1.0273}_{-0.0027}^{+0.0017}$	${1.0276}_{-0.0022}^{+0.0014}$
${\mathrm{log}}_{10}{\kappa }_{{\rm{IR}}}$ (cm2 g−1)	$-{2.23}_{-0.50}^{+0.51}$	$-{1.92}_{-0.61}^{+0.75}$	$-{2.32}_{-0.43}^{+0.57}$
${\mathrm{log}}_{10}{\gamma }_{1}$	$-{1.57}_{-0.81}^{+0.72}$	$-{1.74}_{-0.75}^{+0.77}$	$-{1.76}_{-0.80}^{+0.95}$

Note. Here C.I. corresponds to the medians and the 68% credibility bands around them. For details on the definition for each parameter and priors, see Table 2.

^a The JWST simulations had as underlying true values the medians from our HST retrievals.

Download table as:  ASCIITypeset image

The next experiment we ran consisted of taking the median of all the parameters found from our two-limb retrieval of the HAT-P-41b HST data—i.e., the medians presented in the second column of Table 3—and use the corresponding limb models implied by them to simulate a single JWST NIRISS/SOSS transit observation using the catwoman transit model introduced in this work. As in the last subsection, we used PandExo to estimate the per-integration light-curve precision for this transit observation, with which we simulated 608 53 s integrations, i.e., an 8.95 hr exposure, the recommended time on target for JWST observations, obtained using Equation (3). To generate the transit models at each wavelength, we used the same transit parameters used by Wakeford et al. (2020), with the exception of the radius for each limb, for which we use the model described above. As for limb darkening, we used the ExoCTK (Bourque et al. 2021) limb-darkening tool ¹⁵ to compute the limb-darkening coefficients of the quadratic law for orders 1 and 2 of NIRISS/SOSS.

We tackled the analysis of these simulated light curves in two different ways in order to showcase why performing inference on the limb spectra is the best way to optimally extract information about the limbs of exoplanets. The first consisted of extracting transit spectra from them by simply fitting these light curves with batman light-curve models with the same procedure explained in Section 3.1. Then, we performed two-limb retrievals on this transit spectrum using the same modified CHIMERA atmospheric retrievals described above and performed on the HST data. The results of that experiment are presented in the middle panel of Figure 7, the retrieved pressure–temperature profiles are also shown in the middle panel of Figure 8, and posterior C.I.s for each fitted parameter in the retrieval are presented in the third column of Table 3. As can be seen, these two-limb retrievals much better constrain the different limbs than the HST data but still have a fairly large uncertainty, such that the retrieved limb spectra overlap between the limbs at certain wavelengths. While the posteriors for most parameters are fairly wide (mainly due to the cloudy nature of this particular exoplanet), the temperature posterior distribution is much better constrained than the corresponding one for the HST data; they have uncertainties that are about onefold to twofold smaller for the cold limb and about twofold to fivefold smaller for the hot limb.

Next, we performed the transit fitting using our catwoman model, again using the same procedures as the ones outlined in Section 3.1; this allowed us to extract limb spectra directly from the transit light curves. In order to perform retrievals on these limb spectra, however, we had to slightly modify the log-likelihood of our retrievals because, as explained in the previous subsection, for each wavelength bin, the two limb depths are strongly correlated. We characterize the posterior distributions of the limb depths at each wavelength bin in our simulations as multivariate Gaussians, which are described by covariance matrices Σ_w and mean vectors μ _w = [δ_1,w , δ_2,w ] for each wavelength w. The covariance matrices were estimated in the case of our fits by computing the sample covariance matrix of the posterior samples of the limb depths at each wavelength bin. With this, the log-likelihood we use in our retrievals is given by

$$\begin{eqnarray*}&&\mathrm{ln}{ \mathcal L }=-\displaystyle \frac{1}{2}\sum _{w}2\mathrm{ln}2\pi +\mathrm{ln}| {{\rm{\Sigma }}}_{w}| +{{\boldsymbol{r}}}_{w}^{T}{{\rm{\Sigma }}}_{w}^{-1}{{\boldsymbol{r}}}_{w},\end{eqnarray*}$$

where r _w = μ _w − m _w is the vector of the residuals, with m _w being the vector containing the model limb spectra at the given wavelength bin w.

The results of performing the retrieval on the limb spectra directly using the framework described above are presented in the bottom panel of Figure 7, with the corresponding retrieved pressure–temperature profiles shown in the right column of Figure 8 and posterior C.I.s for each fitted parameter in the retrieval presented in the fourth column of Table 3. As can be observed, the constraints on the limb spectra are much better than the ones obtained from fitting the transit spectra. The corresponding limb spectral models are much better differentiated by our retrievals, allowing us to tightly constrain any features that might arise in them. The temperature posteriors for the limbs are also much better constrained; in particular, there is a two-to-fourfold precision increase on the temperature of the cold limb and a three-to-eightfold precision increase in the hot limb when compared against the HST data. In addition, the correlation between the retrieved limb temperatures is much smaller for this limb retrieval, as can be observed on the posterior distribution presented in Figure 7. This tighter constraint on the temperature can also be visually observed in Figure 8, where the much tighter constraint in the temperature–pressure profile of each limb is illustrated as well. These results showcase that extracting and analyzing limb spectra directly allows us to much more tightly constrain the properties of the limbs than methods that rely on extracting this information from "classic" transit spectra. ¹⁶

Aside from impacting directly on inferences made with atmospheric retrievals, as has been shown above, we would like to highlight that extracting limb spectra directly from transit light curves has an additional, unique benefit: it opens up the discovery space for atmospheric features that might be individual to each of the limbs. This, in turn, has the benefit of allowing a direct comparison against global circulation models (GCMs). The GCM modeling assumptions and implementation details have been actively driven by observations in the past few years. For instance, until recently, clouds and hazes were typically added in a postprocessing stage, not including feedback from these components in the modeling. Recent works, such as those of Roman & Rauscher (2019) and Parmentier et al. (2021), have managed to implement clouds in their GCMs directly, showing, in turn, that this is critical to understanding the overall cloud structure itself and heat redistribution in the planets under study, which has provided a much better match to observed phase-curve observations. We believe the study of limb spectra can open up a similar set of insights providing key advancements in the overall modeling of exoplanetary atmospheres.

From an observational perspective, it is interesting to consider that observations aiming to extract these limb spectra can be much cheaper than, e.g., phase curves. Whereas a single NIRISS/SOSS transit for HAT-P-41b, for instance, would amount to a total science time for a JWST proposal of about 8.95 hr (using Equation (3)), a full phase curve for this particular exoplanet requires over 65 hr. In this sense, extraction of limb spectra could serve as a good diagnostic as to what to expect in a phase-curve observation before performing these expensive observations. Since for fixed planetary and stellar radii, the ingress/egress duration increases with the square root of the semimajor axis (i.e., $\tau \propto \sqrt{a}$), this technique for detecting limb asymmetries might, in turn, be an excellent alternative avenue for studying mornings and evenings of longer-period planets where phase-curve signals are too small to be detectable in a reasonable amount of time.

In Section 3, we discussed how, as predicted by von Paris et al. (2016), small changes in the time-of-transit center can give rise to equally good fits on symmetric models (such as the ones assumed by the batman package, for instance), even if the data are truly arising from an asymmetric transit model (although this limitation can be bypassed if the aim is to retrieve limb spectra by simultaneously fitting transit light curves at different wavelengths; see, e.g., Powell et al. 2019). It is important to note that this, in turn, can give rise to biases in transit times if a symmetric model is used when the data in fact come from an asymmetric model such as the one modeled by the catwoman library. In our simulations, these can give rise to timing offsets of up to 5 s, which is, in turn, within the timing precision that the TESS mission is able to reach and will be for sure within reach of the JWST mission. Thus, care must be taken when searching for small (second-level) timing offsets in these precise transit light curves in the search for, e.g., transit timing variations.

It is important to note that throughout this work, we have assumed that the only alternative model to that of limb asymmetries is that of a symmetric limb in order to explain transit light-curve asymmetries. However, there are other competing effects that might give rise to asymmetric light curves as well. For instance, known stellar effects such as gravity darkening (see, e.g., Ahlers et al. 2020, and references therein) and yet-to-be-uncovered effects/properties such as exoplanetary rings (see, e.g., Rein & Ofir 2019, and references therein) can also give rise to asymmetric light curves. Performing a detailed model comparison study between these effects and the one studied here is outside the scope of this work, but we warn researchers that proper care must be taken when aiming to claim the detection of asymmetric limbs in light of these possibilities. While effects like, e.g., gravity darkening are most likely known at good enough precisions in order to understand when a given transit light curve might be asymmetric due to this effect, or at the very least to put limits on asymmetries generated by it, known unknowns such as exoplanetary rings might be more complicated to rule out. Perhaps the easiest way to constrain this would be through the wavelength dependence of these asymmetries, which we hypothesize should be markedly different in the case of exoplanetary rings and those produced by opacity sources in an exoplanetary atmosphere. Still, it is important to be mindful of these alternative hypotheses when analyzing data on the search for these light-curve asymmetries.

In addition to the above, the very validity of the framework developed here—i.e., fitting transit light curves with a model of two "stacked," rotated semicircles—remains to be put to the test with real data and has plenty of room for improvement as we increase the precision of our measurements. For instance, our model assumes a sharp discontinuity at the poles, whereas some GCMs actually predict smooth transitions at them (see, e.g., Pluriel et al. 2020, and references therein), even suggesting that morning and evening terminators might be themselves asymmetric. While a path in the right direction, our model is just an approximation of reality by construction, as those simplifications were the ones that allowed us to create a modeling scheme that is fast and efficient to apply to real transit light-curve data. We expect that the detection of the signatures of mornings and evenings on actual, precise transit light curves from TESS or JWST could indeed motivate more flexible and accurate models for their shapes guided by GCM modeling (or the data themselves) that could expand on the simple modeling scheme discussed in this work. While these shapes might be complex enough that simple geometrical arguments like the ones made in this work would most likely not be easy to make, making the problems hard to parameterize, we are hopeful that good ideas might flourish in the exoplanet community to make this happen. There currently exists a continuum of light-curve analysis methods ranging from models that have few parameters but are constrained, like ours, to models that are very flexible but have a larger number of parameters, like the shadow imaging technique presented in Sandford & Kipping (2019). We believe expanding our model in the direction discussed above lies in between those methodologies and is thus bounded and therefore approachable. Developing this idea further, however, is outside the scope of this work.

Although the case of calculating the intersection area between two circles can be obtained via elemental trigonometry, the problem of calculating the intersection area between a circle and a (rotated) semicircle is not, in general, as straightforward.

We first note that the problem of finding the intersection area of a circle of radius R and a semicircle of radius R_S rotated by an angle θ with respect to the line that joins the centers of length d is the same problem as the intersection area of a semicircle and a circle rotated by an angle θ with respect to the line that joins the centers. This symmetry argument allows us to put the horizontal axis of this problem in the base of the semicircle, simplifying the notation of the problem. Without loss of generality, we put the origin in the center of the base of the semicircle. This transformed geometry of the problem is shown in Figure 9; here the white hatched area inside the semicircle is the area of interest (i.e., the one that leads to A_S (R, R_S , d, θ)).

Zoom In Zoom Out Reset image size

As is evident in Figure 9, there are three different cases (a, b, and c) we have to take care of in order to find a general formula for A_S (R, R_S , d, θ).

• 1.  
  Case (a), divided into subcases (a1), (a2), and (a3), deals with the problem in which the circle is rotated such that it lies above the semicircle. If we identify the coordinates of the center of the circle as $({x}_{0},{y}_{0})=(-d\cos \theta ,d\sin \theta )$, case (a) deals with the problem in which θ > 0 and, thus, d sin θ > 0. Here the intersection points between the circle and the semicircle have coordinates (a₁, b₁) and (a₂, b₂). The geometry depicted in Figure 9 for this case implies that b₁ = 0. This is because for b₁ > 0, the problem is the same as the intersection of two circles (one of radius R and another of radius R_S ), which has a known analytical solution (see, e.g., Kreidberg 2015). Here the area of interest, A_S , is given by the area of the semicircle ($\pi {R}_{S}^{2}/2$) minus A₁ + A₂ for case (a1) and by A₁ + A₂ + A₃ for cases (a2) and (a3). The different subcases depend, in turn, on the location of the intersection points and the position of (a₃, b₃), the position of the maximum extension of the circle on the x-axis.
• 2.  
  Case (b) deals with the problem in which the circle is rotated such that it lies below the semicircle, i.e., where θ < 0 and, thus, $d\sin \theta \lt 0$. In addition, this case handles only problems in which one intersection of the circle with the semicircle is in its base and the other is on the upper part of the semicircle. Once again, the intersection points between the circle and the semicircle have coordinates (a₁, b₁) and (a₂, b₂). In this case, however, b₂ = 0; the cases in which b₂ ≠ 0 (i.e., when the rightmost intersection is on the upper part of the semicircle) and b₂ = 0 and b₁ = 0 (i.e., in which the leftmost intersection is also on the base of the semicircle) are taken care of by case (c). The area of interest for case (b) is A_S = A₁ + A₂.
• 3.  
  Finally, case (c), divided into subcases (c1), (c2), and (c3), deals with the problem in which the circle is either rotated above or below the semicircle, but where there are two intersections with both either in the base (c1) or in the upper part of the semicircle (c2) or four intersection points (c3) between the circle and the semicircle. In this case, the area of interest, A_S , can be calculated directly via basic trigonometry; thus, here we do not identify the intersection points by coordinates but rather by the red points in order to guide the reader.

Cases (a), (b), and (c) defined above will all be calculated assuming that the center of the circle is to the left of the semicircle. The reason for doing this is that the problem has reflective symmetry with respect to the line that goes through the center of the semicircle and that is perpendicular to its base. As such, for 0 < θ < π/2, we have that A_S (R, R_S , d, θ) = A_S (R, R_S , d, π−θ), whereas for −π/2 < θ < 0, A_S (R, R_S , d, θ) = A_S (R, R_S , d, θ−π). In what follows, we solve each of the cases separately.

A.1.1. Case (a)

Before looking at the integrals that lead to the intersection area between the circle and the semicircle in this case, let us obtain the expressions for the intersection points (a₁, b₁) and (a₂, b₂). To do this, we first note that the point where the circle intersects with the line y = 0 (which is the line that passes through the base of the semicircle) is given by

$$\begin{eqnarray}&&{x}_{\mathrm{int}}=-d\cos \theta \pm \sqrt{{R}^{2}-{d}^{2}{\sin }^{2}\theta }.\end{eqnarray} \tag{ A1 }$$

If the discriminant of this expression is negative (or zero), which in this case implies that $| d\sin \theta | \geqslant R$, then there is no (or one) intersection of the circle with the x-axis. In this case, the intersectional area can either be zero if d ≥ R+R_S or equal to solving the problem of the intersection between two circles of radius R and R_S in the case in which d < R+R_S . If the latter case is true, in turn, the intersectional area reduces to πR² when d+R < R_S (i.e., when the circle is inside the semicircle).

If $| d\sin \theta | \lt R$, then we have two solutions for the intersection points of the circle and the x-axis, x_int. The cases in which both solutions are inside the semicircle, i.e., when ∣x_int∣ ≤ R_S , will be handled in case (c), as depicted in Figure 9. If both solutions are outside the semicircle, i.e., when ∣x_int∣ > R_S , there are two possibilities. If d ≥ R+R_S , the intersection area is zero; in this case, the circle is always to the left of the semicircle. If, on the other hand, d < R+R_S , then because the intersections with the x-axis happen outside the semicircle, there are two options: either the circle intersects twice with the upper part of the semicircle (a case that will be solved in case (c2)) or the circle covers all of the semicircle, in which case the intersectional area is $\pi {R}_{S}^{2}/2$. If θ > 0, then this leaves us, finally, with the problem we will solve for case (a) (if this is not true, then case (b) will apply), in which x_int has both one solution outside and another one inside the semicircle. Because here we are dealing with the cases in which the center of the circle is always to the left of the center of the semicircle, this implies that the intersection point inside the semicircle will always be the rightmost intersection point, i.e., the solution of x_int with the positive sign. In the notation of Figure 9, this gives the intersection points

$$\begin{eqnarray*}\begin{array}{rcl}{a}_{1} & = & -d\cos \theta +\sqrt{{R}^{2}-{d}^{2}{\sin }^{2}\theta },\\ {b}_{1} & = & 0.\end{array}\end{eqnarray*}$$

The intersection point between the upper part of the semicircle and the circle, (a₂, b₂), is obtained by simply equating the equation of the circle (${(x+d\cos \theta )}^{2}+{(y-d\sin \theta )}^{2}\,=\,{R}^{2}$) with the equation of the semicircle, taking it as a full circle (${x}^{2}+{y}^{2}={R}_{S}^{2}$) to begin with. This yields

$$\begin{eqnarray}&&{b}_{2}=-A\sin \theta +\cos \theta \sqrt{{R}_{S}^{2}-{A}^{2}},\end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray}&&{a}_{2}={b}_{2}\tan \theta +\displaystyle \frac{A}{\cos \theta },\end{eqnarray} \tag{ A3 }$$

where

$$\begin{eqnarray*}&&A=\displaystyle \frac{{R}^{2}-{R}_{S}^{2}-{d}^{2}}{2d}.\end{eqnarray*}$$

As we are not dealing with two intersecting circles but rather an intersecting circle and semicircle, we have chosen the largest b₂ that will give rise to the a₂, and therefore the point of intersection, that is on the semicircle.

Finally, an important set of coordinates to define are the ones for (a₃, b₃). As illustrated in Figure 9, these are the coordinates of the maximum value attained on the x-axis by the circle. The coordinates for this point are, evidently,

$$\begin{eqnarray*}\begin{array}{rcl}{a}_{3} & = & -d\cos \theta +R,\\ {b}_{3} & = & d\sin \theta .\end{array}\end{eqnarray*}$$

First, we take on case (a1). This case occurs when the conditions for case (a) are met and the point (a₃, b₃) is outside the semicircle, which in turn implies in this case that b₃ ≥ b₂. To solve it, the strategy to obtain A_S is to compute analytic solutions to the areas A₁ and A₂ depicted in Figure 9 and then subtract these to $\pi {R}_{S}^{2}/2$. First, A₁ is simply the area under the curve of the circle from x = a₁ to a₂. Because in this case, a₁ < a₂ < a₃, we are going to integrate the lower part of the circle; this implies that the equation of the (in this case, semi) circle is simply

$$\begin{eqnarray*}&&y=-\sqrt{{R}^{2}-{\left(x+d\cos \theta \right)}^{2}}+d\sin | \theta | .\end{eqnarray*}$$

Integrating this from x = a₁ to a₂ yields

$$\begin{eqnarray}&&{A}_{1}=\displaystyle \frac{{R}^{2}}{2}{\rm{\Delta }}f+d\sin | \theta | ({a}_{2}-{a}_{1}),\end{eqnarray} \tag{ A4 }$$

where Δf = f(a₁)−f(a₂), with

$$\begin{eqnarray*}\begin{array}{rcl}f(x) & = & \left(\displaystyle \frac{x+d\cos \theta }{{R}^{2}}\right)\sqrt{{R}^{2}-{\left(x+d\cos \theta \right)}^{2}}\\ & & +\arcsin \left(\displaystyle \frac{x+d\cos \theta }{R}\right).\end{array}\end{eqnarray*}$$

Next, we work on obtaining area A₂. This is simply the area under the curve of the semicircle, whose equation is $y=\sqrt{{R}_{S}^{2}-{x}^{2}}$. Integrating this from x = a₂ to R_S yields

$$\begin{eqnarray}&&{A}_{2}=\displaystyle \frac{{R}_{S}^{2}}{2}\left[\displaystyle \frac{\pi }{2}-h({a}_{2})\right],\end{eqnarray} \tag{ A5 }$$

where

$$\begin{eqnarray*}&&h(x)=\arcsin \left(\displaystyle \frac{x}{{R}_{S}}\right)+\displaystyle \frac{x}{{R}_{S}}\sqrt{1-\displaystyle \frac{{x}^{2}}{{R}_{S}^{2}}}.\end{eqnarray*}$$

Thus, using the definitions for A₁ given in Equation (A4) and A₂ given in Equation (A5), the area A_S is given in this case by

$$\begin{eqnarray*}&&\boxed{{A}_{S}=\displaystyle \frac{\pi {R}_{S}^{2}}{2}-\displaystyle \frac{{R}^{2}}{2}{\rm{\Delta }}f-d\sin | \theta | ({a}_{2}-{a}_{1})-\displaystyle \frac{{R}_{S}^{2}}{2}\left[\displaystyle \frac{\pi }{2}-h({a}_{2})\right]}.\end{eqnarray*}$$

Now, we take on case (a2). In this case b₃ < b₂; however, a₂ > a₁. In this case, the area of interest is A_S = A₁+A₂+A₃, as depicted in Figure 9. First, area A₁ in this case is the area of the semicircle from x = −R_S to a₁. Integrating once again the equation of the semicircle ($y=\sqrt{{R}_{S}^{2}-{x}^{2}}$) in this range, one obtains

$$\begin{eqnarray}&&{A}_{1}=\displaystyle \frac{{R}_{S}^{2}}{2}\left[h({a}_{1})+\displaystyle \frac{\pi }{2}\right].\end{eqnarray} \tag{ A6 }$$

Area A₂ in this case can be calculated as the area under the same semicircle between x = a₁ and a₂ minus b₂(a₂−a₁)/2, which is the area of the triangle formed between the points (a₁, b₁), (a₂, b₂), and (a₂,0). Integrating the semicircle between x = a₁ and a₂ and subtracting b₂(a₂−a₁)/2, we obtain

$$\begin{eqnarray}&&{A}_{2}=\displaystyle \frac{{R}_{S}^{2}}{2}\left[h({a}_{2})-h({a}_{1})\right]-\displaystyle \frac{{b}_{2}({a}_{2}-{a}_{1})}{2}.\end{eqnarray} \tag{ A7 }$$

Finally, area A₃ reduces to obtaining the segment of a circle generated by drawing a chord between points (a₁, b₁) and (a₂, b₂). To this end, we ought to know the angle α (in radians) these points make with respect to the center of the circle. This can be easily obtained by the law of cosines to give

$$\begin{eqnarray*}&&\alpha =\arccos \left(1-\displaystyle \frac{{\left({a}_{2}-{a}_{1}\right)}^{2}+{b}_{2}^{2}}{2{R}^{2}}\right).\end{eqnarray*}$$

With this, the area of segment A₃ is thus simply

$$\begin{eqnarray}&&{A}_{3}=\displaystyle \frac{{R}^{2}}{2}\left(\alpha -\sin \alpha \right).\end{eqnarray} \tag{ A8 }$$

Finally, then, using the definition for A₁ in Equation (A6), A₂ in Equation (A7), and A₃ in Equation (A8), we get in this case

$$\begin{eqnarray*}&&\boxed{{A}_{S}=\displaystyle \frac{{R}_{S}^{2}}{2}\left[h({a}_{2})+\displaystyle \frac{\pi }{2}\right]-\displaystyle \frac{{b}_{2}({a}_{2}-{a}_{1})}{2}+\displaystyle \frac{{R}^{2}}{2}\left(\alpha -\sin \alpha \right)}.\end{eqnarray*}$$

Finally, we solve case (a3). In this case, again, b₃ < b₂; however, a₁ > a₂. Here Equation (A6) also applies for A₁, but the upper limit of the integral is in this case a₂ instead of a₁. This implies that in this case,

$$\begin{eqnarray}&&{A}_{1}=\displaystyle \frac{{R}_{S}^{2}}{2}\left[h({a}_{2})+\displaystyle \frac{\pi }{2}\right].\end{eqnarray} \tag{ A9 }$$

To obtain area A₂ in this case, we note that here this is simply the area of the triangle formed by the points with coordinates (a₁, b₁), (a₂, b₂), and (a₂, 0). In this case, thus,

$$\begin{eqnarray}&&{A}_{2}=\displaystyle \frac{{b}_{2}({a}_{1}-{a}_{2})}{2}.\end{eqnarray} \tag{ A10 }$$

Finally, to obtain A₃, we use Equation (A8), which also applies for this case. Using, then, the definition for A₁ in Equation (A9), A₂ in Equation (A10), and A₃ in Equation (A8), we get in this case

$$\begin{eqnarray*}&&\boxed{{A}_{S}=\displaystyle \frac{{R}_{S}^{2}}{2}\left[h({a}_{2})+\displaystyle \frac{\pi }{2}\right]+\displaystyle \frac{{b}_{2}({a}_{1}-{a}_{2})}{2}+\displaystyle \frac{{R}^{2}}{2}\left(\alpha -\sin \alpha \right)}.\end{eqnarray*}$$

A.1.2. Case (b)

Case (b) is similar in many ways to case (a), with the only difference being that now the coordinates of the center of the circle change to $(-d\cos \theta ,-d\sin | \theta | )$, and thus some functions and integration ranges change signs. In this case, the intersection points of the circle with the x-axis are the same as the ones given in Equation (A1), and thus all of the discussion given at the beginning of the previous subsection also applies for case (b). In particular, the intersection points (a₁, b₁) and (a₂, b₂) derived for case (a) are the same for this case.

In this case, the area of interest is the sum of areas A₁ and A₂. The former is the integral of the semicircle circle from x = −R_S to a₂, which is an integral that was already found in Equation (A9). As for area A₂, this is the integral of the upper part of the circle of radius R, i.e., of the function

$$\begin{eqnarray*}&&y=\sqrt{{R}^{2}-{\left(x+d\cos \theta \right)}^{2}}-d\sin | \theta | .\end{eqnarray*}$$

However, the integral of this from x = a₂ to a₁ is exactly the same integral calculated in case (a1), whose result is in Equation (A4), because the integrand there was the same integrand that we have here but multiplied by −1, and the limits of integration there were reversed with respect to the ones we have here (i.e., they went from a₁ to a₂); because inverting the limits of integration is the same as calculating the integral multiplied by –1, both effects cancel out. Thus, area A₂ in our case is area A₁ in case (a1). Thus, for case (b), we have A_S = A₁+A₂, i.e.,

$$\begin{eqnarray*}&&\boxed{{A}_{S}=\displaystyle \frac{{R}_{S}^{2}}{2}\left[h({a}_{2})+\displaystyle \frac{\pi }{2}\right]+\displaystyle \frac{{R}^{2}}{2}{\rm{\Delta }}f+d\sin | \theta | ({a}_{2}-{a}_{1})}.\end{eqnarray*}$$

A.1.3. Case (c)

Case (c) focuses on when $| d\sin \theta | \lt R$; i.e. there are two solutions for the intersection points of the circle with the line y = 0.

More specifically, case (c1) also occurs when ∣x_int∣ ≤ R_S , i.e. when the two solutions for the intersection points are inside the semicircle and the part of the circle above the intersection points is completely enclosed within the semicircle (see Figure 9). This can be quantitatively described by theoretically "extending" the semicircle into a full circle of radius R_S . The coordinates of the intersection of these two circles (setting the center of the circle of radius R_S at the origin) can be found by substituting x²+y² = R_S ² into ${(x+d\cos \theta )}^{2}+{(y+d\sin \theta )}^{2}\,=\,{R}^{2}$ to give the y-coordinates,

$$\begin{eqnarray*}&&y=-A\sin \theta \pm \cos \theta \sqrt{{R}_{S}^{2}-{A}^{2}},\end{eqnarray*}$$

where

$$\begin{eqnarray*}&&A=\displaystyle \frac{{R}^{2}-{R}_{S}^{2}-{d}^{2}}{2d}.\end{eqnarray*}$$

This is similar to Equations (A2) and (A3), except the positions of the circles relative to the origin have been changed slightly.

Therefore, case (c1) applies when ∣x_int∣ ≤ R_S and one of the following occurs.

• 1.  
  There is no solution for the intersection of the two circles. This will occur when ${A}^{2}\gt {R}_{S}^{2}$.
• 2.  
  Both of the y-coordinate solutions are real and negative, i.e. when $-A\sin \theta \pm \cos \theta \sqrt{{R}_{S}^{2}-{A}^{2}}\lt 0$.

Case (c2) occurs when ∣x_int∣ > R_S and there are two intersection points on the curved edge of the semicircle. For this case, the y-coordinates of the intersection between the circle R and the full circle of radius R_S must be positive. Using the same equations as above, this is when $-A\sin \theta \pm \cos \theta \sqrt{{R}_{S}^{2}-{A}^{2}}\geqslant 0$, where A is defined as in case (c1).

Case (c3) occurs when ∣x_int∣ ≤ R_S and there are two further intersection points on the curved edge of the semicircle, making a total of four intersection points. Therefore, case (c3) is when ∣x_int∣ ≤ R_S and $-A\sin \theta \pm \cos \theta \sqrt{{R}_{S}^{2}-{A}^{2}}\geqslant 0$.

To solve case (c1), the points of intersection in the base of the semicircle can be obtained via Equation (A1):

$$\begin{eqnarray*}&&{x}_{\mathrm{int}}^{\pm }=-d\cos \theta \pm \sqrt{{R}^{2}-{d}^{2}{\sin }^{2}\theta }.\end{eqnarray*}$$

The problem is then just calculating the area of the segment A_S , which is a well-known geometric problem with the solution

$$\begin{eqnarray*}&&\boxed{{A}_{S}={R}^{2}\ \arccos \left(\displaystyle \frac{y}{R}\right)-{xy},}\end{eqnarray*}$$

where

$$\begin{eqnarray*}\begin{array}{rcl}x & = & (1/2)({x}_{{\rm{int}}}^{+}+| {x}_{{\rm{int}}}^{-}| ),\\ & = & \sqrt{{R}^{2}-{d}^{2}{\sin }^{2}\theta },\\ y & = & -d\sin \theta .\end{array}\end{eqnarray*}$$

To solve case (c2), the problem can be set up by first theoretically "extending" the semicircle into a full circle of radius R_S and calculating the area of the intersection of the two circles using the equations described in Kreidberg (2015):

$$\begin{eqnarray}&&{A}_{{\rm{int}}}={R}^{2}\ \arccos \ u+{R}_{S}^{2}\ \arccos \ v-(1/2)\sqrt{w},\end{eqnarray} \tag{ A11 }$$

where

$$\begin{eqnarray*}\begin{array}{rcl}u & = & ({d}^{2}+{R}^{2}-{R}_{S}^{2})/(2{dR}),\\ v & = & ({d}^{2}+{R}_{S}^{2}-{R}^{2})/(2{{dR}}_{S}),\\ w & = & (-d+R+{R}_{S})(d+R-{R}_{S})(d-R+{R}_{S})(d+R+{R}_{S}).\end{array}\end{eqnarray*}$$

Using A_int, one can find A_S by subtracting half the area of the "extended" circle of radius R_S to yield

$$\begin{eqnarray*}&&\boxed{{A}_{S}={A}_{\mathrm{int}}-\pi {R}_{S}^{2}/2}.\end{eqnarray*}$$

To solve case (c3), the problem needs to be split up into two parts. The first part involves finding the area of intersection, A_int, of the circle and the semicircle "extended" into a full circle of radius R using the same method from part (c2), and the second involves finding area A₁. As can be seen from Figure 9, once these two areas are found, it is simply a matter of subtracting A₁ from A_int to find A_S .

To find A₁ is a very similar problem to case (c1); the points of intersection along the base of the semicircle can be found using Equation (A1), and then the problem reduces to that of finding the area of a segment, which is A₁ in this case. Following a similar method for case (c1),

$$\begin{eqnarray*}&&{A}_{1}={R}^{2}\arccos \left(\displaystyle \frac{-y}{R}\right)+{xy},\end{eqnarray*}$$

where x and y are defined the same as for case (c1); however, note the change in sign of y, due to the change of orientation of the shapes. As mentioned, A_int is the same as in case (c2) and described in Equation (A11).

Therefore, the total intersection area A_S is given by

$$\begin{eqnarray*}&&\boxed{{A}_{S}={A}_{\mathrm{int}}-{A}_{1}}.\end{eqnarray*}$$

In this paper, we have defined θ as the angle between the base of the semicircle and the line that extends from the base of the semicircle to the center of the circle (of length d). Importantly, θ is defined as positive when extending clockwise, assuming the center of the circle is to the left of the semicircle, which, as explained in Section A.1, due to the symmetry of the problem, can always be achieved by flipping the frame of reference.

Here θ is calculated from parameters that correspond to how the system physically appears in the sky. To explain, we place the center of the star at the origin of a 2D Cartesian coordinate system. Due to the symmetry of the star, the planet is assumed to move from left to right (horizontally) across it so that the dY/dX gradient of the direction of motion of the planet at t₀ (inferior conjunction) is zero (the planet moves along the X-direction only). This coordinate system is different from the frame of reference used in Section A.1, where the frame is rotated so that the center of the base of the semicircle is at the origin and the base lies along the y = 0 axis.

Therefore, θ depends on the angle the semicircle is rotated through with respect to the X-direction in the 2D system defined above, ϕ, the impact parameter, b, and whether the semicircle is originally to the left or right of the center of the circle. The latter is characterized by the time, t, and the time of inferior conjunction, t₀, and by assuming the planet moves from left to right. It is important to note that ϕ will be static if the planet moves in an almost straight line across the star (along the X-direction only). However, due to orbital mechanics, the orbital path may appear curved, and as such, the planet will rotate slightly as it crosses the star. This rotation will perturb ϕ in such a way as discussed in Section A.3.

In this paper, ϕ is defined in the range −π/2 to π/2, from the X-axis to the base of the top semicircle, where a positive angle is defined going in the counterclockwise direction. Once this is used to find the θ for R_p,1, to obtain the θ for R_p,2, one simply has to multiply θ by –1, due to the change in direction.

Table A.2 shows how θ is calculated for the different values of ϕ, b, and t for R_p,1.

t	b	ϕ	θ
≤t0	Positive	$\lt \arccos b/d$	$-\phi -\arcsin b/d$
		$\geqslant \arccos b/d$	$-\pi +\phi +\arcsin b/d$
≤t0	Negative	$\leqslant -\arccos | b| /d$	$\pi +\phi -\arcsin | b| /d$
		$\gt -\arccos | b| /d$	$-\phi +\arcsin | b| /d$
>t0	Positive	$\lt -\arccos b/d$	$-\pi -\phi +\arcsin b/d$
		$\geqslant -\arccos b/d$	$\phi -\arcsin b/d$
>t0	Negative	$\leqslant \arccos | b| /d$	$\phi +\arcsin | b| /d$
		$\gt \arccos | b| /d$	$\pi -\phi -\arcsin | b| /d$

Download table as:  ASCIITypeset image

Using the same coordinate system as described in A.2, if the planet were to move in a straight line across the star, then the planet's impact parameter, b, would stay constant at the value of the impact parameter defined at t₀ (the time of inferior conjunction). From Seager (2010), this is defined as

$$\begin{eqnarray}&&b=\displaystyle \frac{a\cos i}{{R}_{{\rm{star}}}}\left(\displaystyle \frac{1-{e}^{2}}{1-e\sin \varpi }\right),\end{eqnarray} \tag{ A12 }$$

where a is the semimajor axis, i is the inclination of the orbit, e is the eccentricity, ϖ is the longitude of periastron, and R_star is the radius of the star. However, due to the orbital motion of the planet around the star, from the perspective of our X–Y system, the orbital path of the planet is curved across the sky, so the impact parameter becomes a function of time. For the most accurate model, b used in the equations in Section A.2 should be adjusted to

$$\begin{eqnarray}&&b=-r\sin (\omega +f)\cos i,\end{eqnarray} \tag{ A13 }$$

where

$$\begin{eqnarray}&&r=\displaystyle \frac{a(1-{e}^{2})}{1+e\cos f},\end{eqnarray} \tag{ A14 }$$

$$\begin{eqnarray}&&\cos f=\displaystyle \frac{\cos E-e}{1-e\cos E},\end{eqnarray} \tag{ A15 }$$

$$\begin{eqnarray}&&E-e\sin E=\displaystyle \frac{2\pi }{T}(t-{t}_{0}).\end{eqnarray} \tag{ A16 }$$

Here T is the orbital period of the planet, t is the current time of interest, and t₀ is the time of periastron passage. The last equation can be solved numerically using the Newton–Raphson method.

Furthermore, due to the change in gradient of the planet's orbit with respect to this coordinate system, the angle between the base of the semicircle and the X-axis will change slightly as it passes across the star. Because of how the method of deriving θ from ϕ is configured, we need this changing angle at every time step instead of using the original "static" ϕ. The change in angle caused by this movement will be labeled ψ and should be added to the original ϕ to create an "updated" ϕ_new. To derive ψ and make it clearer, b has been relabeled y, where x and y are the coordinates of this 2D system with the center of the star at (0, 0). By differentiating,

$$\begin{eqnarray}&&{d}y=-r\cos (\varpi +f)\cos i\ {d}f.\end{eqnarray} \tag{ A17 }$$

Also explained in Seager (2010), the x-coordinate of the center of the semicircle will move as

$$\begin{eqnarray}&&x=-r\cos (\varpi +f),\end{eqnarray} \tag{ A18 }$$

so

$$\begin{eqnarray}&&{d}x=-r\sin (\varpi +f)\ {d}f.\end{eqnarray} \tag{ A19 }$$

Therefore, ψ can be calculated from

$$\begin{eqnarray}&&\tan \psi =\displaystyle \frac{{d}y}{{d}x},\end{eqnarray} \tag{ A20 }$$

$$\begin{eqnarray}&&\psi =\arctan [\cot (\varpi +f)\cos i].\end{eqnarray} \tag{ A21 }$$

Then, ϕ should be adjusted so that

$$\begin{eqnarray}&&{\phi }_{{\rm{new}}}={\phi }_{{\rm{old}}}+\psi .\end{eqnarray} \tag{ A22 }$$