TransitFit: combined multi-instrument exoplanet transit fitting for JWST, HST and ground-based transmission spectroscopy studies

We present TransitFit, a package designed to fit exoplanetary transit light-curves. TransitFit offers multi-epoch, multi-wavelength fitting of multi-telescope transit data. TransitFit allows per-telescope detrending to be performed simultaneously with transit parameter fitting, including custom detrending. Host limb darkening can be fitted using prior conditioning from stellar atmosphere models. We demonstrate TransitFit in a number of contexts. We model multi-telescope broadband optical data from the ground-based SPEARNET survey of the low-density hot-Neptune WASP-127b and compare results to a previously published higher spectral resolution GTC/OSIRIS transmission spectrum. Using TransitFit, we fit 26 transit epochs by TESS to recover improved ephemeris of the hot-Jupiter WASP-91b and a transit depth determined to a precision of 111 ppm. We use TransitFit to conduct an investigation into the contested presence of TTV signatures in WASP-126b using 180 transits observed by TESS, concluding that there is no statistically significant evidence for such signatures from observations spanning 27 TESS sectors. We fit HST observations of WASP-43b, demonstrating how TransitFit can use custom detrending algorithms to remove complex baseline systematics. Lastly, we present a transmission spectrum of the atmosphere of WASP-96b constructed from simultaneous fitting of JWST NIRISS Early Release Observations and archive HST WFC3 transit data. The transmission spectrum shows generally good correspondence between spectral features present in both datasets, despite very different detrending requirements.


INTRODUCTION
Over the last decade, the study of exoplanetary atmospheres though transmission spectroscopy studies has been a growing and maturing field, seeing a significant increase in the number of surveys targeting transiting exoplanets.Dedicated space-based transit surveys, such as the Kepler Space Telescope (Borucki et al. 2010), and the more recent Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2014) together with ground-based surveys such as the Next Generation Transit Survey (NGTS, Wheatley et al. 2013), have provided an ever-growing list of targets for study.These surveys have contributed significantly to the growing number of confirmed exoplanets: at the time of writing, there are over 5500 confirmed exoplanets, with over 3900 of these exhibiting observable transits 1 .Early surveys, such as the Wide Angle Search for Planets (WASP, Pollacco et al. 2006), were ground-based and therefore most of the exoplanets discovered in the early days of the field were bright enough to allow for ground-based follow-up studies.
Many of the planets discovered by recent surveys such as TESS can be observed with both ground-and space-based telescopes, and as such, the field of transmission spectroscopy is leaving the "targetstarved" era and entering an "asset-starved" era.The limiting factor to exoplanetary studies is now no longer the availability of viable targets, but instead the availability of ground-based facilities to conduct follow up.To adapt to this change, new tools and techniques need to be developed to best utilise available resources.
The identification of these species relies on two stages of retrieval.First, a spectrum must be acquired from measurements of the radius of a planet at different observation wavelengths, and then an atmospheric model can be used to obtain atmospheric parameters for the planet.Obtaining accurate planet-host radius ratios from light curves is a significant challenge as many factors affect the shape of a transit light curve, including atmosphere and orbital parameters of the exoplanet, the behaviour of the host, and instrumental and terrestrial factors.
There are currently a few publicly available codes designed for fitting light curves of transiting exoplanets.PyLightcurve (Tsiaras et al. 2016a) is a complete forward model and retrieval package, which uses an MCMC routine to fit transit light curves.Py-Lightcurve can also simultaneously remove trends from data using a 2nd-order polynomial and offers a variety of limb-darkening laws.exoplanet (Foreman-Mackey et al. 2020) is a toolkit for modelling transit and radial velocity observations of exoplanets, and is built with multi-planet systems in mind.Similarly to PyLightcurve, it uses MCMC methods to fit transit curves, and has a variety of limbdarkening laws, though it does not offer detrending functionality.ExoFastv2 (Eastman et al. 2019) is an IDL package that can perform simultaneous fitting of transit, radial velocity and astrometric observations.It is highly optimized and offers a host of user options, including simultaneous detrending, using an MCMC optimizer.
Approaches to fitting limb-darkening coefficients (LDCs) vary, with some researchers fixing values before retrieval and some instead fitting them as free parameters as part of retrieval.Espinoza & Jordán (2015) discuss this and conclude that it is best to freely-fit LDCs as fixing them can lead to biases of up to 3 per cent in measurements of   / ★ , which can have significant effects on retrieved spectra.Chen 1 The Extrasolar Planet Encyclopedia: http://exoplanet.eu/et al. ( 2018) demonstrated that by using information about the host star, namely temperature, mass (or surface gravity), and metallicity, it is possible to improve the fitting of LDCs and consequently the fitting of transit light curves in general.Whilst fitting LDCs as free parameters is common practise, it is clear that LDCs depend on observation wavelength and on the properties of the host star and therefore it is worth trying to develop tools that can exploit this additional information.
In this paper, we present TransitFit, an open-source Python 3 package for robust multi-wavelength, multi-epoch fitting of transiting exoplanet light curves obtained from one or more telescopes.TransitFit has been developed in response to the fast-growing numbers of available transmission spectroscopy targets, as part of the Spectroscopy and Photometry of Exoplanetary Atmospheres Research Network (SPEARNET), a survey that is employing automated transmission spectroscopy target selection for follow-up by a globally dispersed and heterogeneous telescope network (Morgan et al. 2019).In Section 2, we discuss the implementation of TransitFit, including the approach to limb-darkening, simultaneous detrending, and its handling of large multi-epoch, multi-wavelength datasets.In Section 3 we demonstrate the application of TransitFit to five different situations.We illustrate the use of TransitFit on multi-telescope, multi wavelength, multi-epoch observations of WASP-127 b obtained by SPEARNET.We then look at two examples of the use of multi-epoch TESS data, producing improved ephemeris for WASP-91 b and conducting a sensitive investigation into the presence of contested TTV signatures from WASP-126 b.Finally, we demonstrate how Tran-sitFit can be used to fit complex systematics through analysis of HST observations of WASP-43 b, and the capability of TransitFit to fit JWST NIRISS and HST WFC3 observations of WASP-96 b.We offer our conclusions in Section 4.

IMPLEMENTATION OF TransitFit
TransitFit is an open-source, pure Python 3.x package designed specifically with transmission spectroscopy studies in mind and uses transit observations at different wavelengths and epochs from different telescopes simultaneously to fit transit parameters using nested sampling retrieval.This approach allows lightcurves to be fitted by exploiting coupled information across wavelength and epoch.Tran-sitFit offers a wavelength-coupled approach to LDC fitting, which we discuss in depth in Section 2.1.We also discuss how TransitFit can detrend and normalise light curves simultaneously with fitting other parameters in Section 2.2.Different light curves being fitted simultaneously may benefit from using different detrending functions, especially if obtained from different telescopes.TransitFit can also deal with multi-epoch observations from systems which exhibit transit timing variations (TTVs), although, as highlighted in Section 2.3, there are some current limitations to fitting these systems.
The forward model used to calculate transit light curves in Tran-sitFit is the batman Python package (Kreidberg 2015).We have chosen batman as it is well-established within the community, and the open-source nature of it means that we can easily incorporate it into an open-source retrieval code.Since the choice of forward model is what limits the parameters which can be fitted, using batman means that the physical parameters retrievable by TransitFit are: orbital period, ; the time of inferior conjunction,  0 ; the planetstar radius ratio,   / ★ ; the ratio of semi-major axis of the planet orbit to the radius of the star, / ★ ; the orbital inclination, ; the orbital eccentricity, ; the longitude of periastron for the planet's orbit, ; and a set of limb-darkening coefficients, written as a matrix U, where  ,  is the th LDC for filter .Along with these, TransitFit can calculate detrending parameters and normalisation constants for individual light curves.
Retrieval is conducted using nested sampling routines (Skilling 2004;Skilling 2006) implemented with the dynesty Python package (Speagle 2020).During each iteration of fitting, we use batman Python package to calculate model light curve from the sampled parameters, while simultaneously detrending and normalising the raw lightcurves using the parameters sampled in the same iteration.The likelihood function for the fitting includes an additional errorrescaling term (   ) as described in the online documentation for EMCEE 2 (Foreman-Mackey et al. 2013).This allows for cases where the flux errors might be underestimated.The retrieval returns the highest-likelihood parameter values as the best fit parameters.We define upper limit on the best fit value as the 68.27 quantile of the weighted samples beyond the best fit value, and the lower limit as 31.73quantile of weighted samples below the best fit.
General use of TransitFit involves 3 configuration files, one directing TransitFit to the data to be fitted, one defining the priors, and one defining the observation filter profiles, and a single wrapper function, run_retrieval, which handles under-the-hood code interfacing.In order to be able to denote individual observations and their relationships, each light curve passed to TransitFit is identified by three indices: telescope, filter, and epoch.Observations which come from the same telescope will share a telescope index, and so on.By taking this three-index approach, TransitFit is able to easily fit shared parameters, such as   values for observations made in the same filter.In addition to the telescope, filter, and epoch indices, each observation is given an index referring to a specific detrending model (see Section 2.2).The models that each detrending index refers to is set by the user, and allows the use of more than one detrending model across a set of observations.

