FRECKLL: Full and Reduced Exoplanet Chemical Kinetics DistiLLed

We introduce a new Python 1D chemical kinetic code, Full and Reduced Exoplanet Chemical Kinetics distiLLed (FRECKLL), to evolve large chemical networks efficiently. FRECKLL employs “distillation” in computing the reaction rates, which minimizes the error bounds to the minimum allowed by double precision values (ϵ ≤ 10−15). Compared to summation of rates with traditional algorithms like pairwise summation, distillation provides a tenfold reduction in solver time for both full and reduced networks. Both the full and reduced Venot2020 networks are packaged in FRECKLL as well as a TauREx 3.1 plug-in for usage in forward modeling and retrievals of exoplanet atmospheres. We present TauREx retrievals performed on a simulated HD 189733b JWST spectra using the full and reduced Venot2020 chemical networks and demonstrate the viability of total disequilibrium chemistry retrievals and the ability for JWST to detect disequilibrium processes.


INTRODUCTION
In the last decade, observations from space using mainly the Hubble Space Telescope (HST) and the Spitzer Space Telescope (Spitzer) and from the ground, have allowed to characterise the atmospheric properties of a handful of planets from their transit (Tinetti et al. 2007;Kreidberg et al. 2014;Sing et al. 2016;Sedaghati et al. 2017;Wakeford et al. 2017;Tsiaras et al. 2018;Fisher & Heng 2018;Anisman et al. 2020;Edwards et al. 2021;Gressier et al. 2022;Wong et al. 2022;Saba et al. 2022;Edwards et al. 2022), eclipse (Swain et al. 2008;Crouzet et al. 2014;Line et al. 2014;Haynes et al. 2015;Line et al. 2016;Edwards et al. 2020;Changeat & Edwards 2021;Changeat et al. 2022;Fu et al. 2022) or phase-curve observations (Stevenson et al. 2017;Arcangeli et al. 2019;Changeat et al. 2021;Changeat 2022;Mikal-Evans et al. 2022;Chubb & Min 2022).Due to Corresponding author: Ahmed Faris Al-Refaie ahmed.al-refaie.12@ucl.ac.uk * Paris Region Fellow the low resolution and narrow wavelength coverage of older generation space-based instrumentation, however, degeneracies can often lead to multiple interpretations of exoplanet spectra, depending on model and prior assumptions (e.g.Changeat et al. 2020b).To explore the information contained in these spectra, exoplanet teams have developed sophisticated methods to invert the information content in the spectra of exoplanets.These methods, often called spectral retrieval techniques (Irwin et al. 2008;Madhusudhan & Seager 2009;Benneke & Seager 2012;Line et al. 2013;Waldmann et al. 2015;Min et al. 2020;Mollière et al. 2019;Al-Refaie et al. 2021b) require the evaluation of thousands to millions of forward models, therefore requiring significant computing resources.Often, the computing requirements imply that simplified atmospheric models have to be employed, for instance, by assuming 1-dimensional geometries and other idealized assumptions on the thermal structure, the chemistry and the cloud properties.Since the information extracted from current spectra is low, assumptions are commonly used throughout the literature.These assumptions include isothermal ther-mal structure, constant chemical profiles or equilibrium chemistry, and fully opaque cloud opacities.
In retrieval codes, the chemistry is often recovered using profiles constant with altitude; a single free parameter representing each molecule (e.g., see references above).While not representative of an entire atmosphere, current observations mostly probe small pressure regions where chemical variations remain small.An alternative assumption is thermochemical equilibrium (White et al. 1958;Eriksson et al. 1971), which requires computing the chemistry state by minimizing the Gibbs free energy of the system.Such an assumption has gained popularity due to the reduced degrees of freedom.Furthermore, it often only requires two free parameters for metallicity and the C/O ratio chosen for their natural links to planetary formation and evolution processes.Equilibrium chemistry, however, is a strong assumption with little justification.Given with Tsai et al. (2023) and Dyrek et al. (2024), in which SO 2 , produced by photolyses, has been detected, equilibrium does not represent the underlying chemical processes being detected in exoplanet atmospheres.Furthermore, simulations employing kinetics methods and thus taking into account disequilibrium processes such as mixing and photochemistry have proven that chemical equilibrium is inadequate in many scenarios (e.g.Moses et al. 2011Moses et al. , 2013Moses et al. , 2016;;Venot et al. 2012Venot et al. , 2014Venot et al. , 2020a,b;,b;Molaverdikhani et al. 2019;Tsai et al. 2021;Morley et al. 2017;Mollière et al. 2020;Kawashima & Min 2021).With future telescopes, accurate representation of the chemical processes will be essential to ensure unbiased interpretation of the observations as discussed in the first analyses of JWST data (The JWST Transiting Exoplanet Community Early Release Science Team et al. 2022;Tsai et al. 2023;Dyrek et al. 2024).
In this paper we present the first implementation of a full chemical kinetic scheme into an atmospheric retrieval framework.We use the flexibility of the plugin system in TauREx 3.1 to integrate this new scheme and explore the use of chemical kinetic models in atmospheric retrievals.In particular, we focus our study on quantifying the impact of the equilibrium chemistry assumption in interpreting atmospheres exhibiting disequilibrium processes.Section 2 presents our implementation of the chemical kinetic code and the steps carried out in this work.In Section 3, we present the results of our simulations.Finally, Section 4 discusses our findings and provides the main conclusions of our exploration.

