One of the main endeavors of the field of exoplanetary sciences is the characterization of exoplanet atmospheres on a population level. The current method of choice to accomplish this task is transmission spectroscopy where the apparent radius of a transiting exoplanet is measured at multiple wavelengths in search of atomic and molecular absorption features produced by the upper atmosphere constituents. To extract the planetary radius from a transit light curve it is necessary to account for the decrease in luminosity away from the center of the projected stellar disk known as limb darkening. Physically motivated parametrizations of limb darkening in particular of the quadratic form are commonly used in exoplanet transit lightcurve fitting. Here we show that such parametrizations can introduce significant wavelengthdependent biases in the transmission spectra currently obtained with all instrument modes of the JWST and thus have the potential to affect atmospheric inferences. To avoid such biases we recommend the use of standard limbdarkening parametrizations with wide uninformative priors that allow for nonphysical stellar intensity profiles in the transit fits and thus for a complete and symmetrical exploration of the parameter space. We further find that fitting the light curves at the native resolution results in errors on the measured transit depths that are significantly smaller compared to light curves that are binned in wavelength before fitting thus potentially maximizing the amount of information that can be extracted from the data.
['Exoplanet atmospheres', 'Exoplanet detection methods', 'Exoplanet atmospheric composition', 'Limb darkening', 'James Webb Space Telescope', 'Transits', 'Exoplanet astronomy', 'Exoplanets', '487', '489', '2021', '922', '2291', '1711', '486', '498', 'Astrophysics - Earth and Planetary Astrophysics', 'Astrophysics - Instrumentation and Methods for Astrophysics', 'Astrophysics - Solar and Stellar Astrophysics']
Biases in Exoplanet Transmission Spectra Introduced by Limbdarkening Parametrization
{'Biases in Exoplanet Transmission Spectra Introduced by Limb Darkening Parametrization': "<!-- image -->  \n1 Department of Physics and Trottier Institute for Research on Exoplanets, Universit'e de Montr'eal, Montr'eal, QC, Canada louis-philippe.coulombe@umontreal.ca", 'ABSTRACT': 'One of the main endeavors of the field of exoplanetary sciences is the characterization of exoplanet atmospheres on a population level. The current method of choice to accomplish this task is transmission spectroscopy, where the apparent radius of a transiting exoplanet is measured at multiple wavelengths in search of atomic and molecular absorption features produced by the upper atmosphere constituents. To extract the planetary radius from a transit light curve, it is necessary to account for the decrease in luminosity away from the center of the projected stellar disk, known as the limb darkening. Physically-motivated parametrizations of the limb darkening, in particular of the quadratic form, are commonly used in exoplanet transit light-curve fitting. Here, we show that such parametrizations can introduce significant wavelength-dependent biases in the transmission spectra currently obtained with all instrument modes of the JWST, and thus have the potential to affect atmospheric inferences. To avoid such biases, we recommend the use of standard limb-darkening parametrizations with wide uninformative priors that allow for non-physical stellar intensity profiles in the transit fits, and thus for a complete and symmetrical exploration of the parameter space. We further find that fitting the light curves at the native resolution results in errors on the measured transit depths that are significantly smaller compared to light curves that are binned in wavelength before fitting, thus potentially maximizing the amount of information that can be extracted from the data.', '1. Introduction': "The transit method is arguably the most powerful tool currently at our disposition to detect exoplanets (Charbonneau et al. 2000) and characterize their atmospheres (Charbonneau et al. 2002). With the wide wavelength coverage and unprecedented photometric precision of the JWST, our understanding of exoplanet atmospheres through transmission spectroscopy is now rapidly and constantly evolving (e.g., JWST Transiting Exoplanet Community Early Release Science Team et al. 2023, Gao et al. 2023, Bell et al. 2023, Madhusudhan et al. 2023, Dyrek et al. 2023, Xue et al. 2024, Benneke et al. 2024, Welbanks et al. 2024, Cadieux et al. 2024, and many more).  \nAs a planet transits its host star from an observer's point of view, the measured stellar flux decreases by a factor of D ≈ ( R p /R ⋆ ) 2 , enabling the extraction of the planetary radius R p from this signal, provided that the stellar radius R ⋆ is known. When performed with spectroscopic instruments, this method yields the wavelength-dependent apparent radius of an exoplanet and is highly sensitive to the presence of atomic or molecular absorbers in its atmosphere (Seager & Sasselov 2000).  \nWhile the depth of a transit is, to first order, dictated by the planetary radius, the exact shape of the light curve also depends on the limb-darkening profile of the  \nhost star as well as the orbital parameters of the system (Seager & Mallen-Ornelas 2003). The limb darkening describes the non-uniform luminosity from the stellar disk, which appears brighter along an observer's line-of-sight in the center of the disk than near the limbs, whereas the orbital parameters define the path of the transiting planet (transit chord) across said stellar disk. The effect of limb darkening is highly wavelengthdependent due to the different optical depths and effective temperatures probed when observing the center (head-on) or the limb (slanted geometry) of the disk, which result in large center-to-limb contrast at optical wavelengths that gradually decreases towards the infrared. Because of this, limb-darkening coefficients are generally fitted independently at each wavelength when performing transmission spectroscopy.  \nSeveral parametrizations of the limb darkening have been proposed thus far, such as: linear (Milne 1921), quadratic (Kopal 1950), non-linear (Claret 2000), square-root (Diaz-Cordoves & Gimenez 1992), and exponential (Claret & Hauschildt 2003), with the quadratic limb-darkening law currently being the most widely used in exoplanetary sciences. This parametrization describes the normalized stellar intensity profile I ⋆ ( µ ) /I 0 as a quadratic function of the normalized radial coordinate µ :  \nI ⋆ ( µ ) I 0 = 1 -u 1 (1 -µ ) -u 2 (1 -µ ) 2 , (1)  \nwhere u 1 and u 2 are the quadratic limb-darkening coefficients (LDCs). Due to the flexibility of this parametrization, it is possible for a given set of parameters [ u 1 , u 2 ] to produce a stellar intensity profile where the intensity increases towards the limbs (limb brightening) or reaches negative values, scenarios that are non-physical for stars. To avoid considering such solutions when fitting a transit light curve, constraints and reparametrizations of the LDCs have been proposed (e.g. Burke et al. 2007, Carter et al. 2009, Kipping 2013), with the physically-motivated parametrization of Kipping (2013) (from hereon simply referred to as [ q 1 , q 2 ]) now commonly used for exoplanet transit light curves.  \nExplicitly, the [ q 1 , q 2 ] parametrization is obtained by imposing two physical conditions to the quadratic limbdarkening problem: the stellar intensity profile must be everywhere-positive, and it must decrease monotonically away from the center of the stellar disk (Kipping 2013). This parametrization is especially useful as it significantly reduces the LDCs parameter space and allows for the uniform sampling of q 1 and q 2 values over a simple unit square.  \nEven though the physical assumptions used to define these reparametrizations of the LDCs are valid, we find that a subtle but significant bias can arise from them: light-curve fits are biased towards smaller planetary radii when the limb-darkening effect becomes negligible. This is most relevant for transmission spectroscopy at the infrared wavelengths probed by the JWST (Gardner et al. 2006), the decommissioned Spitzer (Werner et al. 2004), and future missions such as the Ariel (Pascale et al. 2018) and Twinkle (Stotesbury et al. 2022) space telescopes. While we specifically focus on the bias introduced by the [ q 1 , q 2 ] parametrization in this work, we note that any prior or parametrization enforcing physical solutions for the stellar intensity profile are prone to such biases.  \nIn this work, we present, describe, and quantify the bias that arises from physics-driven reparametrizations of the limb-darkening problem when used in exoplanet transit light-curve fitting. In section 2, we describe the methods used to quantify the bias on simulated and real observations. We present the results in section 3, and we discuss the implications for atmospheric inference, as well as present a bias-free reparametrization in section 4. Finally, our conclusions are presented in section 5.", '2. Methods': 'In this section, we describe the methodology used to quantify the bias introduced by physically-motivated limb-darkening parametrizations on the measured transit depth. We begin by characterizing this effect on simulated observations and proceed to quantify this effect on the JWST NIRISS/SOSS observations of WASP-39 b afterward.', '2.1. Light-curve fitting of simulated data': 'We begin this work by performing transit light-curve fitting on simulated observations where the ground truth is known. We simulate transit light curves using the batman python package (Kreidberg 2015) for the case of a planet in a circular orbit ( e = 0, ω = 90 · ) with planet-to-star radius ratio R p /R ⋆ =0.1, semi-major axis a/R ⋆ = 4, impact parameter b = 0, and orbital period P = 1 day. We simulate the transit assuming no limb-darkening ( u 1 = 0, u 2 = 0), which is analogous to transit observations at infrared wavelengths where the bias on the planetary radius is strongest. The simulated transit light curve consists of a thousand equallyspaced integrations spanning times t ∈ [ -0 . 25 , 0 . 25] days, where t = 0 is the mid-transit time. We consider error envelopes for the simulated observations rather than adding white noise to the model to avoid the introduction of randomness in the retrieved parameter constraints. For the fits, we consider seven photometric scatter values ranging from 1 to 1000 ppm with equal spacing in log scale. We note that the values of photometric scatter considered here are for an integration time of 0.72 minutes and thus these values cannot be directly compared to observations that have a different cadence.  \nWe subsequently fit each light curve twice: fitting for R p /R ⋆ ( U [0,0.5]), u 1 ( U [-3,3]) and u 2 ( U [-3,3]) for the standard LDCs case as well as R p /R ⋆ , q 1 ( U [0,1]) and q 2 ( U [0,1]) for the physically-motivated LDCs case. We consider wide priors for u 1 and u 2 to ensure that a diverse range of stellar intensity profiles may be modelled, including limb-brightening scenarios. The light curves are fitted considering the standard chi-square likelihood and the parameter space is sampled using the emcee python package (Foreman-Mackey et al. 2013). We use 4 walkers per free parameter and run the chains for 100,000 steps, tossing the first 60% of the steps as burnin.', '2.1.1. Light-curve fitting of the NIRISS/SOSS observations of WASP-39 b': 'We proceed to explore the impact of the choice of limb-darkening parametrization on the NIRISS/SOSS spectrum of WASP-39b. A transit of WASP-39b was observed on 26 July 2022 using the NIRISS instrument  \nNormalized Flux  \n0  \nI  \n/  \n)  \n(  \ns  \nI  \ne  \nl  \ni  \nf  \no  \nr  \nP  \ny  \nt  \ni  \ns  \nn  \ne  \nt  \nn  \nI  \n0.9910  \n0.9905  \n0.9900  \n0.9895  \n1.1  \n1.0  \n0.9  \n0.00  \n0.25  \n0.50  \n0.75  \n1.00  \n0.04  \n0.02 0.00  \n0.02  \n0.04  \nTime [days]  \nFigure 1. Impact of the choice of limb-darkening on the light curve fit and measured transit depth for the case of no limb darkening and photometric scatter of σ = 1000ppm. Top left: Median stellar intensity profile I s ( µ ) /I 0 (blue line) from the transit fit to the simulated observations using the [ u 1 , u 2 ] limb-darkening parametrization. The 1σ and 2σ confidence intervals are shown by the shaded regions. Bottom left: Median transit model (blue) obtained from the transit fit with the [ u 1 , u 2 ] limbdarkening parametrization. The 1σ and 2σ confidence intervals are shown by the shaded regions. Top middle: Same as Top left considering the [ q 1 , q 2 ] parametrization. Bottom middle: Same as Bottom left considering the [ q 1 , q 2 ] parametrization. Right: Probability density distributions of the measured transit depth minus the true transit depth for the [ u 1 , u 2 ] (blue) and [ q 1 , q 2 ] (red) parametrizations. The transit depth measured using [ u 1 , u 2 ] is centered on the expected value, whereas the one measured using [ q 1 , q 2 ] peaks at a smaller value.  \n<!-- image -->  \n0.00  \n0.25  \n0.50  \n0.75  \n1.00  \n0.04  \n0.02 0.00  \n0.02  \n0.04  \nTime [days]  \n(Doyon et al. 2023) of the JWST in the SOSS mode (Albert et al. 2023) as part of the Early Release Science (ERS) program (Bean et al. 2018), the full analysis of which is presented in Feinstein et al. (2023). WASP-39 b is currently the only exoplanet that has been observed in transmission with all instrument modes of the JWST (Ahrer et al. 2023, Alderson et al. 2023, Feinstein et al. 2023, Rustamkulov et al. 2023, Powell et al. 2024, Carter et al. 2024), making it an ideal test case for the exploration of the bias introduced by limb-darkening assumptions. We further choose the NIRISS/SOSS instrument since it covers a large wavelength range (0.85-2.85 µ m) at a high resolving power (R = 500 - 1400), and therefore should display a transition from light curves with important limb darkening to light curves where the limb darkening becomes negligible and the bias on the planetary radius arises.  \nWe begin from the raw spectroscopic light curves extracted from the order 1 spectral trace ( λ = 0.852.85 µ m) using the NAMELESS reduction (Coulombe et al. 2023, Radica et al. 2023, Lim et al. 2023, FournierTondreau et al. 2023) presented in Feinstein et al. (2023). We consider a linear trend for the systematics model with a normalization factor c ( U [-10 9 ,10 9 ]) and slope v ( U [-10 9 ,10 9 ]). The planetary radius R p /R ⋆ and LDCs are kept free for both limb-darkening parametriza-  \ns, considering the same priors as in section 2.1. We fix the orbital parameters to T 0 = 2459787.556740 BJD, a/R ⋆ = 11.4, and b = 0.4510 following the values measured from the NAMELESS order 1 white-light curve in Feinstein et al. (2023). For both the [ u 1 , u 2 ] and [ q 1 , q 2 ] parametrizations, we fit the light curves at native resolution (i.e., one light curve per detector column) and at a fixed resolving power of R = 80. Each light curve fit is run for 10,000 steps using 4 walkers per free parameter, where we again discard the first 60% of the steps as burn-in.', '3.1. Simulated data': "We find that, because of the physical assumptions from which the [ q 1 , q 2 ] parametrization is derived, the transit fit is unable to fully and symmetrically explore the parameter space around solutions with low limbdarkening (Figure 1). To produce a symmetric ensemble of models around a transit light curve in the limit of no limb-darkening, the model must allow for unphysical scenarios with limb brightening, where the stellar intensity profile instead increases from the center to the limb. Because the [ q 1 , q 2 ] parametrization is unable to produce such models, the fit instead compensates for this by decreasing the planetary radius (Figure 1).  \nRadial Coordinate  \nRadial Coordinate  \nTruth  \n[u1,u2]  \n[q1,q2]  \nProbability Density  \nFigure 2. Magnitude of the bias on the measured transit depth as a function of photometric scatter. The data points show the deviation from the true transit depth, along with the 1σ uncertainty, as a function of photometric scatter for both the [ u 1 , u 2 ] and [ q 1 , q 2 ] parametrizations. We use a symmetric log scale for the y-axis. Contrastingly from the [ u 1 , u 2 ] parametrization, the [ q 1 , q 2 ] parametrization shows an increasing bias of the transit depth when going toward larger scatter. The blue dashed line corresponds to no bias on the transit depth, with the measured transit depths obtained using the [ u 1 , u 2 ] parametrization following that trend. The bias on the transit depth when using the [ q 1 , q 2 ] parametrization closely follows the red dashed line, which traces the trend y = -√ 2 /π σ t,d , where σ t,d is the measured transit depth uncertainty.  \n<!-- image -->  \nThis effect is reminiscent of other astronomical measurements where parameter exploration near a hard boundary or the exclusion of non-physical scenarios results in biased parameter inference. For instance, Lucy & Sweeney (1971) found spectroscopic binaries in circular orbits to be biased towards small non-zero eccentricities due to the rigid boundary at e = 0. This problem is also present in radial velocity fits where, in the limit of low signal-to-noise, negative planetary masses must be allowed in the parameter exploration in order not to bias the measurements (e.g., Fulton et al. 2018, Eastman et al. 2019).  \nThe magnitude of the bias on the measured transit depth introduced through the [ q 1 , q 2 ] parametrization scales with the photometric scatter (Figure 2), and reaches approximately 85 ppm for a photometric scatter of 1000 ppm. While we keep the current investigation of the bias to its simplest by only fitting for R p / R ⋆ and the LDCs, we note that increasing the complexity of the model (e.g., fitting for additional parameters such as b and a / R ⋆ ) can potentially increase the bias as it provides more flexibility to the transit model. We also find that the magnitude of the bias relative to the measured transit depth uncertainty ( σ t , d ) does not vary with the light-curve photometric scatter, and is close to  \n√ 2 /π σ t , d ≈ 0 . 8 σ t , d . This value corresponds to the mean of a half-normal distribution, consistent with the scenario of a measured parameter whose true value is situated directly at the boundary of the prior space.  \nIn the context of transmission spectroscopy, this means that a spectrum that is fitted at a high resolution and then binned to a lower resolution will show systematically lower transit depths than a spectrum obtained by fitting light curves that are already binned to the same low resolution. This is because the highresolution spectrum is more sensitive to the bias since the light curves that are fit have larger scatter. This effect will especially arise at infrared wavelengths where the limb darkening becomes less important. Furthermore, the magnitude of the bias on the measured transit depths for a given spectrum will follow the wavelengthdependent scatter of the light curves, which depends on the instrument throughput and stellar spectrum. The consequence of this effect on JWST transmission spectra has already been observed by May et al. (2023) and Carter et al. (2024), who then recommended to bin the light curves at a lower spectral resolution before proceeding with the light curve fitting (effectively mitigating the effect of the bias by decreasing the scatter in the light curves) and/or to fix the LDCs to avoid this bias.  \nWhile fixing the LDCs to theoretical model predictions leads to spectra that are consistent no matter the resolution at which they are fit, such models can differ significantly from reality and could in turn bias the measured transmission spectrum (e.g. Espinoza & Jord'an 2015, Maxted 2022, Patel & Espinoza 2022, Rustamkulov et al. 2023, Kostogryz et al. 2024, Sarkar et al. 2024). The continuous temporal coverage and high photometric precision of the JWST allows for a precise measurement of the LDCs from the data and the use of theoretical models should preferably be reserved to planets with high impact parameters, in which case the LDCs and planetary radius are strongly correlated (e.g. Roy et al. 2023, Radica et al. 2024). We propose instead to use the [ u 1 , u 2 ] parametrization with wide uninformative priors such that the resulting transmission spectrum is free of any bias no matter the spectral resolution used for the fit, which we demonstrate for real observations in the following subsection.", '3.2. NIRISS/SOSS observations of WASP-39 b': 'We find that the bias arising from the [ q 1 , q 2 ] limbdarkening parametrization significantly changes the NIRISS/SOSS transmission spectrum of WASP-39b redwards of two microns. The spectra obtained by fitting the light curves at native resolution and at a fixed resolution of R = 80 for both the [ u 1 , u 2 ] and [ q 1 , q 2 ]  \nFigure 3. Impact of the choice of limb-darkening parametrization on the measured transmission spectrum. Top left: NIRISS/SOSS spectra of WASP-39b, fit at a spectral resolution of R = 80 (blue points) and at native resolution (black points), obtained using the [ u 1 , u 2 ] parametrization with wide uniform priors. The fit performed at native resolution is binned at a resolution of R = 80 after the light curve fitting. Bottom left: Difference between the fit performed at a resolution of R = 80 and the one at native resolution (blue points) for the [ u 1 , u 2 ] parametrization. The error bars for the residuals are computed in quadrature. A smoothed version of the residuals (blue line), along with the 1σ confidence region, are shown using a linear Savitzky-Golay filter with a window size of 31 points. Top right: Same as Top left for the [ q 1 , q 2 ] parametrization. Bottom right: Same as Bottom left for the [ q 1 , q 2 ] parametrization. The residuals between the transit spectra obtained from the light-curve fits at different resolutions clearly show a systematic deviation away from 0 past 2 µ m in the [ q 1 , q 2 ] case.  \n<!-- image -->  \nparametrizations are shown in Figure 3. When the light curves are fitted using the [ q 1 , q 2 ] parametrization, the native resolution spectrum is systematically lower past 2 µ m than the fit performed on the binned light curves, reaching a discrepancy of up to 200 ppm at the longest wavelengths. This discrepancy between the spectra fit at different resolutions disappears when using the [ u 1 , u 2 ] parametrization.  \nIn the case of NIRISS/SOSS, it is mainly the longest wavelengths that are affected by this bias since this is where the scatter is highest and also where the limbdarkening effect becomes negligible. We note that this effect is not significant only for NIRISS/SOSS but for all instruments onboard the JWST due to their coverage of infrared wavelengths. Furthermore, because the instruments onboard the JWST all have different spectral resolutions and throughputs, using the [ q 1 , q 2 ] parametrization may also lead or contribute to offsets between transmission spectra of the same planet observed with multiple instruments. However, such offsets have also been observed in datasets where the LDCs have been fixed  \n(e.g. Wallack et al. 2024), meaning other causes might also be at play.  \nWe find that the errors on the measured WASP-39 b NIRISS/SOSS transit depths are 20 to 60% lower across the spectral order 1 wavelength range for the spectrum where the light curves have been fitted at the native resolution compared to the one that was directly fitted at a resolution of R = 80 (Figure 4). We measure this difference in precision in each spectroscopic channel via  \nRel . difference = 100 · [ σ R=80 /σ native -1] % ,  \nwhere σ R=80 and σ native are the 1 σ uncertainties on the spectroscopic transit depths obtained when the light curves are binned before and after the fit, respectively. We find that the gain in precision is most notable towards the blue and red ends of the wavelength range (Figure 4), where the flux is lower and thus where there is the most increase in precision to be made. One possible explanation for this increase in precision is that, as the light curves are binned, the ratio of the systematicto-Poisson noise σ sys /σ Poisson increases, forcing the fitted scatter to compensate for this deviation from pure  \nFigure 4. Fitting the light curves at higher spectral resolution produces a more precise spectrum. Relative difference in percent between the errors on the transit depths measured at a resolution of R = 80 and those measured at the native resolution and then fitted afterward. A smoothed version of the difference transit depth precision (black line) is shown using a linear Savitszky-Golay filter with a window size of 31 points. The errors obtained from the lower resolution fit are systematically larger than those obtained at the native resolution, with the difference in precision ranging from 20 to 60% across the wavelength range.  \n<!-- image -->  \nPoisson noise. However, when we bin the spectrum that was fitted at the native resolution, we assume that the uncertainties decrease following Poisson noise, which could explain the improvement in precision compared to a spectrum that is fitted at a lower resolution. Another possibility is that some of the fit parameters, such as the systematics model parameters and limb-darkening coefficients, introduce uncertainties on the measured transit depth that do not scale following Poisson noise behavior, thus proving advantageous to fit the light curves at a higher spectral resolution. Hence, we conclude that fitting the light curves at the native resolution might be necessary to achieve the best possible precision for transmission spectroscopy.  \nWhile we find that a uniform prior U [ -3 , 3] allows for sufficient exploration of the LDCs parameter space for the simulated data as well as the WASP-39 b spectrum presented in this work, we note that this range might need to be adjusted depending on the precision of the observations. One should also ensure that the retrieved LDCs are physical and that there is no statistically significant preference for limb-brightening/unphysical scenarios in the light curve fits.', '4. Discussion': '4.1. Effect on atmospheric inferences  \nWhile the [ q 1 , q 2 ] parametrization bias described above results in measured transit depths that are mostly consistent to one standard deviation with those obtained using the [ u 1 , u 2 ] parametrization (Figures 1, 2, 3), its effect on atmospheric inference can be significant. This is because this bias is repeated over many wavelength bins and its magnitude may depend on wavelength for observations where the scatter vary significantly over the bandpass. This bias can thus lead to unphysical trends in the obtained transmission spectra (e.g., slopes, offsets between instruments or detectors) which cannot be reproduced by atmosphere models and must be accounted for with ad-hoc methods such as fitting for detector offsets at the retrieval stage or error inflation. Furthermore, the amplitude of the bias can reach the same order of magnitude as the atmospheric features we are aiming to probe in transmission spectroscopy. Looking at the NIRISS/SOSS spectrum of WASP-39b specifically, the bias introduced by the [ q 1 , q 2 ] parametrization reaches 200 ppm at the longest wavelengths whereas the atmospheric features themselves have an amplitude of approximately 1000 ppm. Such a systematic downward trend at long wavelengths, affecting practically half of the wavelength range covered by the observations, could impact the inferred atmospheric metallicity, which depends on the relative strength of the water features (Benneke & Seager 2012), as well as the need for non-gray clouds and/or hazes to explain the spectrum. This is even more significant for smaller terrestrial planets, at which point the expected amplitude of the transmission spectroscopy features for a secondary atmosphere are on the level of ∼ 100 ppm and the transit depth precision can reach 25 ppm (e.g., Benneke et al. submitted) Performing light-curve fitting of JWST observations with the [ u 1 , u 2 ] parametrization is thus the best way to ensure that the observed transit spectra are free of these limb-darkening-induced biases before proceeding with atmospheric retrievals.', '4.2. Updating the [ q 1 , q 2 ] parametrization': 'Because the bias towards smaller planetary radii is introduced by excluding limb brightening from the parameter space and considering only limb-darkening scenarios, one must then ask: is it possible to relax the constraints in the [ q 1 , q 2 ] parametrization such that limb brightening is also allowed? We derive such a parametrization in the Appendix and find that one can sample uniformly from a [ r 1 , r 2 ] unit square to produce the following values of the quadratic LDCs:  \nu 1 = 2( r 2 √ r 1 + √ r 1 -1) (2)  \nu 2 = u 1 | u 1 | √ r 1 (1 -r 2 ) -u 1 2 , (3)  \nsuch that the values of u 1 and u 2 allowed result only in positive stellar intensity profiles that either monotonically decrease or increase towards the limbs. We find however that, other than for the advantage of reducing the size of the prior space and circumventing the bias on the planetary radius, this parametrization results in larger uncertainties on the measured transit depth and might also be trickier to explore by sampling algorithms (see details in the Appendix). It therefore seems that the best option is to simply use the standard [ u 1 , u 2 ] parametrization.', '5. Summary and Conclusions': 'We have demonstrated the existence of a bias on the measured planetary radius that is introduced through the use of physically-motivated limb-darkening parametrizations. This bias is most obvious at infrared wavelengths where the limb-darkening effect is low and the assumption of a monotonically decreasing stellar intensity profile hinders a complete exploration of the planetary radius parameter space. Tests on simulated data show that the measured planetary radius is systematically lower when using the [ q 1 , q 2 ] parametrization of Kipping (2013). We also show that this bias increases with photometric scatter explaining the discrepancy between spectra fit at high versus low spectral resolutions observed in May et al. (2023) and Carter et al. (2024). This bias can affect large portions of the wavelengths covered by a transmission spectrum, and lead to significant biases in atmospheric characterization currently obtained with the instruments onboard the JWST.  \nWe find that using the standard [ u 1 , u 1 ] quadratic limb-darkening parametrization with wide uninformative priors prevents this bias and demonstrate this on the JWST NIRISS/SOSS spectrum of WASP-39b. Furthermore, we find that fitting the light curves at the native resolution results in errors on the transit depth that are 20 to 60% lower across the NIRISS/SOSS wavelength range compared to the spectrum directly fitted at a low resolution, showing that it is potentially necessary to fit at the native resolution of the instruments to maximize the precision of the measured spectrum.  \nIn the limit of high-precision optical light curves, where the effect of limb darkening is important, and for which such physically-motivated parametrizations were developed (e.g. Burke et al. 2007, Carter et al. 2009, Kipping 2013), the use of the [ q 1 , q 2 ] parametrization (or similar) could still be advisable, as the limb-darkeningfree solutions should confidently be excluded from the  \nparameter space. Nevertheless, it is unclear exactly at which wavelengths and what level of photometric scatter the bias on the measured transit depth becomes important. As we move to the analysis of infrared light curves with the JWST and other future space telescopes, the use of such methods should be avoided to ensure that transmission spectra measurements are free of biases.', 'Acknowledgments': "We thank the anonymous reviewer for their insightful comments that improved the quality of this work. We would also like to acknowledge Aarynn Carter, Erin May, N'estor Espinoza, and Luis Welbanks for their role as part of the Early Release Science data synthesis effort in characterizing the discrepancy between JWST exoplanet transmission spectra fitted at low vs high spectral resolution. This work is based in part on observations made with the NASA/ESA/CSA/JWST. The data were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127. The specific observations analyzed can be accessed via DOI: 10.17909/efv0-j030. B.B. and P.-A.R. acknowledge funding from the Natural Sciences and Engineering Research Council (NSERC) of Canada. L.-P.C. and P.-A.R. acknowledge support from the University of Montr'eal, and from the Trottier Research institute for Exoplanets (iREx). B.B also acknowledges financial support from the Canadian Space Agency and the Fond de Recherche Qu'eb'ecois-Nature et Technologie (FRQNT; Qu'ebec).  \nSoftware: Astropy (Astropy Collaboration et al. 2013, 2018, 2022), batman (Kreidberg 2015), emcee (ForemanMackey et al. 2013), Matplotlib (Hunter 2007), NumPy (Harris et al. 2020), and SciPy (Virtanen et al. 2020).", 'APPENDIX': "We derive a new parametrization for the limb-darkening, closely following the methodology presented in Kipping (2013), where we relax the constraint of monotonically decreasing stellar intensity profiles that was used to derive the [ q 1 , q 2 ]. To do this, we impose the two following conditions:  \n- A. The normalized stellar intensity profile cannot reach values below 0 and, per symmetry in the case of limbbrightening, must also not go above 2.\n- B. The stellar intensity profile must either monotonically decrease or monotonically increase away from the center of the stellar disk.  \nCondition A can be expressed mathematically as the following relation:  \n-1 ≤ -u 1 (1 -µ ) -u 2 (1 -µ ) 2 ≤ 1 . (1)  \nWe can derive the boundary conditions imposed by this relation by evaluating it at its maxima ( µ = 0 and µ = 1). The µ = 1 limit returns the trivial condition -1 ≤ 0 ≤ 1. At the µ = 0 limit however, we have that the standard LDCs must follow the condition -1 ≤ -u 1 -u 2 ≤ 1, from which we can express the boundaries of this parameter space as u 2 = -u 1 -1 and u 2 = -u 1 +1. This is consistent with the condition derived in Kipping (2013) for the case of a positive-everywhere stellar intensity profile.  \nAs for condition B , it can be expressed mathematically as enforcing that the ratio of the gradient of the stellar intensity profile ( ∂I/∂µ = u 1 +2 u 2 (1 -µ )) between any two positions on the stellar disk must be equal to or larger than zero  \nu 1 +2 u 2 (1 -µ 1 ) u 1 +2 u 2 (1 -µ 2 ) ≥ 0 . (2)  \nBy evaluating this equation at its extrema ( µ 1 = 0, µ 2 = 1 and µ 1 = 1, µ 2 = 0), we find that the above condition results in the two following boundaries: u 1 ≥ 0, u 2 ≥ -u 1 / 2 and u 1 < 0, u 2 < -u 1 / 2. The regions of the parameter space that follow conditions A and B are shown in Fig. 5. We recover the Kipping (2013) condition (upper triangle in left panel of Figure 6) in our new bounded parameter space, with the addition of a second region that is now also allowed. This second region corresponds to cases where the stellar intensity profile increases away from the center of the disk (limb-brightening). We then change our system of coordinates such that the x-axis is u 1 / 2 and the y-axis is u 1 / 2+ u 2 , transforming our two regions into right-angled triangles. In order to transform this parameter space in a unit square from which we can sample uniformly the limb-darkening parameters, we must join these two regions together into a single triangle. We thus rotate the region that corresponds to limb-brightening scenarios about the u 1 / 2 axis using the parametrization w 1 ≡ u 1 / 2, w 2 ≡ sgn( u 1 )( u 1 / 2 + u 2 ). From this, we follow the same method as in Kipping (2013) and use the triangular sampling of Turk (1990) such that ( w 1 w 2 ) T = (1 -√ r 1 ) ⃗ A + √ r 1 (1 -r 2 ) ⃗ B + r 2 √ r 1 ⃗ C , where ⃗ A = ( -1 0) T , ⃗ B = (0 1) T , and ⃗ C = (1 0) T are the three vector positions corresponding to the vertices of the triangle in [ w 1 , w 2 ] space. We can then sample uniformly from the [ r 1 , r 2 ] unit square, the values of which can be converted to the quadratic LDCs via the following relations:  \nu 1 = 2( r 2 √ r 1 + √ r 1 -1) ; u 2 = u 1 | u 1 | √ r 1 (1 -r 2 ) -u 1 2 . (3)  \nHaving defined our updated limb-darkening reparametrization, we run the same comparison shown in Figure 1 but this time using the [ r 1 , r 2 ] parametrization, which is shown in Figure 6. While the [ r 1 , r 2 ] parametrization does not seem to induce a bias on the measured transit depth, we find that it actually results in a larger uncertainty on the measurement compared to simply fitting for [ u 1 , u 2 ]. We also derive a second parametrization by rotating the limb-brightening scenarios about the u 1 / 2 + u 2 axis and find that it produces virtually the same results as the [ r 1 , r 2 ] parametrization. This increase in the transit depth uncertainty when using [ r 1 , r 2 ] is most likely due to the fact that this parameter space is more difficult to explore. The transition from the limb-darkening to the limb-brightening  \nscenarios happens at the [ u 1 = 0 , u 2 = 0] 'bottleneck' (Fig. 5), meaning the amount of prior space that corresponds to solutions with low limb-darkening/brightening is significantly reduced when using [ r 1 , r 2 ] compared to the [ u 1 , u 2 ] parameter space.  \nFigure 5. Left: Region of the [ u 1 , u 2 ] parameter space where conditions A and B are respected. Center: . Reparametrization of the limb-darkening parameter space such that the two triangular regions forming the parameter space following conditions A and B are right-angled triangles. Right: Parameter space joining the limb-darkening and limb-brightening scenarios, obtained by rotating the limb-brightening scenarios about the u 1 / 2 axis.  \n<!-- image -->  \nFigure 6. Impact of the choice of limb-darkening on the light curve fit and measured transit depth for the case of no limb darkening and photometric scatter of σ = 1000ppm. Top left: Median stellar intensity profile I s ( µ ) /I 0 (blue line) from the transit fit to the simulated observations using the [ u 1 , u 2 ] limb-darkening parametrization. The 1 σ and 2 σ confidence intervals are shown by the shaded regions. Bottom left: Median transit model (blue) obtained from the transit fit with the [ u 1 , u 2 ] limbdarkening parametrization. The 1 σ and 2 σ confidence intervals are shown by the shaded regions. Top middle: Same as Top left considering the [ r 1 , r 2 ] parametrization. Bottom middle: Same as Bottom left considering the [ r 1 , r 2 ] parametrization. Right: Probability density distributions of the measured transit depth minus the true transit depth for the [ u 1 , u 2 ] (blue) and [ r 1 , r 2 ] (red) parametrizations. The transit depths measured using [ u 1 , u 2 ] and [ r 1 , r 2 ] are centered on the expected value, but the posterior distribution obtained with [ r 1 , r 2 ] is wider.  \n<!-- image -->", 'REFERENCES': 'Madhusudhan, N., Sarkar, S., Constantinou, S., et al. 2023, Carbon-bearing Molecules in a Possible Hycean  \nAtmosphere. https://arxiv.org/abs/2309.05566  \nMaxted, P. F. L. 2022, Monthly Notices of the Royal Astronomical Society, 519, 3723-3735,  \ndoi: 10.1093/mnras/stac3741  \n- May, E. M., MacDonald, R. J., Bennett, K. A., et al. 2023, The Astrophysical Journal Letters, 959, L9, doi: 10.3847/2041-8213/ad054f\n- Milne, E. A. 1921, Monthly Notices of the Royal Astronomical Society, 81, 361,  \ndoi: 10.1093/mnras/81.5.361  \nPascale, E., Bezawada, N., Barstow, J., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10698, Space Telescopes and Instrumentation 2018: Optical, Infrared, and Millimeter Wave, ed. M. Lystrup, H. A. MacEwen, G. G. Fazio, N. Batalha, N. Siegler, & E. C. Tong, 106980H, doi: 10.1117/12.2311838  \n- Patel, J. A., & Espinoza, N. 2022, AJ, 163, 228, doi: 10.3847/1538-3881/ac5f55\n- Powell, D., Feinstein, A. D., Lee, E. K. H., et al. 2024, Nature, 626, 979, doi: 10.1038/s41586-024-07040-9  \nRadica, M., Welbanks, L., Espinoza, N., et al. 2023, Monthly Notices of the Royal Astronomical Society, 524, 835-856, doi: 10.1093/mnras/stad1762  \nRadica, M., Coulombe, L.-P., Taylor, J., et al. 2024, Muted Features in the JWST NIRISS Transmission Spectrum of Hot-Neptune LTT 9779 b.  \nhttps://arxiv.org/abs/2401.15548  \nRoy, P.-A., Benneke, B., Piaulet, C., et al. 2023, The Astrophysical Journal Letters, 954, L52, doi: 10.3847/2041-8213/acebf0'}