Limb-darkening
TransitFit uses the Limb Darkening Toolkit (LDTk) (Parviainen & Aigrain 2015) to calculate the likelihood of sets of limb-darkening coefficients given planet host characteristics and the wavelength of observations.This allows limb darkening priors to be constructed that tension (though not fix) the wavelength-dependent behaviour of limb darkening in accordance with the Phoenix stellar atmosphere models (Husser et al. 2013).As well as this "coupled" approach, TransitFit can compute limb-darkening coefficients independently for each filter, corresponding to the more traditional "uncoupled" approach that allows for easy comparison with other models.As a third option, TransitFit can use a "fixed" mode where LDC are fit in only one waveband and LDTk is used to compute LDC values for other observed wavebands.This latter mode is offered for use in situations where coupling is desired, but a large number of filters leads to an unreasonable number of parameters to be fitted.However, such a fixed LDC approach can lead to biased results (Espinoza & Jordán 2015).We therefore advocate using TransitFit in coupled mode whenever feasible.

Limb-darkening models
Stellar intensity varies between the centre and the edge of the disk, which often results in the base of a transit curve being rounder than if the stellar disk were a uniform brightness.Typically, these variations 2 https://emcee.readthedocs.io/en/stable/tutorials/line/ in intensity are described by analytical functions   (), where  is the cosine of the angle between the line of sight and the emergent intensity. can also be expressed as  = √ 1 −  2 where  is the unit-normalised radial coordinate on the stellar disk, and as such, all limb-darkening models must be valid for 0 ≤  < 1.
There are a wide variety of limb-darkening models (Claret 2000), and several have been implemented in TransitFit.These are the linear law (Schwarzschild & Villiger 1906) the quadratic law (Kopal 1950), the square-root law (Diaz-Cordoves & Gimenez 1992), the power-2 law (Morello et al. 2017) and the non-linear law (Claret 2000) (5) Each of  0 ,  1 ,  2 , and  3 are the LDCs, which must be fitted simultaneously with other transit parameters during retrieval.The LDCs are, however, dependent upon wavelength, and consequently must be fitted for each filter used in an observation or set of observations.

Constraining limb-darkening coefficients
The most basic approach for fitting limb-darkening coefficients is to sample them independently and find the best fit values, but this can allow unphysical values to be trialled with no penalty.This is an issue which was addressed by Kipping (2013), who imposed two conditions on limb-darkening profiles to ensure that they are physically allowed: (i) The intensity profile must be always positive, or The intensity profile must be monotonically decreasing from the centre to the edge of the stellar disk, meaning that Kipping (2013) showed that by applying these conditions, it is possible to place constraints on allowed values of the LDCs for the two-parameter quadratic and square root laws.To improve the efficiency of sampling within this restricted space, rather than sampling the LDCs { 0 ,  1 }, Kipping (2013) instead reparameterises the laws in terms of the coefficients { ( 0 ,  1 ) ,  ( 0 ,  1 )} ∈ [0, 1].This reparameterisation ensures a sampling efficiency of 100 per cent (i.e.all samples are within the physically allowed region), without the need for checking that sampled values of  0 and  1 follow the imposed constraints.An alternate way of parameterisation for the power-2 law coefficients has been discussed in Short et al. (2019).
We have implemented the parameterisation following Kipping (2013) in TransitFit, and have extended the method to the power-2 law in Equation 4. In the limit  → 0, Condition (i) yields the constraint  0,p2 < 1. (8) Condition (ii) implies that This does not give us anything overtly useful in the limit  → 0 due to the cross terms, however, as  → 1, we see that This places the constraint that the power-2 LDCs must have the same sign, and can be viewed as a 'quadrant limiting' constraint.Mathematically, there is no lower bound on  0,p2 , and there are no bounds at all on  1,p2 but, from a computational perspective, we must place limits on them in order to be able to sample values.Therefore, we can say that  0,p2 ∈   0,p2 , 1 and  1,p2 ∈   1,p2 ,   1,p2 .We implement this by fitting  2 ,  2 ∈ [0, 1] and using the conversions and Provided that   1,p2 = −  1,p2 , it can be shown that this method uniformly samples in all the allowed regions in the  0,p2 ,  1,p2 plane, with 100 per cent efficiency.
It is trivial to also apply this method to the linear law, which places the constraint 0 <  0,l < 1, (13) and this has been implemented in TransitFit.However, this method has yet to be successfully applied to the non-linear law.The current best attempt is by Kipping (2016), where the methodology is extended to the three-term law of Sing et al. (2009), which drops the  1/2 term from the non-linear law in Equation 5. Consequently, TransitFit does not use the Kipping parameterisation to limit sampling of non-linear LDCs to a physically-allowed region of parameter space.

Coupling limb-darkening coefficients across wavelengths
Multiple codes use the Kipping parameterisation to constrain the LDC values to those which are physically allowed, but it is possible to improve the quality of LDC fitting further.All of the currentlyavailable transit-fitting codes fit LDCs for each filter independently.This means that for each filter a transit is observed at, the best-fit LDCs for each filter may not be physically consistent with each other.Parviainen & Aigrain (2015) developed the Limb Darkening Toolkit (LDTk) to allow researchers to address this problem, but we have been unable to find a publicly available transit fitting code that makes use of LDTk.LDTk uses the library of PHOENIX stellar atmospheres and synthetic spectra (Husser et al. 2013) to estimate the likelihood of a set of stellar LDCs for a given set of observation filters.Using this, we have given TransitFit the functionality to couple LDCs across multiple filters, which can then be fitted simultaneously.This allows for the refinement of the limb-darkening physics included in transmission spectroscopy studies by ensuring that the retrieved LDC values are statistically tensioned across filters in a manner consistent with stellar atmosphere models.The filter profiles used by TransitFit can be either uniform, box filter profiles, which may be suitable to represent individual spectroscopic channels, or user-supplied filter profiles for broadband photometric studies.In the case where a specific filter profile cannot be obtained, we recommend using the equivalent width of a filter.TransitFit is distributed with filter profiles for the Johnson-Cousins UVBRI set, and the SLOAN-SDSS u'g'r'i'z' set, as well as profiles for Kepler and TESS.All of these were obtained from the SVO Filter Profile Service3 (Rodrigo et al. 2012;Rodrigo & Solano 2020).
Since there may be edge cases where this coupling is not desired, for instance where host information is unavailable, or computational limitations in the case of fitting a very large numbers of observation filters, such as spectroscopic observations, TransitFit offers three modes of LDC fitting: Independent: This is the traditional approach of fitting LDCs for each filter separately.TransitFit still uses the Kipping parameterisations laid out in Section 2.1.2,but LDTk is not used to couple LDCs across filters.
Coupled: Using the Kipping parameterisations, each LDC is fitted as a free parameter, with LDTk being used to estimate the likelihood of sets of LDCs, using information on the host star and the observation filters.TransitFit also provides the functionality to use the uncertainty multiplier from LDTk.
Single filter: When fitting with multiple wavebands, the number of parameters required to be fitted can increase dramatically when using the coupled mode.Consequently, we have provided a method of only freely fitting the LDCs of one filter, and using LDTk to extrapolate LDC values for the remaining filters.For the -th coefficient of a filter  ,  ,  , this extrapolation is calculated by where   is the sampled value of the -th LDC in the actively fitted filter, and ⟨ ,  ⟩ and ⟨  ⟩ are the maximum likelihood values initially suggested by LDTk.