Description of chemical kinetic model
As opposed to thermochemical equilibrium models, which predict the chemical state of a planet's atmosphere by minimising the Gibbs free energy of the system, chemical kinetic models necessitate integrating the system of differential equations representing each considered reaction until a steady state is reached.The continuity equation (Equation 1) describes the temporal evolution of the abundance of each species i, considering a one-dimensional plane-parallel atmosphere.
where n i , P i and L i are the number density (cm −3 ), production rate and loss rate of species i (cm −3 .s−1 ), z is the vertical coordinate of the atmosphere, and ϕ i is the vertical flux for species i which has the form of a diffusion equation given in Equation 2.
Here, D i is the molecular diffusion coefficient (cm 2 s −1 ), H i is the scale height (km) and y i the mixing ratio for species i, T is the temperature (K) and K zz is the eddy diffusion coefficient (cm 2 s −1 ).While, rigorously, thermal diffusion should be included in this equation (Drummond et al. 2016), it has been found that it was negligible compared to the other terms of Equation 2 (Venot 2012).We thus don't include it.
At t = 0 s, an initial abundance is set.FRECKLL can accept any initial atmospheric abundance, either user-supplied or from an external code.As the default, the system is initialised with the abundance of each species assumed to be at thermochemical equilibrium.This initial state is computed using the ACE code (Agúndez et al. 2012) with supplied or user-defined NASA polynomial thermochemical coefficients and, subsequently, Equation 1 is evolved using a stiff ODE solver such as VODE (Brown et al. 1989) or DLSODES from the ODEPACK (Hindmarsh 1983) package until steady state is achieved or a user-defined condition is reached.Metallicity, elemental ratios (e.g.C/O and N/O ratios) can be set to determine the initial abundance produced by ACE.In addition to those parameters, the model also allows for the definition of the eddy diffusion parameter K zz , given as a constant value or layer-by-layer.For our test cases, we employ the full Venot et al. (2020a) network.This network is based on the Venot et al. (2012) network with updates to its methanol chemistry and includes 108 species, 1906 reactions and 55 photodissociations (see Table A.1).The network, including all the photolyses data can be found on the website of the ANR EXACT 1 .In addition, Venot et al. (2020a) provides a reduced network consisting of 44 species and 582 reactions and omits photolysis reactions.
For a 130 layer atmosphere, FRECKLL takes roughly 4-5 minutes to reach a steady-state, using the full chemical network and roughly 30 seconds on the reduced, significantly speeding up convergence.