Detrending and normalisation of light curves
Transit light curves are sensitive to a variety of factors which stop the out-of-transit baseline being flat.These can range from host variations through to internal reflections within the telescope.These trends must be removed in order to obtain accurate parameters from observations.Additionally, many transit models, including batman, normalise the out-of-transit baseline to a flux of 1.In some cases, detrending and normalisation is conducted before any further analysis, but ideally, detrending and normalisation coefficients should be fitted simultaneously with other model parameters to ensure that light curve features of interest are not inadvertently removed.
In order to ensure that detrending does not bias our measurement of the transit depth we impose a constraint on the detrending function such that it conserves the relative flux at the epoch of mid-transit,  0 , (i.e we conserve the transit depth).We refer to detrending functions which meet this criterion as "depth-conserving."The enforcement of depth conservation provides an important constraint for normalising the detrending function.Similar relative photometry conservation is also at the heart of other differential photometry methods, such as difference imaging (Alard & Lupton 1998).Note that since  0 is itself a parameter that is being fitted for, the detrending correction and the  0 determination are correlated and so are determined by TransitFit simultaneously.
TransitFit offers functionality to simultaneously detrend and normalise light curves during retrieval.Built into the package are depth-conserving, -th order detrending functions, and the user can supply their own custom functions if more complicated detrending is required.
The -th order functions are calculated by writing the detrended flux values of a time series t, D (t) as where F (t) are the raw flux values and d (t) is some detrending function.We place the constraint of depth conservation upon the detrended light curves such that which gives us the constraint that By applying this conservation of transit depth at t 0 , we place constraints on the 0-th component (intercept) of the detrending function.
In the case of a linear detrending function, where applying Equation 17yields from which we can cast the 0-th order term  in terms of the other parameters to give which can be substituted into Equation 18, resulting in This can be generalised to -th order (for  > 0) detrending functions given by where   are the detrending coefficients and the exponent of the time series is bit-wise.This method allows us to also fit a normalisation constant without falling foul of degeneracy between the scaling due to the normalisation constant and the shift that a freely-fitted 0-th order detrending term introduces.0-th order detrending cannot be applied due to this degeneracy, and we assume that these light curves are detrended, but not necessarily normalised.In the case of a user-defined detrending function, we strongly recommend following the above depth-conservation procedure in order to avoid the risk of degenerate solutions.

Dealing with systems exhibiting TTVs
The basic implementation of TransitFit assumes that there are no transit timing variations (TTVs) within multi-epoch observations and fits one value of  0 , assumed to be consistent across all epochs.In the event that a system does exhibit TTVs, this method will fail to produce an accurate result.Consequently, in these cases, TransitFit takes a slightly different approach.
(i) First, we consider each filter separately.We run retrieval on all the curves in this filter, using all the data to fit   and limb-darkening coefficients.However, we fit a separate  0 for each observation epoch within the filter, and cannot directly fit a period, , in this mode, which must instead be provided.
(ii) Using the results from these single-filter retrievals, we detrend and normalise each light curve and then use the retrieved  0 values to produce a phase-folded light curve for each filter.The observation times t for each light curve are folded to give t ′ , where where accounts for the offset caused by the different  0 values for each epoch.By choosing a value for  0,base this term ensures that all the light curves are centred on  0,base .
(iii) Fit the folded light curves using the standard TransitFit approach, coupling LDCs where required.
As stated above, when allowing for the presence of TTVs, Tran-sitFit cannot fit for , which must be provided.For consistency, we recommend first running TransitFit on data assuming that no TTVs are present, in order to obtain an appropriate value for .TransitFit cannot automatically detect TTVs in data, and must be instructed explicitly to allow for them.In the case where TTVs are present but TransitFit is not allowing for them, the retrieved results will be incorrect.Additionally, TransitFit does not solve the system dynamics associated with any present TTVs, as the purpose of TransitFit is to produce robust transit fitting for the purposes of transmission spectroscopy studies.

Batched retrieval: fitting a large number of parameters
As with any retrieval algorithm, increasing the dimensionality of the parameter space leads to instability in the nested sampling routines and can lead to inaccurate results.Since TransitFit is anticipated to be used in transmission spectroscopy studies, where many tens, or even hundreds of light curves may need to be fitted, we have provided a solution to this in the form of "batched" retrieval.
In this mode, the user can specify the maximum number of parameters for TransitFit to fit at one time.The light curves are then grouped by observation filter and split into multi-filter batches, where the number of parameters being fitted in each batch is less than the user-set limit.The batches are calculated to try and ensure that filters are present in multiple batches, which results in coupling between them.The exception to this is in the case that one filter has a high enough number of observations in it that the number of parameters required exceeds the user-set limit.In this case, this filter is fitted independently and does not benefit from any coupling.In these cases, we recommend using the "folded" mode.After retrieval has been run on all the batches, a set of summary results are produced by calculating a weighted mean of all parameters.

Folded retrieval: Producing folded light curves
With the launch of large surveys such as TESS, many exoplanets have multiple-epoch observations in a single filter.TransitFit can make use of these through a two-step retrieval process.In the first step, TransitFit runs a retrieval on each filter independently, and uses the results to produce a phase-folded light curve for each filter.In the second step, TransitFit runs a standard multi-wavelength retrieval using the batched algorithm above.This mode of retrieval allows for the production of high-quality folded light curves from nondetrended data, as well as providing a method where observations from long-term, single-waveband surveys such as TESS can be easily combined with single-epoch observations at multiple wavelengths, such as from ground-based spectrographic follow-up.

APPLICATIONS OF TransitFit
TransitFit was initially designed for use in spectroscopy studies, but also be applied to temporal studies, either in updating ephemeris of planets, or in studying systems for TTVs, which can be indicative of other planets in a system.
In this section, we demonstrate the application of TransitFit in four different scenarios, illustrating the impact of using LDTk to inform LDC fitting.First, we will discuss the fitting of multiwavelength, ground based photometric observations of the lowdensity hot Neptune WASP-127 b (Lam et al. 2017), using previously unpublished data acquired from the SPEARNET network of telescopes.We then move to applying TransitFit to TESS observations of the warm Jupiter WASP-91 b (Anderson et al. 2017) and provide updated ephemeris and orbital parameters for the system.Third, we analyse TESS observations of the hot Jupiter WASP-126 b (Maxted et al. 2016), a system which contentiously exhibits TTVs (Pearson 2019;Maciejewski 2020), and use TransitFit to show that there is no statistically significant evidence of TTVs within 180 transits observed by TESS.Fourth, we analyse a single spectroscopic channel of the HST observation of WASP-43 b made by Kreidberg et al. (2014) to demonstrate the capability of TransitFit to handle complex systematics through custom detrending functions.Finally, we fit the JWST and HST observations of WASP-96 b simultaneously to generate a transmission spectrum and show the capability of TransitFit to work with JWST observations.For all five systems, we assume circular orbits and use TransitFit to fit the global parameters of ,  0 , / ★ , and , as well as the filter-specific   / ★ and LDC values.We use an uncertainty multiplier of 10 for the LDC values, which is based on comparison of results from different grids in ExoCTK (Bourque et al. 2021) and ExoTiC-LD (Grant & Wakeford 2022).

Broadband photometric observations of WASP-127 b
SPEARNET is a prototype transmission spectroscopy survey which is utilising a globally-distributed network of heterogeneous optical telescopes, the locations of which are shown in Figure 1.It was conceived to anticipate and address the challenges that the transition into the so-called "asset-starved" era poses, primarily by designing tools which allow for increased utility of resources, both before (Morgan et al. 2019) and after (Hayes et al. 2020) transit observations.Tran-sitFit was conceived as part of the SPEARNET suite of tools to handle transit data from non-homogeneous observations, to facilitate transmission spectroscopy studies in the asset-staved era, as time on larger telescopes is becoming ever-more competitive and studies will have to frequently rely on data taken from a combination of telescopes.
As part of the network operation, Morgan et al. ( 2019) developed a metric for ranking candidates for observation, effectively pairing targets with telescopes in a way which maximises the signal-tonoise of the observations.The motivation behind this metric is to remove the multiple unquantifiable biases in manual transmission spectroscopy target selection.Since the selection function is known, it is possible to make population-corrected statistical statements based on observations that are guided by the metric.In Table 3 of Morgan et al. (2019), we show that WASP-127 b is consistently ranked in the top three targets for a variety of telescopes when known planet masses are included in the metric calculations, and as such it has become a target of interest for SPEARNET.
With a density of 0.07±0.01 Jup (Lam et al. 2017), WASP-127 b is one of the lowest-density planets so far discovered, and occupies the 'short-period Neptune desert' (Mazeh et al. 2016), which is notable since most planets with its characteristics are not expected to survive due to photo-evaporation (Haswell et al. 2012).Its low density also makes WASP-127 b an idea target for transmission spectroscopy due to its large atmospheric scale height, and several studies have been completed, with potential detections of water (Chen et al. 2018;Skaf et al. 2020), and statistically significant detections of sodium, potassium, and lithium (5, 3, and 4 respectively, Chen et al. 2018).No significant evidence for helium in the upper atmosphere has been found (dos Santos et al. 2020) and it has been proposed that this is due to unfavourable photo-ionisation conditions.
The approaches to LDC fitting in these previous studies differ.Skaf et al. (2020) fix the LDCs for all spectral channels at the white-light values predicted for WASP-127 using the quadratic law of Claret (2000).Chen et al. (2018) also use a quadratic limb-darkening law, but instead find the highest likelihood values for each channel and fit using a Gaussian prior of width 0.1, sourcing the initial predictions from the Kurucz ATLAS9 stellar atmosphere models (Kurucz 2017).
Using the SPEARNET telescope network, we have obtained six photometric light curves in four different wavebands from five transits of WASP-127 b, including the first published transits observed in the  ′ -band.The three telescopes used in these observations were: The 2.4m Thai National Telescope (TNT): Located at the Thai National Observatory (TNO), the TNT observations of WASP-127b were conducted using ULTRASPEC (Dhillon et al. 2014), which uses a 1024×1024 pixel high-speed frame-transfer EMCCD camera with a field-of-view of 7.68 × 7.68 arcmin 2 .The dead time between exposures on this setup is 14 ms.
A 0.7 m telescope at Gao Mei Gu observatory (TRT-GAO) The TRT-GAO is also part of the Thai Robotic Telescope Network, and is located at Gao Mei Gu observatory in Lĳiang, China.Observations were taken using an Andor iKon-L 936 2048 × 2048 CCD with a field of view of 20.9 × 20.9 arcmin 2 .
The 0.6 m PROMPT-8 telescope Located at the Cerro Tololo Inter-American Observatory (CTIO) in Chile, PROMPT-8 is a 0.6 m telescope operated by the Skynet Robotic Telescope Network.Imaging is conducted on this telescope using a 2048 × 2048 pixel CCD camera with a scale of 0.624 arcseconds/pixel.
Figure 1 shows the location of the primary telescopes in the SPEARNET network, with the telescopes used in this study of WASP-127 b highlighted.The precise details of the observations taken are given in Table 1.
The images obtained by SPEARNET were calibrated using IRAF routines along with astrometric calibration using Astrometry.net(Lang et al. 2010).In order to obtain the light curve, aperture photometry was carried out using sextractor (Bertin & Arnouts 1996) using an adaptive scaled aperture based on the seeing in an individual image.Reference stars were chosen to have magnitudes similar to WASP-127 (|Δ | < 3) and no intrinsic variation.The observed time in BJD TDB with flux ratio between WASP-127 b and reference stars with their error are shown in Table 2.
The obtained light curves, shown in 2(a), were then run through TransitFit, using the 'batched' mode, a 2nd-order detrending function, using quadratic limb-darkening model, and both 'coupled' and 'independent' LDC fitting approaches so as to able to identify the improvement from using filter and host parameters to inform the likelihood of LDC values.For the coupled LDC approach, we adopted the stellar parameters of Lam et al. (2017), namely  = 5620 ± 85 K,  ★ = 1.39±0.03 ⊙ ,  ★ = 1.08±0.03 ⊙ , log  = 4.18±0.01cgs, and [Fe/H] = −0.18± 0.06 dex.For the priors, we take the values from Chen et al. (2018).The value of  0 was scaled using  to match the time span of the raw lightcurves, to reduce extrapolation during fitting.
The resulting detrended light curves from the coupled LDC fitting are shown with the best-fit model in Figure 2 2020) have all been analysed separately, there is an intrinsic offset between all the data.We have compared the data with our retrieved spectrum in Figure 3.The wavelength positions of the SPEARNET observations in Figure 3 are derived by weighting the relevant filter profile by the spectral energy distribution (SED) of WASP-127, predicted using the PHOENIX models.
Looking at Figure 3, from coupled LDC fitting approach, the value of   / ★ is somewhat different than the independent LDC fitting method, with a discrepancy of around 8 per cent in the  ′ and  bands.This difference is larger than the 3 per cent bias that Espinoza & Jordán (2015) find can be introduced from not fitting LDCs at all, and clearly illustrates the impact that using host characteristics and filter profiles can have when fitting spectroscopic and photometric measurements.
This could also be due to the incompleteness of the observation, as in the case of -band and  ′ -band, which do not have a complete transit observed.When using the coupled LDC approach, the fact that this observation does not have a complete transit may have an effect on the parameters retrieved for all the filters.Table 3 shows the physical results obtained by each of the two retrievals described above.
Looking at Figure 3, it may be noticed that the transit depth precision obtained from the TNT observations are roughly comparable to those obtained by Skaf et al. (2020) from their HST observations.At first glance, this may appear surprising, as one would instinctively expect that space-based observations would produce higher-quality results than those made from ground-based observatories.However, this fails to account for the impact of the high on-sky efficiency of observations made using ULTRASPEC as well as the difference between broadband imaging and spectroscopy.Assuming that calculations are photon-noise limited, we expect the signal-to-noise ratio (S/N) on the retrieved flux value to go as √︁   , where   is the number of photons collected.Using this, S/N for the transit depth can be estimated by √︁  ,t 14 +  ,base where  ,t 14 and  ,base are the number of photons in-transit and in baseline respectively.In the source-dominated limit   can be estimated by where  int is the on-sky integration time,  is the host apparent magnitude and  0 the instrument zero-point magnitude.
For TNT in  ′ the zero-point magnitude is  0,TNT = 25.255 and with 42,680 in-transit observations with a mean exposure of 0.38 s this gives  ,t 14 ,TNT ≃ 1.7 × 10 10 integrated over all intransit exposures of WASP-127 ( ′ = 10.0).With 26,088 baseline observations, we also get  ,base,TNT ≃ 1.0 × 10 10 .Skaf et al. (2020) observed 38 in-transit observations using the HST Wide Field Camera 3 (WFC3) G141 grism with a mean exposure of 95.782 s.Using the WFC3IR spectroscopic exposure time calculator6 , and adopting a G5V host spectrum for WASP-127 normalised to  = 10.15 mag, we find a photon count integrated over all in-transit exposures of  ,t 14 ,HST ≃ 3.4 × 10 9 at 1.4 m for a spectral resolution of 70.The corresponding baseline count gives us  ,base,HST ≃ 3.2 × 10 9 .This gives a rough ratio of S/N for TNT to S/N for HST, of √︄  ,t 14 ,TNT +  ,base,TNT  ,t 14 ,HST +  ,base,HST ≃ 4.2.
Clearly, the approximate equality here masks the fact that we are comparing the throughput of a broadband filter on TNT to the sensitivity of a single spectral bin from the G141 grism.