The importance of Numerical stability
Floating points, being approximations of real numbers, inherently carry errors with each operation.Numerical stability involves understanding, managing, and minimizing these errors to ensure the accuracy and reliability of computations.This problem is not exclusive to Python or the particular packages highlighted in this paper.These challenges permeate various computational platforms, even extending to compiled languages like FORTRAN.Chemical kinetics, by nature, presents inherent stiffness in its equations, stemming from the wide-ranging magnitudes of chemical timescales and abundances.Integration requires stiff ODE algorithms such as Backwards differentiation formula, Rosenbrock and backwards Euler methods which can vary the time steps over large orders of magnitude.Additionally, an overlooked aspect involves the computation of sums.With double precision we generally expect the upper bound of relative errors from rounding to be ϵ m ≈ 10 −16 .Assuming a function implemented with algorithm f (x) and the true function f (x) the relative error for an algorithm ϵ is computed as: Here, the true function is f (x) = n x i where summation is performed at infinite precision.We must also consider the condition number C: which represents the intrinsic sensitivity of summation.
For naive summation such as the inbuilt python sum function, the error is bounded as: where n is the number of elements.The error for pairwise summation used by the numpy.sumfunction is bounded by: 1 https://www.anr-exact.cnrs.fr/Generally, well-conditioned problems are those where C ≈ 1, one such case is where all values are non-negative (i.e x i > 0).For 10,000 elements, the error from naive summation is ϵ ≤ 10 −12 and for pairwise we expect an error ϵ ≤ 10 −15 .
The problem comes when dealing with extensive magnitudes and a mixture of negative and non-negative values.To illustrate, let us take an array of values x = [10 16 , 10 30 , 1, 5, 10, 10 4 , −10 30 , −10 16 ], where x = 10016.Attempting to use the native sum we get: >>> sum ( [1e+16, 1e+30, 1, 5, 10, 10000.0,-1e+30, -1e+16]) -7638326771712.0 the error is in the order of ϵ ≈ 10 8 .Pairwise summation performs a little better: >>> numpy.sum( [1e+16, 1e+30, 1, 5, 10, 10000.0, -1e+30, -1e+16]) 0.0 here ϵ = 1.The problem is ill-conditioned with a condition number of C = 10 26 which is extremely large.This arises from catastrophic cancellation where precision limits for floating points mean dbl(10 30 + 1) = 10 30 where dbl is an operation under double-precision.For kinetics calculations, this is problematic, as summing production and loss rates for a molecule can suddenly become zero, change sign or magnitude.In the Jacobian, these appear as sudden discontinuities and can cause stiff ODE methods to oscillate at certain times and continually reduce the time-step, which will drop the integration efficiency.To avoid these problems in FRECKLL, we instead employ the K-fold summation method (Ogita et al. 2005).This method first performs an error-free transformation of an array: The twosum computes the resultant floating point sum and residuals from the summation.For each element i, the result is stored at i and the residual at i − 1.For an array x the algorithm produces a resultant array y where y = vecsum(x) which has the property: assuming infinite precision.This is often referred to as 'distillation' (Kahan 1987).For low condition numbers, elements i = 1...n − 1 will be zero and i = n will contain the resultant sum (i.e y n = n x i ).For higher condition numbers this takes the form n−1 y i + y n = x i .Distillation has the effect of reducing the condition number and indeed, applying it to our original array x = [10 16 , 10 30 , 1, 5, 10, 10 4 , −10 30 , −10 16 ], the condition number falls to C ≈ 10 9 .As distillation preserves the original sum and reduces the condition number, we can apply it K −1 times until we reach our desired condition number before applying the summation; this is K-fold summation in its essence: >>> kfold ( [1e+16, 1e+30, 1, 5, 10, 10000.0,-1e+30, -1e+16], K=2) 10016.0Applying this algorithm for K = 2 we indeed get the correct result.The error bounds for K = 2 are ϵ K=2 ≤ 10 4 .Increasing to K = 4 gives an error bounded ϵ K=4 ≤ 10 −15 which is the maximum possible with double precision.K-fold summation is significantly slower than numpy.sum;10,000 elements takes roughly 7-10x longer than numpy.However, as we will demonstrate, the increase in precision greatly benefits convergence.We employ K-fold summation in computing the production and loss rates of molecules.To maximise computational efficiency and precision we combine the rates into a single array.For molecule i and reaction r we combine the production rate P i r and loss rate L i r into a total molecule rate R i .If we have p production reactions and l loss reactions then: The total rate for molecule i is given as: where k is our K-fold function.We can rewrite Equation 2as: Additionally we also include convergence criteria to numerically determine if steady-state has been reached.The criteria is the same as the one used in VULCAN (Tsai et al. 2017).Given timesteps m and m + 1 we compute the following: Our criteria for steady state is therefore: where δ and η are our criterion parameters for the relative change and relative change over time, respectively.We demonstrate the benefit of K-fold summation by solving a benchmark system.We compute a HD 209458 b model between 10 −5 -10 2 bar, using the thermal and vertical mixing profiles from Venot et al. (2020a) displayed in Figure 1.The model consists of 130 layers, 108 molecules, 1906 reactions and 55 photodissociations.For the actinic flux, as HD 209458 is a G0 star, we use the UV spectral irradiance of the Sun (Thuillier et al. 2004) scaled to correspond to the radius and effective temperature of HD 209458.We evolve the system using the VODE solver with a relative tolerance of 10 −3 and absolute tolerance at 10 −25 until t = 10 10 s with steady-state occurring at t = 10 8 s.For this case we set δ = 0 and η = 0, disabling the criteria.We will also assess the number of function evaluations (evaluation of Eq. 1) and the number of times the jacobian matrix is evaluated.
Figure 2 shows the initial abundances at equilibrium computed using ACE and the final steady state solution achieved by FRECKLL at t = 10 10 s.Using pairwise summation implemented by numpy.sumtakes roughly 128 minutes to evolve until t = 10 10 s requiring 467335 function evaluations and 2179 jacobian evaluations.The issue is that the solver is unable to take larger time steps, especially at t ≈ 10 8 s where ∆t ≈ 10 4 s.As discontinuities appear more often, the solver is forced to use small timesteps in order to ensure smoothness in the function.This effect becomes more pronounced as the system approaches steady state as the solver has difficulty integrating | ∂ni ∂t | below a certain threshold.Tracing this during integration estimates this threshold to be around 10 −10 .Solving the same system using K-fold summation (K = 4) until t = 10 10 s takes 5 minutes, 2,682 function evaluations and 158 jacobian evaluations.As we stated previously, K-fold summation is significantly slower than piecewise summation but we manage to gain a 25x reduction in solver time as well as a 175x reduction in function evaluations.The improved precision means that | ∂ni ∂t | can reach 10 −15 and the solver is choosing larger time-steps that skip from 10 8 s−10 10 s.In fact, solving further to 10 12 s takes only 10 extra function evaluations.
There is always a trade-off between raw performance and precision.When dealing with stiff non-linear systems, convergence can be hampered by the underlying precision of algorithms.It is sometimes easy to forget that summation is also an algorithm and not an intrinsic feature of computation.We demonstrate that choosing a slower, more precise summation algorithm can lead to significant performance gains from faster convergence.We also like to note that the summation algorithm presented here is not exclusive to chemical kinetics and can be applied to other ill-conditioned problems.