Multi-epoch study of WASP-91 b with TESS
WASP-91 b is a 1.34  Jup warm Jupiter with a 2.8 day orbit around a K3 host star (Anderson et al. 2017).With a radius of 1.03  Jup , WASP-91 b is a smaller example of a hot Jupiter, and has not been the subject of further studies since its discovery.TESS observed WASP-91 b in Sectors 1, 27, and 28, capturing a total of 50 158 photometric data points covering 26 transits, the data which we acquired from the Barbara A. Mikulski Archive for Space Telescopes (MAST) portal7 .We use 120 s cadence data, and PDC_SAPFLUX values are taken as flux.Since these data contain vast numbers of observations outside of transit, and there are multiple offsets and other trends within the data which cannot be modelled with a simple polynomial, we estimate the transit duration  14 , assuming a circular orbit with a 90 degree inclination, using and discard any data which are more than 2.5 14 away from the midtransit times predicted by the ephemeris given in Anderson et al. (2017).This results in 26 transits, capturing a combined total of

Gaomeigu Observatory Lijiang
Regional Observatory for the Public Regional Observatory for the Public Regional Observatory for the Public Thai National Observatory Table 1.Details of the SPEARNET observations of WASP-127 b, obtained using the Thai National Telescope (TNT), telescope from the Thai Robotic Telescope network at Gao Mei Gu (TRT-GAO), and the PROMPT-8 telescope at the Cerro Tololo Inter-American Observatory.Since the observations were taken, TRT-TNO has been upgraded to a 1 m aperture telescope.9 789 photometric points, which can then be individually normalised and detrended with low-order polynomials.This splitting of data is provided in TransitFit through the split_lightcurve_file function.

Date
After splitting the full data into individual transit observations, we use the "folded" mode of TransitFit to fit these data, using both the "coupled" and "independent" LDC fitting modes, assuming a circular orbit and using the quadratic limb-darkening model from Equation 2 and a second-order detrending polynomial.The values from Anderson et al. (2017) were used as priors for our fitting.For the "coupled" mode, we inform LDC fitting using the stellar parameters found by Anderson et al. (2017), given in Table 4, and the TESS filter profile given on by the SVO Filter Profile Service8 (Rodrigo et al. 2012;Rodrigo & Solano 2020).
Figure 4 shows the fitting process for the "coupled" run of Tran-sitFit at various stages.Figure 4 a(i) and a(ii) show the raw TESS observations which clearly exhibit various offsets and long-term trends.In order to reduce the impact of these, we split the data into individual transits using the approach described above, and Figure 4(b) shows the resulting raw data for the first transit.In the "folded" fitting mode, each transit is normalised and detrended, and Figure 4(c) shows the first transit after this first stage of processing.We have overplotted the final best fit model to help illustrate this step, but it should be noted that this model is calculated from all the transits, not just this single epoch.Once normalisation and detrending has been fitted for each transit, all the light curves are folded together and this curves is analysed to retrieve the final best-fit model.Figure 4(d) shows this folded light curve, along with the final best-fit transit model and residuals, which have an rms of 0.0017.We also show the same folded data binned to a cadence of two minutes, to demonstrate the improvement in observation precision when compared to the single-transit TESS observations like the one shown in Figure 4(c).An average of nearly 25 photometric points are contained within each of the bins, which through naïve root- statistics suggests a maximum improvement in precision of factor 5. The rms of the residuals of this binned light curve is 0.0004, which is an improvement of factor 4.8.The similar- ity of these two factors suggest that the improvement due to folding is near maximal, and thus the binned data are not systematics-limited.
We present the results from both runs of TransitFit alongside the results from Anderson et al. (2017) in Table 4, which are all consistent with each other.Precise ephemerides are required for accurate study of TTVs, and updating ephemeris by applying TransitFit to planets within TESS data will prove invaluable to future surveys.
The uncertainties on the two TransitFit retrievals are generally comparable, with the notable exception of the LDCs.Through using LDTk to calculate LDC likelihoods, we see upto an order-ofmagnitude increase in the precision of LDC values, which demonstrates the impact of introducing host parameters and filter information into transit-fitting routines.

TTV analysis of WASP-126 b
Orbiting a type G2 star with a period of 3.2888 ± 0.0008 days, WASP-126 b (Maxted et al. 2016) is a 0.28 ± 0.04  Jup hot Jupiter which has been identified as potentially exhibiting TTVs.Through Bayesian -body simulation coupled with machine learning analysis of Sectors 1-3 of the TESS observations of WASP-126 b, Pearson (2019) showed that there was evidence of a TTV signal with amplitude of ∼ 1 minute and a period of ∼ 25 days, which could be attributed to a non-transiting planet with   = 0.202 ± 0.077  Jup on a 7.63 ± 0.17 day orbit, dubbed WASP-126 c. Maciejewski (2020) studied the TESS observations from sectors 1-13 and found that when the extra sectors were included, the TTV signal was not present at a statistically significant level.
Here we use the ability of TransitFit to account for TTVs to analyse the TESS observations of WASP-126 b from sectors 1-13, 27-34, 36-39, 61, and 63 to further investigate the presence of TTVs within these data, and to produce the most up-to-date values for the planetary and orbital parameters of WASP-126 b.As with the analysis of WASP-91 b, we use PDC_SAPFLUX values for flux, and we discard all data that is more than 2.5 14 from the mid-transit times predicted by Maxted et al. (2016), which gives 180 individual transits.We then run analysis using TransitFit in both "coupled" and "independent" LDC fitting modes, using the quadratic limb-darkening law, Table 3.Output parameters from using TransitFit to run retrieval on photometric observations of WASP-127b, in both coupled and independent LDC fitting mode.For comparison, we provide the orbital parameters obtained by Chen et al. (2018) and Skaf et al. (2020).We have converted the Chen et al. (2018) 4. Output parameters for retrieval using TransitFit in "folded" retrieval mode on TESS observations of 26 transits of WASP-91 b, using both "coupled" and "uncoupled" LDC fitting approaches, alongside those obtained by Anderson et al. (2017).The data were fitted using a 2nd-order detrending model and the quadratic limb-darkening law given in Equation 2, and assuming zero orbital ellipticity, using the profile of the TESS filter available on the SVO Filter Profile Service to calculate LDC likelihoods when using LDTk.The TransitFit retrieval uses host parameters fixed at those obtained by Anderson et al. (2017), which are used in tandem with the filter profile to inform the LDC likelihood fitting.a 2nd order detrending polynomial, assuming a circular orbit, and using results from Maxted et al. (2016) as priors.When using LDTk to calculate LDC likelihoods, we use the host parameters given in Maxted et al. (2016), which are presented in Table 5, and the TESS filter profile from the SVO Filter Profile Service.As discussed in Section 2.3, we first run analysis of the data assuming that there are no TTVs, and we use these results to provide priors to parameterise the ephemeris for the analysis where we allow for the presence of TTVs.Since TransitFit requires a fixed period to be provided when considering TTVs, this step should always be used to ensure complete consistency of results.The priors used in this initial step are based on the ephemeris of Maxted et al. (2016).The resulting ephemerides from this initial analysis are given in Table 5, and we use these results when allowing for the presence of TTVs, fixing the period at the values given and using the retrieved  0 values as the mean of a Gaussian prior with a width of 0.007 days.
The orbital parameters of WASP-126 b retrieved by TransitFit when allowing for TTVs are given in Table 5, for both "coupled" and "independent" LDC fitting modes.We present these alongside the results of Pearson (2019), Maciejewski (2020), andMaxted et al. (2016) for comparison.We find that the results from both runs are generally consistent but, as in the WASP-91 b analysis, the uncertainties on the LDCs for the "coupled" run are smaller.
We present the O-C plots for both modes of the TransitFit analysis in the left plots of Figure 5, with the associated Lomb-Scargle periodograms on the right.The top row are the results for the "coupled" LDC run, and the bottom row are the results for "independent" LDCs.The solid horizontal lines on the Lomb-Scargle periodograms represent the false alarm probabilities of 10, 5, and 1 per cent from bottom to top, calculated using astropy routines (Astropy Collaboration et al. 2013, 2018).For the "independent" LDCs, the O-C plot The raw data from the first observed transit, chosen purely as an illustration; c: The first observed transit after the first stage of "folded" mode retrieval.These data have been detrended and normalised, the error bars have been scaled using the retrieved   , and have the final best-fit model overplotted; d: the final phase-folded light curve, best-fit transit model and residuals, containing 9 789 photometric data points, and, in red, the same data binned to a 2 minute observation cadence.The error bars have been scaled using   , and the corresponding  2 value has also been displayed.has a reduced chi-squared value of   =0.87, whilst the "coupled" LDC O-C gives   =0.85.We find that there are no periodicities of statistical significance within the O-C data for either the "coupled" or "independent" LDC runs, and consequently conclude that there is no evidence of TTVs that would be indicative of a second planet in the WASP-126 system, in agreement with the findings of Maciejewski (2020).
The approach to fitting taken by TransitFit differs to both Pearson (2019) and Maciejewski (2020), and consequently this result can be taken to be an independent verification of the findings of Maciejewski (2020).Pearson (2019) uses a simultaneous 2nd order detrending polynomial but this is is not explicitly constructed to be depth-conserving.They do, however, use LDTk in their handling of LDCs, but not to directly inform the likelihood of trial parameters.Instead, Pearson (2019) uses LDTk to find the highest likelihood LDC values for the host parameters and TESS filter and fixes the values here.We note however that the host parameters used by Pearson (2019) for this do not exactly match those of (Maxted et al. 2016), as indicated in Table 5, and it is unclear where the alternative value of host metallicity originates from.Maciejewski (2020) fit physical parameters to the light curves after independently detrending the raw data.We note that they freely fit the LDCs without using the Kipping (2013) parameterisation built into TransitFit, and do not use any host information to inform the likelihoods.This is reflected in the larger LDC uncertainties.Maciejewski (2020) compares their final LDC values to those predicted by bi-linearly interpolating the LDC tables provided by Claret & Bloemen (2011) for the Cousins R and I bands and the Sloan Digital Sky Survey z band and then averaging the results to approximate the TESS filter, which suggests values of  0 = 0.30 and  1 = 0.28.These predicted values are closer to the results of the TransitFit "independent" LDC run, which most closely resembles the approach of Maciejewski ( 2020), but are in significant tension with the results from the "coupled" LDC run.
We therefore conclude that, by taking a different approach to the investigations of both Pearson (2019) and Maciejewski (2020), we have been able to independently verify the result of Maciejewski (2020) and find no statistically significant TTV signatures and, consequently, no evidence of a second planet in the WASP-126 b system.

HST observations of WASP-43 b
WASP-43 b, discovered by Hellier et al. (2011), is a 2.034 ± 0.052  Jup hot Jupiter with a radius of 1.036 ± 0.019  Jup , orbiting a type K7V star with a period of 0.813 days (Gillon et al. 2012).The proximity of the planet to the host star makes it an ideal candidate for emission spectroscopy, and it is one of the few exoplanets to have observed emission phase curve data (Stevenson et al. 2014).Other studies have shown that WASP-43 b possesses a strong day-night temperature contrast (Stevenson et al. 2014;Gandhi & Madhusudhan 2018) and consequently strong equatorial jets in the atmosphere (Kataria et al. 2016).
Since the first detection of an atmosphere by Murgas et al. (2014), multiple transmission spectroscopy studies have been completed, with detections including water (Kreidberg et al. 2014), carbon monoxide, carbon dioxide (Feng et al. 2016), aluminium oxide (Chubb et al. 2020), and multiple hydrocarbon hazes (Helling et al. 2020).Emission spectroscopy studies have also been utilised (Stevenson et al. 2014;Kreidberg et al. 2014;Blecic et al. 2014;Stevenson et al. 2017), and from these it has become apparent that a difference in the abundance of water exists between the day and night sides of the planet.
Several of these studies, including those of Tsiaras et al. (2018) and Chubb et al. (2020), make use of the Kreidberg et al. (2014) observations from HST, taken as part of observing program 13467 (Bean 2013).These observations were taken using the WFC3 G141 grism over the wavelength range of 1.135-1.642m, and include three full-orbit phase curves, three primary transits, and two secondary eclipses.In this section, we demonstrate the application of TransitFit to the observations in the 1.135-1.158m waveband.We limit ourselves to this single waveband in this paper, as we are using TransitFit to conduct an in depth study of observations of WASP-43 b from a wide range of sources (SPEARNET, in prep).
The observations which include a transit are shown in the top plot in Figure 6, normalised to a median value of 1.For the purposes of this discussion, we distinguish between visits and orbits.A visit is a single one of these observations, and can be directly translated to an epoch within TransitFit.Since HST is an orbital satellite, there are times during a visit where the source cannot be seen, which leads to a temporal striping effect.Each one of these stripes is referred to as an orbit.We also define  visit , the time elapsed since the first exposure in a visit and  orb , the time elapsed since the first exposure in an orbit.
It can clearly be seen that, in addition to the observation-long trends, there are two complex systematics which affect the observations.These are the alternating offset introduced by the upstream/downstream effect (McCullough & MacKenty 2012), and ramp-up effects in observations made in each orbit (Berta et al. 2012;Agol et al. 2010).Following Kreidberg et al. (2018), we assume that the observed flux,  obs is given by where  sig is the astrophysical signal and  sys is the signal from the systematics that need to be removed.We model these using where  = 1 for forward scans  for reverse scans (30) and ,  1 ,  2 , , and  are all detrending coefficients.We note that Kreidberg et al. (2018) multiply  by a normalisation coefficient, but we exclude this as it would be degenerate with the normalisation constant fitted by TransitFit.
Implementing the detrending model in Equation 29as a custom detrending model within TransitFit, we analyse the observations using the folded mode , "coupled" LDCs, and priors centred around (Kreidberg et al. 2014) results.The phase-folded lightcurve, residuals, and binned residuals are shown in the bottom plot in Figure 6, where it is clear that the systematics have been removed.In Table 6 we give the retrieved orbital parameters alongside those of Kreidberg et al. (2014).The temporal results are broadly consistent, but there is a discrepancy between the retrieved values of  and   .We suggest that this could be caused by the slightly different approach used by Kreidberg et al. (2014).The Kreidberg et al. (2014) detrending model is similar to our Equation ( 29), but excludes the quadratic  2  2 visit term from the visit-long trend, and it is not fitted simultaneously with their transit model.Additionally, there are only two free parameters used in the Kreidberg et al. (2014) transit model;   / ★ and a linear LDC.Their period is fixed to that of Blecic et al. (2014), whilst their values for , , and  0 are determined by fitting white-light observations.In our analysis all of the transit and de-trending parameters are fitted simultaneously by TransitFit.Further investigation into these effects are being made as part of the aforementioned SPEARNET study, but we have included this section to demonstrate that TransitFit is capable of fitting light curves that exhibit systematics as complex as those seen in observations made by HST.