FORWARD MODELS
FRECKLL includes a plugin for TauREx 3.1 (Al-Refaie et al. 2021b,a) for generation of synthetic spectra and retrievals using the chemical kinetic code.We demonstrate its forward modelling capabilities by simulating HD 189733 b with parameters taken from the literature which are given in Table 1.Note that we simulate HD 189733 b with a constant K zz of 4×10 8 cm 2 s −1 .This is a fairly large value for the eddy diffusion coefficient, suggesting strong vertical mixing in the atmosphere in this scenario.For simplicity, the temperature profile for those simulations is modelled using an isothermal profile, even if GCM models predict variations with altitude, as well as with longitude and latitude (e.g.Drummond et al. 2020).In the model, we include absorption using the ExoMol line-lists (Tennyson & Yurchenko 2012;Chubb et al. 2021;Tennyson et al. 2020) Cox (2015).The atmosphere is modelled in plane-parallel geometry with 100 layers spaced between 10 bar and 10 −5 bar in log space.The actinic flux used for HD 189733 is the same as the one used in Venot et al. (2012), originally produced by Ignasi Ribas (priv.comm.), it comprises of spectra from X-exoplanets (Sanz-Forcada et al. 2011), FUSE and HST data of ϵ Eridani and PHOENIX data (Hauschildt et al. 1999) for the spectral regions 0.5-90 nm, 90-330 nm and 330+ nm respectively.We solve the kinetics with a relative tolerance of 10 −3 and we reduce the absolute tolerance to 10 −20 for speed without harming the precision of the retrieval as molecules below this density are not spectrally visible.We set the convergence criteria to δ = 10 −4 and η = 10 −4 and maximum integration time of t = 10 30 s.On a 2.3 GHz Quad-Core Intel Core i5, the single-core combined runtime (kinetics and radiative transfer) for the reduced and full networks are 23 seconds and 3.2 minutes, respectively.The chemical profiles from the reduced and full networks are presented in Figure 3 with corresponding transmission spectra given in Figure 4.
We observe large differences between the two chemical networks from the chemistry predictions.In particular, while the predictions at the bottom of the atmosphere are consistent, large differences in the predicted abundances can be seen for the top of the atmosphere (below 0.1 bar).Those differences are likely due to the inclusion of reactions for photochemistry in the full network.We note, in particular, that the abundances of CH 4 and NH 3 decrease very rapidly for pressures above 10 −3 bar, while the C 2 H 2 profile is significantly affected.This translates in large differences in the observed spectrum at the wavelengths that are probing those altitudes.For  instance, the 2.3 µm and 3.6 µm methane features in Figure 4 and highlighted in Figure 5 are muted in the full network scenario (blue plot) compared to its reduced counterpart (orange plot) as this molecule is strongly photolysed in the upper atmosphere, with differences of the order of 200 ppm.As these differences are an order of magnitude greater than our simulated noise of JWST (Gardner et al. 2006) and the ≈ 20 ppm noise floors of Twinkle (Edwards et al. 2019) and Ariel (Tinetti et al. 2018(Tinetti et al. , 2021)), it is may be possible to infer, at least the presence of, methane photodissociation processes in the spectral data from these observatories.