Combined fitting of JWST and HST observations of WASP-96 b
WASP-96 b is a hot Jupiter with a radius of 1.20 ± 0.06  Jup and mass of 0.48 ± 0.03  Jup orbiting a G8 type star with a period of Table 5.The planetary and orbital parameters of WASP-126 b derived using TransitFit in both "coupled" and "independent" LDC fitting modes on observations from sectors 1-13, 27-34, 36-39, 61, and 63 of the TESS mission.We present these alongside the values found by Pearson (2019), Maciejewski (2020), and Maxted et al. (2016).All fits assume zero orbital ellipticity, use a quadratic limb-darkening model, and unless otherwise stated assume the same host parameters as given by Maxted et al. (2016).A total of 94,021 datapoints were fitted, and it gave a value of   = 1.17 × 10 −04 and  2 = 89, 956.47 for the "coupled" mode; and   = 1.21 × 10 −04 and  2 = 89, 898.40 for the "independent" mode.
.17 ± 0.08 a These values were derived assuming that no TTVs were present.b These values are those predicted by LDTk, using host parameters from Maxted et al. (2016) with [Fe/H] = −0.06.c The value of  = 2.8493819 ± 0.0000013 days provided in Maciejewski (2020) appears to be a typo as it exactly matches the period given for WASP-100 b. d The time standard is not specified.
3.4 days (Hellier et al. 2014).It has been a target for several atmospheric studies.An analysis of VLT FORS2 observations by Nikolov et al. (2018) showed an abundance of sodium in its atmosphere while predicting a cloud-free atmosphere at the limb.This analysis was reiterated in Nikolov et al. (2022) while confirming oxygen signatures in its atmosphere using HST WFC3 and Spitzer IRAC data.Further, McGruder et al. (2022) also confirmed that the WASP-96 b atmosphere showed minimal aerosol content at its terminator.However, Samra et al. (2022) has predicted a presence of a cloudy atmosphere following an analysis of VLT, Spitzer, and HST data.The atmosphere of WASP-96 b was shown to have features of water in its spectrum by Yip et al. (2020) using HST observations, while also inferring a lack  Kreidberg et al. (2014).These data were all obtained using the 1.135-1.158m waveband of the WFC3 G141 grism.The data from the first orbit of each visit and the first exposure of each orbit have been discarded and are not displayed or used in the subsequent analysis.The data are displayed here in raw form, normalised to a median value of 1 and phase folded using the TransitFit-retrieved ephemeris given in Table 6.The scan direction used when obtaining the data introduces a systematic offset between the forward scan data (in dark grey) and the reverse scan data (in light grey).Bottom: The phase-folded light fitted curve obtained from the above data using TransitFit with a wavelength-coupled, quadratic limb darkening law and the custom detrending function given in Equation 29.The error bars have been scaled using the retrieved value of   .The best-fit light curve is over-plotted in black and the residuals are also displayed.We also report the  2 corresponding to the fitting for a total of 469 datapoints.(2023) and JWST data analysis tools9 .The images were corrected for bad pixels identified by DQ flags, by replacing them with the median value of the 9x9 pixel grid centred on them.After applying a flatfield correction, we traced the order 1, order 2, and order 3 spectra.The trace was then smoothed out by fitting it with a Chebyshev polynomial.The spectra of all three orders were masked out within 20 pixels from their traces and a noise and background correction was applied by removing the column median of the residual image.The traces of the spectra and the background-corrected versions are shown in Figure 7.With an aperture radius of 15 pixels, we extracted lightcurves from order 1 in 850-2800 nm range and from order 2 spectra in 600-850 nm range only, leaving out the noisier order 3 spectrum.The wavelength integrated lightcurve of this spectrum had a median noise of 720.1 ppm as compared to 77.0 ppm for order 1 and 164.0 ppm for order 2. The lightcurves were binned in wavelength bins of 25 nm, using inverse variance weighting.
HST data from G141 and G102 grisms were extracted using Iraclis 10 (Tsiaras et al. 2016b; Tsiaras et al. 2016c) in bins similar to that of Yip et al. (2020).A total of 88 lightcurves obtained from JWST observations and 38 lightcurves from HST observations were fitted using TransitFit to generate a transmission spectrum.We use a 2nd-order detrending function for JWST lightcurves and a custom detrending function for HST lightcurves as with the analysis of WASP-43 b described in Section 3.4.The limb darkening coefficients were fitted in the coupled mode and using quadratic limb-darkening model.Planetary parameters and host parameters from Hellier et al. (2014) were used as priors and inputs, respectively, and are presented in Table 7.The "batched" mode was used for fitting these lightcurves, and the average of results weighted by inverse variance, were generated.Given the large number of lightcurves involved, there is a possi- bility of inter-batch variability in the wavelength-independent parameters.Consequently, the wavelength-dependent parameters might not result in the best-fit when used with the final results for wavelengthindependent parameters.In order to reduce this discrepancy, we run the "batched" mode again.To generate priors for this run, we take the results from the first run, and calculate inverse variance weighted results from all batches leaving one batch at a time.The union of these results is taken as the prior for the wavelength-independent parameters in the second run, which gives us the final results.
Figure 8 shows the HST and JWST lightcurves after fitting , along with the best-fit model, with the best-fit values listed in Table 7, while the generated transmission spectrum is shown in Figure 9.The TransitFit spectrum appears to be in good agreement with the analysis of Yip et al. (2020).

CONCLUSIONS
We have presented TransitFit, a new open-source code for fitting exoplanetary transit light curves using nested sampling routines.Tran-sitFit has been designed for transmission spectroscopy surveys employing multiple telescopes, and allows coupling of limb-darkening coefficients across observation wavelengths by utilising information on the host star and the LDTk Python package (Parviainen & Aigrain 2015).
TransitFit has been developed in anticipation of a new "assetstarved" era of transmission spectroscopy studies, where limited observational time and resources mean that studies will frequently have to combine data of various quality, wavelength coverage, and sources.One such example of this is SPEARNET, a survey which is using a heterogeneous distributed network of small-to mid-sized ground-based telescopes to conduct atmospheric studies of transiting exoplanets.
Using TransitFit and observations from the SPEARNET telescope network, we have presented analysis of new data of the hot-Neptune WASP-127 b, which includes the first  ′ -band observations of the planet.We have shown that introducing a wavelength-coupled approach to LDC fitting can result in changes as large as 8 per cent in the retrieved value of   / ★ , or 17 per cent in measured transit depth.
We have demonstrated the application of TransitFit in more temporal-focused studies, analysing TESS observations of 26 transits of WASP-91 b to produce updated planetary ephemerides.This will prove invaluable in analysis of planets in the TESS catalogue, allow- ing for easy searches for TTV signatures.We have used TransitFit to analyse 180 transits of WASP-126 b observed by TESS and have found no statistically significant evidence for the presence of TTVs proposed by Pearson (2019).
We have also shown how TransitFit can be used in situations where observations display trends that cannot be adequately modelled with a low-order polynomial.We fitted observations of WASP-43 b from a single wavelength channel of HST and found that the complex systematics can be adequately removed.Further analysis of the full set of HST observations for WASP-43 b, in conjunction with observations from other telescopes, is being conducted as part of a separate study (SPEARNET, in prep).Moreover, we have already used TransitFit to analyse the combined ground-based, HST and TESS observations of HAT-P-26 b (A-thano et al. 2023).The analysis found the presence of TTVs and H 2 O dominated atmosphere on the planets, which agrees with previous studies.
We used TransitFit to construct a combined transmission spectrum of WASP-96 b from simultaneous fitting of JWST NIRISS and HST WFC3 data.In general we found good correspondence between the HST and JWST datasets.
(b).The corresponding corner plot showing posteriors are shown in Figure A1 in Appendix A. We also show the resulting transit depths from both the coupled and independent LDC mode retrievals alongside the Hubble Space Telescope (HST) spectrum from Skaf et al. (2020), and the Gran Teliscopio Canarais (GTC) and Nordic Optical Telescope (NOT) spectral data and best fit model from Chen et al. (2018) in Figure 3. Since our data, and those of Chen et al. (2018) and Skaf et al. (

Figure 1 .
Figure 1.The location and size of the telescopes currently within the SPEARNET network.Telescopes with a red circle are those used to take the observations of WASP-127 b discussed in Section 3.1.(Image credit: NARIT)

Figure 2 .
Figure 2. (a): The raw data from the SPEARNET observations of WASP-127 b outlined in Table1, normalised using the median values for each observation and phase folded using the model obtained in (b).We have binned the observations to a 2-minute cadence, with the raw observations shown in grey.For clarity, error bars on data points have been excluded and an average error bar of the raw data points has been provided as a black coloured bar to the right of each curve.(b): The detrended light curves and associated best-fit models obtained using TransitFit in "batched" mode with a coupled quadratic LDC fitting and simultaneous normalisation and 2nd-order detrending.Each curve has an arbitrary offset from a normalised baseline of 1 and has been phase-folded to centre  0 at a phase of 0.0.The best fit transit model from the retrieval is over-plotted.We show the average of the errors scaled using retrieved value of   , as a black coloured point on the right of each curve.(c): The residuals after model is subtracted from the fitted lightcurves.(d): The histogram of the residuals along with the respective  2 values.
13 a The time standard is not specified .bAnderson et al. (2017) use the non-linear limb-darkening law but do not provide coefficients to compare with.

Figure 3 .
Figure 3.The transmission spectrum of WASP-127 b, obtained by fitting ground-based observations from the SPEARNET telescope network using TransitFit in coupled (orange) and independent (blue) LDC fitting modes with a second-order detrending function.From left to right the filters are  ′ ,  ′ , , and  ′ , and we have labelled the independent LDC depths for quick reference.Also shown are the GTC/NOT observations (after applying corrective offset) reported in Chen et al. (2018, light grey) and the HST reported in Skaf et al. (2020, dark grey).For the SPEARNET data points, the horizontal positions and error bars are derived from the SED-weighted filter profiles, with the data point positioned at the weighted mean of each profile and the error bars reflecting the standard deviations.For the other data sources, horizontal error bars represent the complete bandwidth of each spectroscopic filter.

Figure 4 .
Figure 4.The full data processing of the TESS observations of WASP-91 b using the "folded" retrieval mode of TransitFit, with LDCs informed by LDTk.a(i): the raw observations from TESS Sector 1; a(ii): the raw observations from Sectors 27 and 28.We have excluded the error bars from these two plots for clarity; b:The raw data from the first observed transit, chosen purely as an illustration; c: The first observed transit after the first stage of "folded" mode retrieval.These data have been detrended and normalised, the error bars have been scaled using the retrieved   , and have the final best-fit model overplotted; d: the final phase-folded light curve, best-fit transit model and residuals, containing 9 789 photometric data points, and, in red, the same data binned to a 2 minute observation cadence.The error bars have been scaled using   , and the corresponding  2 value has also been displayed.

Figure 5 .
Figure 5. Left: The O-C plots of 126 transits of WASP-126 b observed by TESS, obtained from results of using TransitFit in "coupled" (top) and "independent" (bottom) LDC modes whist allowing for the presence of TTVs.The values of  0 in this case had a median error of 4.97×10 −4 in the coupled mode and 4.90×10 −4 in the independent mode.The dotted lines represent the uncertainty on the predicted ephemerides derived from using TransitFit without allowing for TTVs.Right: The associated weighted Lomb-Scargle periodograms for each of the O-C plots.The horizontal grey lines represent the false alarm probabilities of 10, 5, and 1 per cent from bottom to top.

Figure 6 .
Figure 6.Top: Observations of six transits of WASP-43 b taken using the Wide Field Camera 3 (WFC3) on HST byKreidberg et al. (2014).These data were all obtained using the 1.135-1.158m waveband of the WFC3 G141 grism.The data from the first orbit of each visit and the first exposure of each orbit have been discarded and are not displayed or used in the subsequent analysis.The data are displayed here in raw form, normalised to a median value of 1 and phase folded using the TransitFit-retrieved ephemeris given in Table6.The scan direction used when obtaining the data introduces a systematic offset between the forward scan data (in dark grey) and the reverse scan data (in light grey).Bottom: The phase-folded light fitted curve obtained from the above data using TransitFit with a wavelength-coupled, quadratic limb darkening law and the custom detrending function given in Equation29.The error bars have been scaled using the retrieved value of   .The best-fit light curve is over-plotted in black and the residuals are also displayed.We also report the  2 corresponding to the fitting for a total of 469 datapoints.

Figure 7 .
Figure 7. Top: The traced spectra for Order 1 (red), 2 (blue), and 3 (brown) of the GR700XD grism are shown overplotted on the median image from WASP-96 b observations.Bottom: The median of the image after noise correction.

Figure 8 .
Figure 8. Left column (a): a sample of fitted lightcurves from HST observations of WASP-96 b, with the best-fit model overplotted in black.The error bars on all data points are plotted in grey and are scaled using the retrieved value of   .The corresponding wavelength bin for the lightcurves is indicated.The right column (b) shows the same for JWST data.TransitFit computed the best-fit model from both HST and JWST datasets simultaneously.The  2 values for JWST lightcurves were in the range of 269.38-294.91 with an average of 280.08.For HST lightcurves, the  2 values were in the range of 38.07-96.33 with an average of 53.89 for the data from G141 which had 44 data points.The data from G102 had 40 integration points and resulted  2 values in the range of 38.46-81.09with an average of 59.22.

Figure A1 .
Figure A1.A posterior plot showing distribution of parameters samples across the range of priors.We use gaussian priors for ,  0 , and / ★ ; and uniform priors for   / ★ .The priors for inclination angle () are from a gaussian distribution clipped at a maximum of 90 degrees.These samples were taken for WASP-127 b fitting in 'coupled' mode, and are plotted without weights.The LDCs, normalisation factor, and detrending coefficients are not shown for simplicity.The best fit values and corresponding error limits are shown as the titles of histograms.

Table 2 .
The photometric data of WASP-127 b.

Table 6 .
(Pontoppidan et al. 2022) parameters of WASP-43 b obtained by TransitFit using the folded mode and wavelength-coupled quadratic LDC fitting on 1.135-1.158mHSTobservations,alongsidethevaluesusedand obtained byKreidberg et al. (2014).The only free parameters in theKreidberg et al. (2014)model are   / ★ and a linear LDC, the final value of which is not quoted.3 in its atmosphere.As part of Early Release Observations of the JWST, WASP-96 b was observed in NIRISS Single Object Slitless Spectroscopy (SOSS) mode consisting of 280 integration points, to demonstrate the capabilities of the instrument(Pontoppidan et al. 2022).In this section, we show the transmission spectrum of WASP-96 b generated from a TransitFit analysis applied simultaneously to this JWST data and to archive WFC3 HST data.For the extraction of the lightcurves, the Stage 1 images from the NIRISS observations were processed following Feinstein et al.

Table 7 .
Hellier et al. (2014)ital parameters of WASP-96 b derived using TransitFit analysis of the JWST mission.The values are compared with the values retrieved fromHellier et al. (2014), which were also used as priors for TransitFit.This value is in UTC.The corresponding prior for  0 was checked to be in-transit for the raw lightcurve. a