RETRIEVALS
We now evaluate the performance and biases introduced (1) when using chemical kinetics rather than equilibrium and (2) when using two different chemical kinetics networks.The simulated spectra to be fit against in the retrievals make use of the same methodology as the preceding section but are convolved with the JWST instrument response for one transit of HD 189733 b with NIRISS GR700XD and one with NIRSpec G395M.The error bars are obtained using the ExoWebb instrument simulator (Edwards et al. in prep), which is based upon the radiometric model from Edwards & Stotesbury (2021).Normally HD 189733 would saturate the NIRISS instrument, which necessitates moving the star to 27 parsecs to prevent non-linearity in the detector response.We utilize the same priors for all cases described in Table 2.The K zz parameter, in particular, is fitted in the reduced and full networks cases to uniform priors between 10 3 -10 12 cm 2 .s−1 inclusive in log-space.For benchmark purposes and to provide comparisons with our previous works (Al-Refaie et al. 2021b,a), we high-

Reduced chemistry
In order to scrutinize the potential limitations and efficiencies of reduced chemical kinetic networks, we compared similarly to Venot et al. (2020a) but used a retrieval framework to characterize the applicability of a reduced chemical network.This test involved using a full chemical kinetic network with photodissociation disabled to simulate synthetic spectra.We then employed a reduced chemical kinetic network for the retrieval process.Our goal was to evaluate whether the simplified network could accurately capture information of the key reactions and constituents in the atmosphere, even with the inherent simplifications it encompasses.Figure 6 depicts the corresponding best-fit spectra.At a glance, the best-fit spectra match well with the simulated, which corroborates the similar results in Venot et al. (2020a).As illustrated in Figure B.1, the posteriors are well-defined and largely lie on or near the truth, with nominal values differing by 6-10%.It's expected that the posterior mass does not align exactly with the truth as, fundamentally, the models are not exactly the same.Observing the molecular profiles in per atmospheric pressures, respectively, likely from missing reactions.HCN is of note as it's one of the targeted molecules in Venot et al. (2020a).The abundances near the top of the atmosphere are higher in the full network, and that is likely due to the additional 29 HCNproducing reactions with many of their highest reaction rates in these layers for which photodissociation would act as a sink.Finally, for C 2 H 4 , its profile is retrieved with high precision.This is due to its involvement in the production of H and H 2 as either a reactant or in byproducts.In particular, its reaction pathways of reactions rates in the order of 10 5 cm −3 s −1 and tightly constrain its profile along the atmosphere.Overall, consistent with previous studies, the reduced network possesses similar chemical information on the composition of the atmosphere against the full chemical network without photodissociation.

Viability of full chemical kinetic retrievals.
Due to their long solve times, full chemical kinetic networks have generally not been used in retrievals.As FRECKLL significantly speeds up convergence, we will assess the viability of using chemical kinetic networks in retrievals.We will assess computational viability and the associated biases using such models in retrievals.In a realistic scenario, we would not possess full information about the atmosphere and its complex processes.However, we can replicate this by using the full network with photodissociation as a proxy for our complex atmosphere and utilize the full network in retrievals to represent perfect knowledge of the system and the reduced and equilibrium chemistry as a means to quantify how our assumptions influence the retrievals.This is the same methodology as Al-Refaie et al. (2021a) This is the first time a full disequilibrium kinetic retrieval, including vertical mixing and photochemistry, is attempted.Due to the large number of samples evaluated in atmospheric retrievals, numerical stability across the full range of parameters explored is key, highlighting the importance of the improvements described in the Methodology section.
The The observed spectrum as well as the results of these retrievals are shown in Figure 7.The kinetic runs took 8 and 24 hours to complete for the reduced and full networks, respectively.The equilibrium chemistry run took around 20 minutes.Posteriors are provided in Figure D.1, and the chemistry profiles are provided in Figure E.1.From the inspection of the best-fit spectrum, we observe that, as expected, the full network matches the observations.Additionally, the recovered free parameters are close to the chosen true value (see Figure D.1), and the recovered chemical profiles match the inputs within the uncertainties (see Figure E.1).Notably, the K zz has a well-defined posterior and is retrieved well by the full network on the simulated JWST observations, implying that disequilibrium processes are observable with JWST.This was also shown in previous works (e.g.Greene et al. 2016;Blumenthal et al. 2018;Molaverdikhani et al. 2019;Drummond et al. 2020;Venot et al. 2020b).The shape of the posterior displays a skew towards lower values, with a distinct boundary observed around 10 8 cm 2 .s−1 .This trend can be attributed to the fact that lower K zz values tend to favour chemical reactions, which generally yield equilibrium-like species profiles.Importantly, K zz predominantly affects the mid to upper atmosphere, where it takes precedence over reaction rates and whose change in chemical composition isn't as easily probed from the JWST spectra.This assertion is supported by the extensive range of the posterior, spanning over an order of magnitude and is further corroborated by the molecular profiles shown in Figure E.1.In these profiles, the uncertainties in the upper atmosphere are notably larger compared to those in the lower atmosphere, which are generally less responsive to variations in K zz .Compared to the retrieval in the previous subsection, photochemical processes (generally not dependent on temperature) are more directly coupled to K zz as their reaction rates depend on the number density of species transported into this region.The temperature profile also displays bi-modality, which suggests a degree of degeneracy, even with perfect knowledge of the atmospheric  processes.This bimodality is small with a difference between peaks of about 40 K around the truth value but suggests that higher resolution and/or lower errors on spectra are needed to fully remove degeneracies.
For the reduced network, the best-fit spectrum is at most wavelengths able to reproduce the observations, but we observe large discrepancies at certain bands.For instance, the 3.6 µm methane band is not well fitted by the reduced network, which, in our case, predicts too much methane.With an absence of photolysis reactions, the model lacks enough flexibility to compensate and deduce the correct abundance of CH 4 without affecting the other important molecules of the atmosphere.This is confirmed in Figure E.1.In most cases, the true value in this atmosphere is outside the 1σ predictions of the reduced retrievals.The metallicity, for instance, is found to be about Z = 6 when the input metallicity was solar (Z = 1).This is problematic as this could lead to incorrect interpretations, especially as such a parameter is commonly used to link atmospheric composition to planetary formation ( Öberg et al. 2011;Moses et al. 2013;Madhusudhan et al. 2016;Line et al. 2021).
We find similar results regarding the chemical equilibrium run.The retrieved parameters are most of the time outside the true values by more than 1σ.To compare the recovered metallicity again, assuming equilibrium chemistry, it is found to be about Z = 32.In general, we find that attempting to recover information content by modelling and simplifying complex processes in an atmosphere introduces strong biases in relatively unconstrained retrievals.The recovered chemical profiles, as seen in Figure E.1 present large departures from the input, with the main molecules being often different by more than two orders of magnitude.The predictions for both equilibrium and reduced runs are overconfident and do not reflect the raw information content in the spectrum.This is due to the assumptions (equilibrium chemistry or pre-selected list of reactions) introduced in those models that do not capture the essential physics in the atmosphere (Changeat et al. 2019(Changeat et al. , 2020a;;Al-Refaie et al. 2021a), in particular photodissociation.Similarly, one-dimensional retrievals, in general, may exhibit these types of biases as well (Feng et al. 2016;Caldas et al. 2019;Taylor et al. 2020;Changeat & Al-Refaie 2020;MacDonald et al. 2020;Skaf et al. 2020;Pluriel et al. 2020).While one might argue for the exclusive use of full chemical networks in retrievals as the most realistic, these models also face limitations.They depend heavily on the UV fluxes from the parent star and makes similar assumptions on what species are (and are not) present in the atmosphere.They are also highly dependent on our knowledge of chemical kinetics.Although considerable advances have been made recently (e.g.Veillet et al. 2024), certain reactions, and even certain couplings between elements (C-S or S-P for instance) are still insufficiently constrained.The approximations or assumptions made in chemical networks can, therefore, induce strong biases.In contrast, it has been shown that simple models such as constant profiles are better at broadly retrieving the chemical composition parameters such as Z and C/O (Al-Refaie et al. 2021a) as they more directly fit the spectral shape of species in the retrieval.A similar approach for K zz could be employed by using nlayer parameterized molecular profiles (Changeat et al. 2019), which could extract information on how species are distributed along the atmosphere.Emission spectra may also be employed in conjunction with transmission spectra to better constrain the temperature profile.Prospective work should look to understand better what information is needed to constrain particular classes of chemical models, especially when involving photodissociation.This is imperative in the era of JWST as it has already observed photolysis processes in exoplanets (Tsai et al. 2023).In summary, while complex models can theoretically provide a detailed understanding of atmospheric processes, their effectiveness is significantly diminished without sufficient data, especially in the context of wide and unconstrained retrievals.Adding information into the retrieval, such as constraining parameters or introducing more spectral data, will enhance a chemical model's ability to extract detailed information and reduce biases and degeneracies.

CONCLUSION
In this work, we introduce FRECKLL, a cuttingedge tool designed for the rapid and stable computation and retrieval of exoplanet chemical kinetics.Central to FRECKLL's efficiency is its distillation algorithm, which significantly enhances convergence to a steady state and allows for solving of large complex chemical networks in minutes.By integrating FRECKLL with TauREx 3 through its plugin system, we have for the first time successfully coupled chemical kinetics with retrievals, facilitating disequilibrium retrieval using a comprehensive kinetic network with photodissociation.We have shown that the use of strong assumptions about chemical composition (equilibrium, reduced or full networks) in retrievals could considerably bias the interpretation of observations, and we caution the reader in their use in exoplanet retrievals without significant constraints.However, this work paves the way for a new type of retrieval.If used with care, it could help improve our knowledge of exoplanetary atmospheres, particularly in the era of new telescopes beginning with the JWST.We understand the importance and benefits of open-source sharing within the academic commu-nity.However, as of the current state, FRECKLL is not yet available for public use.Our decision to delay the open-source release is grounded in ensuring the tool is user-friendly and free from potential pitfalls that might arise from its current intricacies.We're dedicated to refining the codebase, enhancing its documentation, and addressing any existing issues to make it robust and accessible.Posterior distributions derived from the retrieval analysis of simulated JWST spectra for HD 189733 b, where the spectra were generated using the full chemical network from Venot et al. (2020a).Retrievals were performed using the full chemical network, the reduced network, and the equilibrium chemistry as described by (Agúndez et al. 2012).Parameters guiding these simulations and retrievals can be found in Table 1.In the displayed posteriors, the blue and red curves represent the results from retrievals using the full and reduced networks with FRECKLL, respectively, while the green curve showcases results from the equilibrium chemistry-based retrieval.The light blue line indicates the true values as referenced in Table 1.
Values on the top of the posterior are from the full network retrieval.

Figure 3 .
Figure 3. Vertical abundance profiles for the main constituents of HD 189733 b computed with FRECKLL using the Full Venot2020 network (top) and the Reduced Venot2020 network (bottom).Dotted lines are equilibrium molecular profiles.Planetary and star parameters used are from Table1.

Figure 6 .
Figure 6.Simulated JWST observations of HD 189733 b (blue) without photodissociation with retrieval best-fit models (orange) for the reduced network.

Figure 7 .
Figure 7. Simulated JWST observations of HD 189733 b (blue) with retrieval best-fit models (orange) for the full network (top panel), reduced network (middle panel) and equilibrium (bottom panel).

B.
Figure B.1 shows the posterior distributions for the simulated JWST spectra retrieved with the reduced chemical network.

Figure
Figure D.1.Posterior distributions derived from the retrieval analysis of simulated JWST spectra for HD 189733 b, where the spectra were generated using the full chemical network fromVenot et al. (2020a).Retrievals were performed using the full chemical network, the reduced network, and the equilibrium chemistry as described by(Agúndez et al. 2012).Parameters guiding these simulations and retrievals can be found in Table1.In the displayed posteriors, the blue and red curves represent the results from retrievals using the full and reduced networks with FRECKLL, respectively, while the green curve showcases results from the equilibrium chemistry-based retrieval.The light blue line indicates the true values as referenced in Table1.Values on the top of the posterior are from the full network retrieval.

Table 1 .
Planetary and parent star parameters from Addison et al. (2019) used to generate the simulated HD 189733 b JWST transit spectra.
* we've elected to move the star further away to prevent saturation of the JWST NIRISS instrument.Kzz includes the log10 value in brackets.

Table 2 .
Retrieval priors and ranges.The log prefix describes fitting the ranges in log-space Table A.1.The 55 photodissociation reactions, their associated cross-sections and quantum yields included in the Venot et al. (2020a) chemical network.