Lightning: An X-ray to Submillimeter Galaxy SED-Fitting Code With Physically-Motivated Stellar, Dust, and AGN Models

We present an updated version of Lightning, a galaxy spectral energy distribution (SED) fitting code that can model X-ray to submillimeter observations. The models in Lightning include the options to contain contributions from stellar populations, dust attenuation and emission, and active galactic nuclei (AGN). X-ray emission, when utilized, can be modeled as originating from stellar compact binary populations with the option to include emission from AGN. We have also included a variety of algorithms to fit the models to observations and sample parameter posteriors; these include an adaptive Markov-Chain Monte-Carlo (MCMC), affine-invariant MCMC, and Levenberg-Marquardt gradient decent (MPFIT) algorithms. To demonstrate some of the capabilities of Lightning, we present several examples using a variety of observational data. These examples include (1) deriving the spatially resolved stellar properties of the nearby galaxy M81, (2) demonstrating how X-ray emission can provide constrains on the properties of the supermassive black hole of a distant AGN, (3) exploring how to rectify the attenuation effects of inclination on the derived the star formation rate of the edge-on galaxy NGC 4631, (4) comparing the performance of Lightning to similar Bayesian SED fitting codes when deriving physical properties of the star-forming galaxy NGC 628, and (5) comparing the derived X-ray and UV-to-IR AGN properties from Lightning and CIGALE for a distant AGN. Lightning is an open-source application developed in the Interactive Data Language (IDL) and is available at https://github.com/rafaeleufrasio/lightning.


INTRODUCTION
The light emitted from a galaxy contains a plethora of information about many physical properties of the system, ranging from its star-formation history (SFH) and dust content to the presence of an active galactic nucleus (AGN) and the properties of its supermassive black hole (SMBH). These properties are key to our current understanding of how galaxies and SMBHs formed and evolved, and, thus, the methods for deriving them from spectral energy distributions (SEDs) have been the Corresponding author: Keith Doore kjdoore@uark.edu focus of substantial work (e.g., Silva et al. 1998;Devriendt et al. 1999;Dale et al. 2005;Groves et al. 2008;Noll et al. 2009;Ciesla et al. 2015;Iyer & Gawiser 2017;Leja et al. 2018;Shanks et al. 2021). The overarching process of deriving the physical properties from an SED is known as SED fitting (see Walcher et al. 2011 andConroy 2013 for recent reviews). At its core, this process consists of fitting a model (e.g, stellar population synthesis with dust attenuation) to the observed SED. Once a best-fit model is determined using the chosen statistical inferencing method, it can be used to infer the physical properties of the observations. Numerous SED fitting codes currently exist today for the modeling and inferencing of galaxy properties from their SEDs. Some of the more widely cited codes include CIGALE (Burgarella et al. 2005;Boquien et al. 2019), Prospector (Johnson et al. 2021), MAGPHYS (da Cunha et al. 2008), BAGPIPES (Carnall et al. 2018), and FAST (Kriek et al. 2009); also see Pacifici et al. (2022) for a more comprehensive list. Each code was designed to help solve unanswered problems unique to their developers. Therefore, each code is unique and comes with its own set of advantages and disadvantages.
Initially, SED fitting codes were developed to use maximum-likelihood statistical inferencing methods (e.g., linear and non-linear optimization) to estimate galactic properties from optical to infrared (IR) observations (e.g., SEDfit, Sawicki & Yee 1998;Sawicki 2012;STARLIGHT, Cid Fernandes et al. 2005; VESPA, Tojeiro et al. 2007). These codes typically model the observations using simple stellar population (SSP) models (e.g., Fioc & Rocca-Volmerange 1997;Bruzual & Charlot 2003;Conroy et al. 2009;Eldridge et al. 2017) with simple parametric SFHs (e.g., exponentially declining) and attenuation. The advantage of these codes is that they are fast, simple to use, and return reliable bestfit models. However, the major drawback is that they can have difficulties computing accurate uncertainties on the physical parameters, since these parameters can be highly correlated and usually have non-Gaussian likelihoods. This difficulty is compounded as additional model components are included to account for more complex physical processes within galaxies: for example, non-parametric SFHs (see Carnall et al. 2019 andLeja et al. 2019 for overviews on the differences between parametric and non-parametric SFHs), dust emission (e.g., Draine & Li 2007;da Cunha et al. 2008;Casey 2012;Dale et al. 2014), dusty torus emission from an AGN (e.g., Fritz et al. 2006;Nenkova et al. 2008;Stalevski et al. 2012Stalevski et al. , 2016, and nebular emission (e.g., Ferland et al. 1998Ferland et al. , 2013. To estimate more accurate uncertainties, the next generation of SED fitting codes were developed to use a gridded Bayesian statistical inferencing method (e.g., CIGALE, Burgarella et al. 2005;Boquien et al. 2019;FAST, Kriek et al. 2009;MAGPHYS, da Cunha et al. 2008). This method estimates galactic properties and their uncertainties by gridding parameter space, fitting the corresponding gridded models to the observations, and then weighting the models by their goodness of fit. The advantage of this method is that it can account for parameter degeneracies and non-Gaussian likelihoods, while still being computationally fast. However, this computational speed excludes the time to create the grid of models, which increases exponentially with the number of parameters. Therefore, sampling of the full posterior distribution becomes intractable in a reasonable amount of time for complex models with many parameters unless parameter space is sparsely sampled.
In order to better sample the parameter space of complex models, new SED-fitting codes were developed to use Bayesian sampling statistical inferencing methods which utilize Markov Chain Monte Carlo (MCMC) and/or nested sampling algorithms (Skilling 2004). This approach (e.g., BAGPIPES, Carnall et al. 2018;BayeSED, Han & Han 2012BEAGLE, Chevallard & Charlot 2016;Prospector, Johnson et al. 2021) has the advantage of efficiently sampling parameter space of complex models to generate posterior distributions of parameters, while taking into account any prior information on the parameters. Additionally, models can be changed without any computational cost unlike the gridded Bayesian methods, which require the entire grid of models to be recomputed. However, the disadvantage of the Bayesian sampling approach is that sampling the posterior distribution can take significantly longer computational times on a per-SED basis. Some next generation SED fitting codes are trying to bridge the gap between parameter estimation and computational speed using machine learning. For these codes (e.g., Lovell et al. 2019;mirkwood, Gilda et al. 2021), machine-learning models are trained to learn the relationship between input observations and inferred properties using synthetic galaxy SEDs generated by cosmological simulations. The major advantage of these codes is that, once trained, fitting a new input SED is incredibly fast and derived parameters can be more accurate than the fully Bayesian approach (Gilda et al. 2021). However, in their current state of development, machine-learning SED-fitting codes come with a few serious drawbacks. The first is that they have an overreliance on theoretical models to explain how real galaxies should appear. Since they typically utilize a variety of cosmological simulations, it can be difficult to create a complete training set that is truly representative of all observed galaxies (Genel et al. 2014;Schaye et al. 2015;Somerville & Davé 2015). Additionally, overtraining (which can lead to overfitting) can occur when an appropriate test set or cross-validation set is not utilized. This can lead to galaxy property estimates that have high precision, which is a direct result of overfitting rather than a correct uncertainty estimate. Finally, these codes cannot handle missing observations that are commonly present in typical SEDs without retraining. This comes at a significant computational cost if the sample that is to be fit has a variety of observations. Motivated to have a computationally fast yet fully Bayesian code, we developed the SED-fitting code Lightning. Originally developed to derive the SFHs needed to empirically calibrate the X-ray luminosity function (XLF) of X-ray binary (XRB) populations Lehmer et al. 2017Lehmer et al. , 2020Lehmer et al. , 2022Gilbertson et al. 2022), Lightning has since been utilized to check for enhanced star formation and AGN activity in protoclusters (Monson et al. 2021), model local analogs to high-redshift galaxies (Motiño Flores et al. 2021), investigate the inclination-dependence of derived SFHs (Doore et al. 2021, and provide evidence in favor of density wave theory (Abdeen et al. 2022). Written in IDL, the newest updates to Lightning now allow for modeling of photometric SEDs from the Xrays to submillimeter using efficient MCMC algorithms to fit physical models that account for any combination of stellar, dust, and AGN emission. In Section 2, we describe these physical models and their dependencies. The statistical inferencing methods that fit these models to input SEDs are described in Section 3. In Section 4, we demonstrate the capability of Lightning by applying it to a variety of examples. Finally, we summarize and discuss future planned additions for Lightning in Section 5.
Lightning is an open-source, well-documented, and publicly available SED fitting code available at https:// github.com/rafaeleufrasio/lightning. Since Lightning has been in development over the past several years, we include the references to past works that first described each feature (i.e., models and statistical inferencing methods) and the motivation for their implementation in Table 1. In Sections 2 and 3, we reiterate the details from these past works and clarify all assumptions for replicability.

PHYSICAL MODELS
In this section, we describe the variety of physical models available in Lightning to account for any combination of stellar, dust, and AGN emission. To help clarify the free parameters corresponding to each model, we give a description of the parameters, their units, and their allowed range in Table 2.

Simple Stellar Populations
The intrinsic UV-to-IR stellar emission in Lightning is generated using the SSPs from the population synthesis code PÉGASE (Fioc & Rocca-Volmerange 1997). The SSPs, which are generated assuming the Kroupa (2001) initial mass function (IMF), are instantaneous bursts of star formation normalized to a unit star formation rate (SFR) of 1 M yr −1 . We allow for a wide range of   Monson et al. (2023) metallicity options when generating the stellar populations (0.001, 0.004, 0.008, 0.01, 0.02, 0.05, and 0.1 in terms of Z) along with the option to include the nebular extinction and emission from PÉGASE (see Section 2.4 of Fioc & Rocca-Volmerange 1997). While the nebular extinction affects SEDs of all ages, nebular emission is only added to stellar populations with ages <50 Myr. We make this simplifying assumption since there is minimal ionizing flux, which causes the nebular emission, for populations with ages >50 Myr (Smith et al. 2002;Byler et al. 2017). In Figure 1, we show some example SSPs used in Lightning at different ages for a metallicity of Z = 0.02 (i.e., solar metallicity).

Star Formation History
To model complex SFHs while remaining computationally fast, Lightning assumes a simple nonparametric SFH. Continuing with the original description in Eufrasio et al. (2017), we define this to be a piece-wise constant SFH, where the free parameters for the SFH are the SFRs (ψ j ) within the user-defined age bins. This is given in analytical form as a function of stellar age t by where t j and t j+1 are the respective lower and upper age bin boundaries of the jth bin. The advantage of normalizing by SFR, versus stellar mass like other SED fitting  Calzetti et al. (2000) Atten.  codes, is that any bias toward rising SFHs is prevented, while still allowing for bursty SFHs .
To compute the intrinsic, rest-frame composite stellar spectrum for the jth age bin,L ν,j (ν) 1 , the SSPs are integrated over the age bin after interpolating the SSP evolution to a common time grid using a user-defined time step (e.g., 0.5 Myr). The total composite stellar spectrum for all ages,L ν (ν), is then given bỹ whereL ν,j (ν) is by construction normalized per unit SFR, specifically 1 M yr −1 .
In Figure 2, we show composite stellar spectra for a metallicity of Z = 0.02 using the default set of age bins in Lightning (0-10 Myr, 10-100 Myr, 0.1-1 Gyr, 1-5 Gyr, 5-13.6 Gyr). These age bins were chosen such that 1 When symbolizing intrinsic emission (i.e., no attenuation), we include a tilde over the variable to clarify that it is intrinsic. Emission variables without a tilde signify attenuation/absorption has been applied.
the youngest age bin encapsulates the stellar population able to emit the majority of the ionizing flux, while the second bin includes the stellar population which generates the remaining bulk of the UV emission. Finally, the last three age bins were selected to have similar bolometric luminosities as the second bin in the case of a constant SFH.

X-ray Binary Model
We include stellar X-ray emission from compact object binaries in Lightning with a power-law spectral model given byL where we assume a photon index of Γ = 1.8 and cutoff energy of E cut = 100 keV. We set the normalization of the power-law by its rest-frame 2-10 keV luminositỹ L X , which we model to include contributions from both high-mass XRB (HMXB) emission and low-mass XRB (LMXB) emission: The PÉGASE SSP SEDs used in Lightning at various ages after ZAMS for a metallicity of Z = 0.02 and initial mass of 1 M . For an age of 1 Myr, nebular emission lines can be clearly seen, while the older displayed ages do not have any lines due to our simplifying assumption. Additionally, as the population ages, the overall bolometric luminosity decreases, and the peak wavelength of the emission shifts from the UV into the NIR.
Our models ofL HMXB X andL LMXB X are calculated using the empirical parameterizations with stellar age from Gilbertson et al. (2022): (t) = −0.24(log 10 (t) − 5.23) 2 + 32.54, andL LMXB X M (t) = −1.21(log 10 (t) − 9.32) 2 + 29.09, whereL X /M has units of erg s −1 M −1 and t is the stellar age in yr. While studies have shown that the luminosity of HMXBs depends on metallicity (e.g., Lehmer et al. 2021), we do not currently implement any metallicity dependence in our X-ray binary model. The agedependent relationship from Gilbertson et al. (2022) was derived for galaxies with metallicities ranging from 0.40-1.16 Z . For larger metallicities,L X /M may thus be slightly overestimated for t 100 Myr, while for lower metallicites it may be similarly underestimated.
To calculate the scaling parameterL X , we first derive the HMXB and LMXB contribution from each age bin utilizing Equations 5 and 6. We calculate the stellar mass of each bin (M ,j ) as where M coeff ,j is the coefficient that converts SFR into stellar mass. We then calculateL X /M at the mean age bin 1: 0 -10 Myr age bin 2: 10 -100 Myr age bin 3: 0.1 -1 Gyr age bin 4: 1 -5 Gyr age bin 5: 5 -13.47 Gyr Example composite stellar spectrum for a metallicity of Z = 0.02 using the default set of age bins in Lightning. Similarly to the SSPs in Figure 1, the nebular emission lines can be clearly seen in the youngest bin, while the older bins lack any obvious emission features. stellar age of each SFH bin, multiply by the stellar mass in the bin, and sum the contributions from each bin to derive the total contributions from the HMXBs and LMXBs, which are then incorporated into Equation 4 to determineL X and finally the X-ray luminosity spectrum.

AGN UV-to-IR Models
Lightning uses the SKIRTOR library of UV-to-IR AGN SED templates (Stalevski et al. 2012(Stalevski et al. , 2016, which consist of a broken power law accretion disk component and a clumpy two-phase dusty torus that reprocesses light from the accretion disk into the NIR. Our implementation uses the default accretion disk power law from the SKIRTOR library, where The full SKIRTOR templates have six free parameters: τ 9.7 , the edge-on optical depth of the torus at 9.7 µm; p, the power law index for the radial dust density gradient; q, the power law index for the polar dust density gradient; ∆, the opening angle of the dusty cone of the torus; R, the ratio of the torus' inner and outer radii; and i AGN , the inclination angle from the pole to the line of sight. To simplify the SKIRTOR models and allow Examples of the SKIRTOR AGN emission model generated by Lightning for a range of inclinations with τ9.7 = 7 normalized by LAGN. The dashed and dasheddotted lines for the iAGN = 90 • example are models where τ9.7 = 3 and τ9.7 = 11, respectively. We only show the variation in τ9.7 for the edge-on example, since τ9.7 is the edge-on optical depth of the torus. Therefore, changing its value for face-on to moderately inclined viewing angles minimally effects the model. us to sample the parameter space and interpolate between models, our implementation in Lightning only allows for τ 9.7 and cos i AGN to be free parameters. We fix p = 1 and q = 0 (i.e., there is no polar dependence of the dust density) as in Stalevski et al. (2016), and fix ∆ = 40 • based on their findings of typical covering factors in the range of 0.6 − 0.7. At the moderate covering factor of sin 40 • ≈ 0.64 that we assume, R has only a small impact on the luminosity of the torus, so we fix R = 20. To implement these simplified models, we linearly interpolate the original gridded models both in τ 9.7 and cos i AGN -space for the desired free parameter value. Then if no X-ray AGN model is used, the UV-to-IR AGN model is scaled using its anisotropic integrated luminosity, L AGN , as another free parameter. Examples of these UV-to-IR AGN emission models for different inclinations are shown in Figure 3.
We note that the UV-optical light that escapes the torus is subject to being attenuated by the host galaxy dust (see Section 2.3). When energy balance is enabled (as discussed in Section 2.4), we integrate this attenu-ated light over all lines of sight 2 and add it to the bolometric luminosity of the attenuated stellar light used to scale the dust model. However, fully energetic selfconsistency is not maintained by the AGN model. The ISM is assumed to be opaque to the ionizing Lymancontinuum radiation from the AGN, and the ionizing flux from the AGN does not currently contribute to the nebular emission component, which is built into our PÉGASE models.

AGN X-ray Models
X-ray observations can place powerful constraints on the bolometric luminosity of AGN, and as such they are very useful in AGN SED modeling. In Lightning, we provide two different models to generate the X-ray emission from the AGN component.
Since X-ray spectra of AGN are often empirically described by power-law models, we provide this as a basic option in Lightning. The equation for the power-law is given in Equation 3, where we use a fixed photon index of Γ = 1.8 and cutoff energy of E cut = 300 keV. To connect the X-ray model to the UV-optical AGN component, we use the Lusso & Risaliti (2017)L 2 keV −L 2500 relationship, which is an empirically calibrated relationship between the intrinsic monochromatic luminosities of luminous AGN at 2 keV and 2500Å. Therefore, the power-law X-ray and UV-to-IR AGN models are both scaled using L AGN as the only free parameter. We note that while Lusso & Risaliti (2017) show that their relationship has a dispersion of ∼0.21 dex, we do not currently provide any increased flexibility by allowing for a deviation from the best-fitL 2 keV −L 2500 relation. Instead, one of our general philosophies in Lightning is to provide model flexibility with physical parameters wherever possible.
To account for the scatter in theL 2 keV −L 2500 relationship in a physically-motivated way, we provide an implementation of the qsosed X-ray AGN models from Kubota & Done (2018), which reproduces the Lusso & Risaliti (2017) relationship, including its scatter, as a function of SMBH mass and Eddington ratio. These models include an accretion disk component and two Comptonizing regions, which produce the X-ray spectrum. Since these models include optical emission from the accretion disk, the relationship betweenL 2500 and L 2 keV is encoded in two model parameters: M SMBH , the  SMBH mass, andṁ, the Eddington ratio of the AGN. Kubota & Done (2018) show that their models reproduce the Lusso & Risaliti (2017) relationship, explaining the dispersion around the relationship due to variance of M SMBH andṁ among AGNs. Thus, when this model is used, we normalize the UV-to-IR AGN model (see Section 2.2.1) to the sameL 2500 as the X-ray model. This allows the free parameters of M SMBH andṁ to set the normalization of the entire X-ray-to-IR AGN model and encapsulate the variation in theL 2 keV −L 2500 relationship, thereby replacing L AGN as a free parameter. Due to its physically-motivated nature, we set the qsosed model to the default X-ray AGN model in Lightning and recommend its usage. However, we note that this model is most appropriate for luminous, high-accretion rate systems (log 10ṁ ranges from −1.5 to 0.3) and is not appropriate for low-luminosity AGN and Comptonthick AGN, the latter of which require more careful and complicated X-ray modeling with reflection components. We show some examples of the connected qsosed and SKIRTOR UV-to-IR AGN models in Figure 4. While the qsosed model extends to optical wavelengths and is used to normalize the UV-to-IR SKIRTOR model, we limit it to λ < 2 nm rather than directly joining the two models across the extreme-UV wavelength range. Therefore, the light gray line segments in Figure 4 show a linear interpolation between the two models for visualization purposes.
This implementation of the connection between the X-ray, optical, and IR AGN emission is a half-step to a fully energetically self-consistent model, in which the X-ray and IR AGN spectra are generated from the same assumed torus model and accretion disk spectrum. For steps toward connecting X-ray and IR spectral models using the same torus, see e.g. Tanimoto et al. (2020).

Dust Attenuation Models
To account for the variety of attenuation laws between and within galaxies, we include several prescriptions for attenuation in Lightning. For the UV-to-IR, these include the original and a variable form of the Calzetti et al. (2000) attenuation curve and the inclinationdependent attenuation curve from Doore et al. (2021). For the X-rays, these include the tbabs absorption model (Wilms et al. 2000) and the Sherpa atten model from Rumph et al. (1994). 3 2.3.1. Calzetti et al. (2000) Attenuation We implement the commonly used Calzetti et al. (2000) attenuation curve as the default in Lightning. The general details and format of this attenuation curve as implemented in Lightning are presented in Section 3.2 of Eufrasio et al. (2017). To briefly summarize these details, we used the Calzetti et al. (2000) curve as linearly extrapolated by Noll et al. (2009) at 1200Å to extend to the Lyman limit (912Å). To allow for more flexibility, we include the optional variable slope and 2175Å bump feature modifications as presented in Noll et al. (2009). Reformatting to use the optical depth, rather than attenuation in magnitudes, this variable diffuse dust attenuation curve is defined as where τ DIFF is the optical depth of the diffuse dust at wavelength λ, τ DIFF,V is the V -band (0.55 µm) normalization, k (λ) is the Calzetti et al. (2000) attenuation curve, D(λ) is the functional Drude profile parameterizing the 2175Å bump feature, and δ is the parameter controlling the variable slope. The Drude profile is defined as where we assume a bump strength of E b = 0.85 − 1.9δ following the results of Kriek & Conroy (2013), a bump FWHM of ∆λ = 350Å, and a central bump wavelength of λ 0 = 2175Å. Finally, we also include an optional birth-cloud attenuation component given by which is only applied to the youngest defined SFH age bin (ψ 1 ). 4 Combining the diffuse dust and birth-cloud components, the effective optical depth of the attenuation curve for a given age bin j is given by where δ j1 is the Kronecker delta (not to be confused with the variable slope parameter δ). Therefore, the variable Calzetti et al. (2000) attenuation has up to three free parameters (τ DIFF,V , δ, and τ BC,V ) that define the shape of the curve.
In Figure 5, we show several modified Calzetti et al. (2000) attenuation curves as colored lines to compare with the base Calzetti et al. (2000) attenuation curve (i.e., δ = 0, τ BC,V = 0, and no UV bump feature), which is shown in black. The variation in color corresponds to a change in δ (left panel), while different line styles correspond to variations in τ BC,V (right panel). This comparison clearly shows how the variable slope parameter δ affects the slope of UV attenuation. Additionally, the inclusion of the birth cloud attenuation increases the amount of attenuation at all wavelengths.

Inclination-dependent Attenuation
To allow for more accurate attenuation in disk galaxies, an inclination-dependent attenuation curve is also included in Lightning, the description of which is presented in detail in Section 4.3 of Doore et al. (2021). To give a brief description, the prescription is based on the Tuffs et al. (2004) inclination-dependent attenuation curves, as updated by Popescu et al. (2011), which assume that disk galaxies are comprised of three components, a young thin disk, an old thick disk, and an old dustless bulge. The equation defining the curves in Tuffs et al. (2004) was restructured to depend on the intrinsic properties of these three components rather than their observable properties. This resulted in the attenuation curve being defined by ∆m λ = −2.5 log r 0,old 1 + B/D 10 Here, ∆m λ is the attenuation in magnitudes at a given wavelength λ. r 0,old is the fraction of intrinsic flux density from the old stellar components (i.e., thick disk and bulge) compared to the total intrinsic flux density. B/D is the intrinsic bulge-to-thick disk ratio. ∆m disk are the attenuations from the diffuse dust, parameterized by fifth-order polynomials in Popescu et al. (2011), which are functions of inclination, i, for tabulated values of the B-band face-on optical depth as seen through the center of the galaxy, τ f B , and wavelength for the disk, thin disk, and bulge, respectively. Finally, F is the birth cloud clumpiness factor of the thin disk, and f λ is a tabulated function of wavelength that provides F its wavelength dependence.
The restructuring of the original Tuffs et al. (2004) equation to intrinsic properties was intentional so that the non-parametric SFH in Lightning could be used to Calzetti et al. (2000) attenuation curves generated by Lightning, normalized by τDIFF,V . The black curve gives the base Calzetti et al. (2000) attenuation curve (i.e., δ = 0, τBC,V = 0, and no UV bump feature). The variation in color in the left panel corresponds to a change in δ with fixed τBC,V = 0, while different line styles in the right panel correspond to variations in τBC,V with fixed δ = 0. effectively eliminate r 0,old as a free parameter. This was done by making r 0,old a binary parameter, where each SFH age bin is given its own value for r 0,old based on its age. A value of 0 indicates that the given SFH bin is to be considered part of the young stellar population (e.g., t < 500 Myr), while a value of 1 considers it part of the old stellar population (e.g., t > 500 Myr). Therefore, with r 0,old tied to the SFH ages, the other four parameters (i, τ f B , B/D, and F ) are the free parameters defining the shape of the inclination-dependent attenuation curve. Examples of these inclination-dependent attenuation curves can be found in Figures 7 and 8 of Doore et al. (2021).

X-ray Absorption
In Lightning, the X-ray absorption is modeled using one of two user-selected X-ray absorption models: the tbabs model with the default Wilms et al. (2000) abundances or the Sherpa atten model which combines cross-sections from Morrison & McCammon (1983) and Rumph et al. (1994). The tabulated curves used in Lightning were generated with Sherpa 5 v4.13, and normalized to a line-of-sight HI column density of N H = 10 20 cm −2 . At energies larger than 10 keV, these models produce negligible absorption; for convenience we set the optical depth to 0 at > 12 keV.
The chosen X-ray absorption is first applied to the SED model in the observed frame to account for Galac-5 https://cxc.cfa.harvard.edu/sherpa/ tic absorption by the Milky Way, with the Galactic N H being a user input for each galaxy. Further intrinsic absorption is then applied in the rest frame on the stellar binary population and the AGN emission model, if applicable. If there is no X-ray AGN emission, the N H of the stellar population becomes a free parameter in the model. We show some examples of the X-ray stellar model with varying levels of absorption and SFHs in Figure 6. The inclusion of X-ray absorption primarily impacts lower energy photons as higher energy X-rays are less likely to be absorbed, with the intensity of the absorption increasing with N H .
However, when the model includes an X-ray AGN component, absorption of the stellar X-ray emission is not a free parameter, and is instead linked to the V -band attenuation via This scaling was chosen to be the average of observed Milky Way N H − A V relationships (Predehl & Schmitt 1995;Nowak et al. 2012). In this case, the N H value of the nuclear region becomes the free parameter in the model, as we expect the X-ray emission from the AGN . Example X-ray stellar spectra with and without X-ray absorption generated by Lightning. The used absorption model is the tbabs model with Wilms et al. (2000) abundances. The variation in color in the left panel corresponds to a change in the SFH parameters ψj with fixed NH = 1 (×10 20 cm −2 ). (The age bins for these ψj values are the default set of age bins in Lightning as discussed in Section 2.1.2.) The variation in lines styles in the right panel correspond to changes in NH with fixed ψj = [1, 1, 1, 1, 1] M yr −1 . The blue line in the left panel shows an example of a spectrum dominated by young stars (i.e. HMXBs), while the red line shows a spectrum dominated by old stars (i.e., LMXBs).
to be the dominant component of the X-ray spectrum in most cases. 6 We note here that while N H is allowed to increase above 10 24 cm −2 in our implementation, our X-ray emission models are not currently appropriate for Comptonthick sources, which typically require reflection components that are not included in Lightning.

Dust Emission Models
We model IR dust emission in Lightning using the Draine & Li (2007) model. To briefly summarize, the model details how the dust mass, M dust , is exposed to a radiation field intensity, U . Analytically, this is given by Equation 23 in Draine & Li (2007), where the first additive component is a delta function at the minimum radiation field intensity U min and the second component is a power-law of slope α derived between U min and U max , the maximum radiation field intensity. The parameter γ in each additive component dictates the fraction of the dust mass exposed to the power law compared to the delta function. Additionally, the polycyclic aromatic hydrocarbon (PAH) index, q PAH , is a hidden parameter in the model and defines the strength of the PAH emission.
Rather than scaling the model with M dust , we instead use the total integrated IR (TIR) luminosity L TIR (i.e., the bolometric luminosity of the dust model), which is proportional to M dust . The reason we make this substitution is because of the mechanism generating the dust emission. Simply put, some fraction of UV-to-NIR emission in a galaxy is attenuated by dust. This attenuated energy is conserved via radiation from dust particles at longer wavelengths, mainly across the mid-to-far IR. To account for this conservation of energy, it is expected that the bolometric luminosity of the attenuated light should be equal to the TIR luminosity. This energy conservation (often termed "energy balance" in the SED fitting community) is optional when fitting with the dust emission model in Lightning. Therefore, when fitting with energy balance, the dust model has five free parameters: α, U min , U max , γ, and q PAH ; and when energy balance is not assumed, L TIR becomes an additional free parameter used for normalization. In Examples of the Draine & Li (2007) dust emission model generated by Lightning, normalized by LTIR. The parameters Umax and α are fixed to 3 × 10 5 and 2, respectively. The variation in color in the left panel corresponds to a change in Umin with fixed γ = 0.01 and qPAH = 0.02. The different lines styles in the middle panel correspond to variations in γ with fixed Umin = 1 and qPAH = 0.02. As for the right panel, the different symbols on the dotted lines correspond to different values of qPAH with fixed Umin = 1 and γ = 0.01. Decreasing the value of Umin can be seen to shift the peak of the emission to shorter wavelengths, while qPAH independently increases the intensity of the PAH emission features.

put parameters. In these examples and as a default in
Lightning, we fix U max = 3 × 10 5 and α = 2. We make this simplifying assumption as recommended by , since they found that the dust model is not very sensitive to these two parameters and most observed IR SEDs are well reproduced by U max = 10 6 and α = 2. We note the discrepancy between the fixed value of U min in Lightning and that recommended by . The reason for the discrepancy comes from how Lightning computes the dust emission model from the publicly available data. To allow for a variable α, the δ-functions of U must be used. However, the largest available δ-function is U = 3 × 10 5 . Therefore, U max is limited to this largest available value, since extrapolating to U = 10 6 would add unwanted uncertainty to the model.

Observational Information
To keep Lightning computationally fast, we restrict it to modeling only photometric observations, since spectroscopic observations would require additional modeling assumptions. However, any combination of narrowto broadband flux densities that have been corrected for Galactic extinction can be used to make up the UV-tosubmillimeter input SED. Additionally, any number of uniform sensitivity top-hat energy bands can be used for X-ray observations, whose measurements can be in the form of either net counts or fluxes.
Since the models in Lightning are in rest-frame luminosity units, an observed distance indicator is required to convert the SED fluxes to the same luminosity units as the model. The distance indicator can simply be a luminosity distance, or it can be a redshift, where the luminosity distance is calculated from the redshift via the user chosen cosmology. We note that when using a redshift that the assumed age of the Universe, t age (z), will decrease as redshift increases. To account for any effect this will have on the SFH age bins, we have designed Lightning to automatically adjust the user-defined age bins such that upper bin boundaries of t j+1 > t age (z) will be updated to t j+1 = t age (z), and bins with lower bin boundaries of t j > t age (z) are completely omitted from the SFH.

Loss Function
In order to fit a given model to any data, a loss function, which determines how well the model fits the data, is required. In Lightning, we implement a χ 2 loss function given by where L is the likelihood probability of the model, L obs ν,f is the observed luminosity in filter f as derived from the input flux,L mod ν,f is the model photometry in filter f , and σ total,f is the total uncertainty associated with filter f . The model photometry is derived by integrating the observed-frame model spectrum through the corresponding filter f usinḡ where T f (λ) is the filter transmission function in units of energy, and L mod ν (λ) is the model spectrum. The total uncertainty consists of both the observed uncertainty, which is the Gaussian uncertainty of the measurement plus any additional fractional calibration uncertainty, and the model uncertainty added in quadrature or where σ obs,f and σ mod,f are the observed and model uncertainties, respectively. We include a model uncertainty component in Lightning to account for any oversimplification of models and potential systematic uncertainties in the models themselves (Charlot et al. 1996;Percival & Salaris 2009;Conroy et al. 2009. The model uncertainty is computed simply as a user-defined fraction of the model photometry that is constant for each filter or where σ frac mod is the user-defined fractional model uncertainty. Therefore, by including model uncertainties, our formulation of the χ 2 loss function accounts for both the observational uncertainty and the inherent uncertainty of the models being used.
When fitting X-ray data in terms of counts, we note that Equation 16 will not be applicable, since the Poissonian nature of the counts is not compatible with uncertainty formulation. Instead, we calculate the χ 2 contribution of the X-ray counts as where χ 2 X is X-ray counts contribution to the total χ 2 , N obs e is the observed net counts (i.e., background subtracted) in energy band e, N mod e is the model net counts in energy band e, and σ N,e is the approximate count uncertainty in energy band e. Since the approximate count uncertainty can be computed in a variety of ways depending on the overall count rate, we allow for the user to either input their own pre-computed count uncertainties or use one of the two built-in methods in Lightning. The first method simply sets the count uncertainty equal to the square root of the counts (i.e., σ N,e = N obs e ), since this is the Gaussian approximation in the high count regime. The other uses the upper uncertainty of the Gehrels (1986) approximation given by which is more appropriate for data in the low count regime. User-supplied uncertainties may be used when more flexibility is required, e.g., when the background contributes strongly to the uncertainty.

MPFIT Algorithm
To allow for the determination of a rapid best-fit solution to an SED without sampling the parameter posteriors, we have added a new maximum-likelihood inferencing method to Lightning. The method utilizes the MPFIT code (Moré 1978;Markwardt 2009), which consists of the gradient-descent Levenberg-Marquardt algorithm used to solve nonlinear least-squares problems. We chose the MPFIT implementation since it allows for several necessary constraints in Lightning, such as fixing parameters and setting parameter bounds.
When searching for the best solution, one drawback to gradient-decent algorithms like MPFIT is that they are prone to getting stuck in local minima in χ 2 space before reaching the global minimum. To mitigate this drawback, we have implemented the algorithm to run a user-defined number of "solvers", where each solver is a fit to the SED using different starting locations in parameter space (see Section 3.5). By running multiple solvers from different starting locations, we can compare the solver solutions. If the majority of the solvers have converged to the same solution, then we can be confident that this is likely the best solution and global minimum. More specifically, we require that (1) at least 50% of the solvers have χ 2 − 4 ≤ χ 2 best , where χ 2 best is the minimum χ 2 value of all solvers, and (2) all of these solvers satisfying criterion (1) have free parameter values that are within a 1% difference of the best fit. 7 Finally, once convergence to the best solution has been confirmed, the best-fit solver is considered the solution to the SED fit.

Bayesian Inferencing
To allow for a high quality sampling of the posterior to an SED model, we implement two Bayesian inferencing methods in Lightning. Both methods utilize MCMC algorithms to sample the posterior distributions of the free parameters, while incorporating prior distributions of these parameters.

Prior Distributions
In terms of analytical priors, we only include two basic options in Lightning: truncated uniform and normal (Gaussian) priors. We implement by default some truncated priors, since practically all possible free parameters have at least a lower bound. Therefore, limiting the 7 The chosen difference in χ 2 of 4 and the 1% difference are arbitrarily chosen from our test fits to well behaved SEDs. We include the χ 2 result and the parameter values for each solver in the output such that a user can further evaluate these cutoffs at their discretion.
prior range to comply with the physically allowed values is required. Additionally, both analytical priors are implemented in each parameter's sampled space (e.g., a parameter sampled in log space has a uniform or normal prior in log space). Besides the analytical priors, we include the option for a user to input a prior of any shape in tabulated form. The only restriction for the shape of these priors is that they conform to a parameter's physically allowed range. Otherwise, any shape is allowed, which can be useful for creating complex prior distributions for a given parameter. However, we do note that no free parameters can be linked together via the prior, tabulated or analytical. We exclude user-specified parameter linking in Lightning, which is common in other SED fitting codes, to minimize the computational complexity and increase computational speed. Additionally, we automatically link commonly correlated parameters (e.g., τ DIFF,V and N H ), which would potentially be linked by the user.
With the priors specified, the posterior probability is calculated using where P post is the posterior probability and P prior is the prior probability. 8 Since the goal of the MCMC algorithms is to sample P post , the loss function from Section 3.2 is updated for these algorithms to include the prior information such that we are minimizing and sampling − log(P post ) space rather than χ 2 space.

Adaptive MCMC Algorithm
The original MCMC algorithm adopted in Lightning was implemented and discussed in Doore et al. (2021). The algorithm is an adaptive version of the standard Metropolis-Hastings algorithm (Metropolis et al. 1953;Hastings 1970) created by Andrieu & Thoms (2008). The algorithm simply adjusts the proposal density distribution to achieve an optimal acceptance ratio. Additionally, this adjustment of the proposal density is vanishing, meaning the adaptiveness decreases with each subsequent iteration. Therefore, after many iterations the adaptiveness is insignificant and the algorithm is practically equivalent to the standard Metropolis-Hastings algorithm.
Similar to the MPFIT algorithm, a single chain of the adaptive MCMC algorithm is prone to getting stuck in local minima in − log(P post ) space. Once stuck, it can take more than the user-specified number of trials to escape and move towards the global minimum. 9 To confirm if a chain reached the vicinity of the global minimum, we have designed the adaptive MCMC algorithm to run a user-defined number of independent chains in parallel, where each chain is initialized using different locations in parameters space (see Section 3.5). By running parallel chains from different starting locations, we can compare the ending segment of each chain to see if they converged to a single best solution. To test for convergence, we have Lightning automatically perform the Gelman-Rubin test (Gelman & Rubin 1992;Brooks & Gelman 1998) and its multivariate version from Brooks & Gelman (1998) on the ending segment of the chains, whose length is user-defined. If the ending segments from the parallel chains results in acceptable Gelman-Rubin statistics (i.e., R ≤ 1.2), then a user can confidently concluded that convergence has been reached and the posterior well sampled by the algorithm.
Since the ending portion of the full parallel chains, which samples the posterior, is the main interest of a user, Lightning automatically post-processes the full chains to create a final post-processed chain portion. For the adaptive MCMC algorithm, this post-processed chain portion is determined as follows. First, each parallel chain has a user-defined number of initial trials discarded as the burn-in phase. Next, these truncated chains are thinned by a user-specified thinning factor (i.e., only every n elements of each chain is kept). Finally, the thinned and truncated chain containing the highest posterior probability is selected, and the ending segment of this highest probability chain, whose length is user-defined, is kept as the sampled posterior.
The reason for selecting the highest posterior probability chain is based on the assumption that convergence may not have been reached. If convergence was reached, then all parallel chains will have very similar distributions, and it does not matter which one is selected for use. However, if convergence was not reached, selecting the chain with the highest probability, guarantees that the chain with the best solution is selected.

Affine-Invariant MCMC Algorithm
To more quickly sample the potentially skewed posteriors of the free parameters, we have recently added an affine-invariant MCMC algorithm to Lightning (Monson et al. 2023), which is our default algorithm. This algorithm uses an ensemble of samplers to adjust the pro-posal density distribution and sample the posterior distribution. The ensemble consists of multiple chains (or walkers) that are run in parallel and allowed to interact with one another so that they can adapt their proposal densities. For our implementation in Lightning, we use the affine-invariant "stretch move" method as presented in Goodman & Weare (2010), which was shown to more quickly sample non-Gaussian and skewed posteriors compared to the Metropolis-Hastings MCMC algorithms.
Unlike the MPFIT and adaptive MCMC algorithms, the ensemble nature of the affine-invariant MCMC algorithm typically prevents it from getting stuck in local minima in − log(P post ) space, since each walker is initialized using different locations in parameters space (see Section 3.5). However, it is still important to confirm that the ensemble has converged to a stationary solution and to quantify any potential sampling error effects, since each walker in the ensemble is not independent. To test for this convergence, we have Lightning automatically perform an autocorrelation analysis on the ensemble (see Weare 2010 andForeman-Mackey et al. 2013 for detailed discussions on autocorrelation analyses). We briefly summarize the idea and methods of the analysis as applied in Lightning below.
Autocorrelation, in the context of MCMC algorithms, is how correlated a sample is with previous samples from the same walker or chain. Specifically in Lightning, we look at the integrated autocorrelation time, which is a measure of the average number of iterations between independent samples. If the autocorrelation time is large, then the samples in the ensemble are likely highly correlated and contain few independent samples. To confirm the ensemble has enough independent samples and to quantify the Monte Carlo error, we have designed Lightning to calculate the autocorrelation time and check for convergence following the methods of Foreman- Mackey et al. (2013). Their methods recommend that the MCMC algorithm runs for a number of iterations at least ∼50 times the integrated autocorrelation time in order for one to trust that the autocorrelation time estimate is accurate and convergence has been reached. A factor fewer than ∼50 can cause the autocorrelation time to be underestimated, which could result in a highly correlated sampling with few independent samples. Therefore, we have Lightning check if the user-defined number of iterations is large enough to have an accurate autocorrelation time estimate (i.e., autocorrelation time ≥50) and flag fits that do not.
Similar to the adaptive MCMC, the ensemble is automatically post-processed by Lightning to produce a final post-processed chain portion. For the affine-invariant MCMC algorithm, the post-processed chain portion is determined as follows. First, each walker in the ensemble has its burn-in phase discarded, where the burn-in phase is a number of iterations from the beginning of the chain equal to either twice the longest autocorrelation time of any parameter in the ensemble (double the autocorrelation time typically encapsulates the entire burn-in phase) or a user-defined value. Next, the truncated ensemble is then thinned by a thinning factor, which is either half the longest autocorrelation time of any parameter in the ensemble (iterations at half the autocorrelation time typically give fully independent samples) or user-defined. Then, if a walker in the thinned and truncated ensemble is classified as stranded (we explain how we classify walkers as stranded below), they are removed from the ensemble. Finally, the nonstranded ensemble is flattened element-wise into a single chain and the ending segment of this flattened chain, whose length is user-defined, is kept as the sampled posterior.
We designed Lightning to classify walkers as stranded if they have an acceptance fraction less than a userspecified number of standard deviations below the median ensemble acceptance fraction. Due to the boundaries of the free parameters, the affine-invariant MCMC can have trouble accepting moves of walkers separated from the ensemble (typically at higher − log(P post ) values) when the ensemble is near a boundary. This results in the walkers becoming stranded and having very low acceptance rates, since they are failing to have any proposal jumps accepted. With enough iterations, these walkers will eventually have a jump that rejoins them with the ensemble. However, only a finite number of iterations are allowed for this to occur. Therefore, once the specified iteration limit has been reached, any stranded walkers that may remain need to be removed, since they would add faulty samples to the final sampled posterior. We have found that the most effective automatic method for correctly selecting stranded walkers is to compare each walker's acceptance fraction with that of the median of the ensemble. Those that are classified as stranded with abnormally low acceptance fractions compared to the rest of the ensemble are consistently considered stranded when using more robust and manual visualization methods.

Algorithm Initialization
All three of the current algorithms in Lightning require initial starting values for each free parameter. To select these starting values, Lightning randomly samples the user-specified prior distribution of each parameter independently. (For the MPFIT algorithm, there are technically no priors, since it is not a Bayesian inferencing method. However, we assume uniform "priors" for all parameters for the purpose of initialization.) Additionally, since a given prior may have a much larger range than what would constitute an appropriate starting range, we allow the user to limit the initialization to a specified range within the prior range. Therefore, each initialization of the algorithms (i.e., each unique solver, chain, or walker) is initialized randomly from a potentially limited range of the corresponding prior.

Derived Quantities
After fitting an SED with the chosen algorithm, we designed Lightning to automatically do additional postprocessing to derived typical quantities of interest. In terms of physical properties that are not free parameters in the model, the total and individual model component spectra and photometry with and without attenuation are derived. When using an MCMC algorithm, Lightning allows for a user to select and keep the model spectra for a specified fraction of best-fit elements in the final post-processed chain portion. This is useful for deriving and quoting model uncertainties on new simulated photometric data points as well as showing model uncertainties in the spectra when plotting. Additional physical properties are also derived, such as the stellar mass and L TIR , if their corresponding model component is included in the total model. One other important quantity of interest that Lightning automatically derives is a p-value from a goodness-of-fit test for each SED fit. Goodness-of-fit tests are often overlooked in SED fitting, but they are important for determining if the chosen model can acceptably model the data and whether the uncertainties are trustworthy. For the MPFIT algorithm, we use a χ 2 goodness-of-fit test to derive the p-value using the χ 2 and degrees of freedom as given by the MPFIT algorithm. We caution against using this p-value to reject the null hypothesis that the chosen model can acceptably model the given SED. Since the effective number of free parameters is lower than the actual number (due to degeneracies and covariances between parameters), the number of degrees of freedom is likely higher than what is given by MPFIT. Therefore, this p-value can be underestimated and can lead one to falsely conclude that the model is not acceptable for the given SED.
As for the MCMC algorithms, since they are Bayesian methods, Lightning performs a posterior predictive check (PPC; Rubin 1984;Gelman et al. 1996) to derive a p-value for the chosen model. A PPC is a goodness-of-fit test that uses the model itself to estimate the distribution from which the p-value is derived. The model can be considered an accurate description of the data if it can generate simulated (replicated) data that is statistically identical to the actual data. By considering the replicated data as data that could have been measured, the PPC tests whether the model encapsulates all of the variability in the actual data.
In terms of the practical application in Lightning, a PPC takes the following steps. First, replicated sets of model photometry are randomly selected from the posterior distribution (i.e., samples are bootstrapped from the post-processed chain portion with the chance of selection being weighted by their posterior probability). Then, the replicated sets of photometry are randomly perturbed by a Gaussian distribution with a variance corresponding to the respective total uncertainty. Next, likelihood probabilities are calculated (see Section 3.2) by comparing the actual observations and the perturbed replicated photometry with the unperturbed replicated photometry. Finally, the p-value is computed as the fraction of corresponding likelihood probabilities for the perturbed replicated photometry that are smaller than the likelihood probabilities for the actual observations. 10

EXAMPLE APPLICATIONS
In this section, we present different example applications of Lightning to give interested users ideas of its capabilities. In Figure 8, we show composite postage stamp images of the galaxies used in each example, along with a brief description of the example topic. We note that, within the "examples" sub-directory of the Lightning GitHub repository and online documentation 11 , we provide the scripts required to run the examples and generate their presented figures. For interested users, we recommend following along in the online documentation in addition to the below text to get further supplementary details.

Property Maps of M81
To show the power and versatility of Lightning when applied to a nearby galaxy, we fit the spatially resolved UV-to-FIR photometry (SED map) of the nearby spiral galaxy, M81 (NGC 3031). To generate the SED map, we gathered 23 publicly-available photometric images ranging in wavelength from GALEX FUV to Herschel 350 µm. We then preprocessed these images for fitting by (1) subtracting foreground stars, (2) convolving each  image to a common 25 PSF, (3) re-binning them to a common 10 pixel scale, (4) estimating the background to update the photometric uncertainties, (5) correcting each bandpass for Galactic extinction, and (6) combining the images into a data cube which contains the pixelby-pixel SEDs. Since any data pre-processing steps are done external to running Lightning, we exclude the intricate details on our image pre-processing methods as the focus of this example is on the application of Lightning rather than the creation of SED maps. Interested readers can refer to Section 2 of Eufrasio et al. (2017) for a detailed description of the pre-processing steps involved.
To model the SEDs, we used a stellar population with solar metallicity (i.e, Z = 0.02) and SFH age bins of 0-10 Myr, 10-100 Myr, 0.1-1 Gyr, 1-5 Gyr, and 5-13.6 Gyr. The stellar emission was attenuated using the modified Calzetti et al. (2000) curve with the 2175Å bump feature and excluding any birth cloud attenuation. Finally, the dust attenuation was set to be in energy balance with the Draine & Li (2007) dust emission model, where we fixed U max = 3×10 5 and α = 2, as recommended by  and discussed in Section 2.4. To fit the model to each SED, we utilized the affine-invariant MCMC algorithm, which we ran for 10 4 trials with 75 walkers, assuming 5% model uncertainty. For all free parameters, we implemented uniform priors over either the entire available range or a broad range of values if the available range is infinite. The values defining the priors associated with each free parameter in the model, along with the limited initialization ranges, are listed in Table 3.
With the described model and algorithm, we used Lightning to fit a subset of the SEDs within the SED map, assuming a luminosity distance to M81 of 3.5 Mpc as given in Dale et al. (2017). The subset included all SEDs that were inside the 25 mag arcsec −2 B-band Table 3. Summary of Parameters Used in the M81 Example.

Parameter
Prior Function Initialization Range Note-U (a, b) indicates a uniform distribution from a to b. Fixed parameters have their value listed in the initialization range column.
isophotal ellipse as given by HyperLeda 12 (Makarov et al. 2014), which limited our example to the general extent of the optical emission of the galaxy. The fitting process took 65.1 hours total to fit all 6972 SEDs using one 32-core, 2.1 GHz CPU on the Arkansas High Performance Computing Center (i.e., 0.30 core-hours per SED), with Lightning automatically running each SED fit in parallel to maximize CPU usage. To give a general sense of fitting speed, other Bayesian sampling SED fitting codes typically take 1-100 core-hours per SED depending on the complexity of the chosen model, which means Lightning fits 1 order of magnitude faster than other codes. Once Lightning is finished fitting, it automatically post-processes the fit to each SED as described in Sections 3.4.3 and 3.6 and combines the results into a single file. When generating the final post-processed chain portion, we manually set the length of the burn-in phase From these maps, it can be seen that the attenuation and recent star formation is concentrated in the spiral arms, while the stellar mass is concentrated in the bulge. This age stratification is clear in the sSFR map, which uses a different color scheme to highlight younger populations in blue and older populations in red. and thinning factor for consistency rather than having Lightning automatically determine them for each SED from the autocorrelation times. We chose a burn-in length of 8000 trials and a thinning factor of 250, which is significantly larger than the overall maximum autocorrelation time of 178. Therefore, each element in the final chain portions should be uncorrelated. Finally, we only keep the final 250 elements of each chain so that we can derive reasonable median and 16th and 84th percentile ranges, while minimizing the total memory of the single post-processed file.
With the post-processed results, we first checked that all SEDs converge to stationary solutions as determined by their autocorrelation times derived from the full chain (see Section 3.4.3). After confirming convergence, we then mapped the derived quantities for each SED back to its associated pixel to generate maps of the spatially resolved properties of the galaxy. In Figure 9, we show a variety of these derived spatially resolved properties, with the values of each pixel being the median value from the corresponding posterior distribution. The only exception to this is the image of the p-values in panel (f ), which only has a single value rather than a distribution. From this image, it can be seen that the vast majority of pixels are well fit by the model, with poor fits mainly being associated with locations where foreground star subtraction occurred. One stand out result from these property maps is that the younger, star-forming population, which is more highly obscured, is concentrated in the spiral arms as seen in the SFR (SFR of the last 100 Myr) and V -band optical depth maps. Additionally, the older, more massive population can clearly been seen to Resulting SFH with uncertainty ranges for the outer and inner regions as the blue and orange lines, respectively. From the SEDs and SFHs, it can be seen that the outer region has comparatively higher UV emission and lower NIR emission, which is distinguished by the fit as an overall younger population compared to the inner region. reside in the bulge region, where the stellar mass is high and SFR is low, resulting in a low specific star formation rate (sSFR; sSFR = SFR/M ).
To further demonstrate how the spatially resolved results could be used to estimate properties for regions of the galaxy, we separate the galaxy into two parts, the outer and inner regions, with the outer region being defined as the pixels between the 25 mag arcsec −2 B-band isophotal ellipse as given by HyperLeda and one half of the 20 mag arcsec −2 K s -band isophotal ellipse as given by Jarrett et al. (2003), and the inner region comprised of pixels within one half of the 20 mag arcsec −2 K s -band isophotal ellipse. In the left of Figure 10, these regions are shown overlaid on the convolved and re-binned SDSS color image as the blue and orange ellipses, respectively. By summing the results of each pixel within each region, a total SED and SFH can be made as shown in the middle and right of Figure 10, respectively. These results show, as inferred from the maps in Figure 9, that the outer region has comparatively higher UV emission and lower NIR emission, which is distinguished by the fit as an overall younger population compared to the inner region.
Finally, to give a sense of how the maximum likelihood and Bayesian algorithms in Lightning compare, we refit the B-band isophote subset of the SED map after swapping the MCMC algorithm with the MPFIT algorithm. Using 20 solvers to test for convergence, fitting took 19.1 hours to fit all 6972 SEDs on one of the same 32-core, 2.1 GHz CPUs (i.e., 0.09 core-hours per SED), which is 3.4 times faster than the MCMC algorithm. In Figure 11, we show the same property maps as in Figure 9, except for the best fit as determined by MPFIT. As expected, there are variations in the results between algorithms, especially in the smoothness of the SFR and, subsequently, sSFR maps. This loss of smoothness is primarily due to the usage of the best fit rather than medians of the posterior, since the median of the posterior is rarely the same as the best-fit solution. For example, the central region can be seen to have best-fit solutions that are suggesting zero recent star formation (the colorbar is truncated due to the log-scale and to match the colorbar ranges in Figure 9). These zero values from MPFIT in low SFR regions are not unexpected, since SFR posterior distributions from the MCMC algorithm generally provide upper limits, with the maximum (f) Figure 11. Panels (a -e) are the same as Figure 9 except the values correspond to the best-fit values derived from the MPFIT algorithm. The p-values in panel (f ) are estimated from the a χ 2 goodness-of-fit test. The colorbar minimum and maximum ranges have been truncated to match those in Figure 9. Comparing these maps to those in Figure 9, it can be seen that some of these maps, particularly the SFR and sSFR, are not nearly as smooth, due to the use of best-fit values rather than the posterior median.
probability coinciding with zero. To show how this variation in best fit and median affects the total SFR of the galaxy, each pixel in the SFR map can be summed to give a total SFR. This results in SFR = 0.22 M yr −1 and SFR = 0.44 M yr −1 for the MPFIT and MCMC algorithms, respectively, a difference by a factor of two. In contrast, the stellar mass, in general, is better constrained by SED fits, and the best fit is typically closer in value to the median of the posterior. Summing the pixels for the masses gives M = 4.05 × 10 10 M and M = 4.26 × 10 10 M for the MPFIT and MCMC algorithms, respectively, a difference of only 5%. Thus, the best-fit mass map from MPFIT appears nearly identical to the median mass map from the MCMC.

Deep Field AGN
Here, we demonstrate how the AGN module in Lightning can be used with and without X-ray data by fitting the SED of J033226.49−274035.5, an X-ray detected AGN in the Chandra Deep Field South at z = 1.03. For both fits, the UV-to-NIR photometry were retrieved from the Guo et al. (2013) CAN-DELS catalog, which covers the U -band to Spitzer IRAC 8.0 µm. From this data, we excluded the VLT/VIMOS U and HST/ACS F435W bands from our SED, due to potential contamination by broad-line emission, since J033226.49−274035.5 is classified from optical spectra as a type 1 AGN (Szokoly et al. 2004). We additionally included the FIR data (Spitzer MIPS 24 µm to Herschel SPIRE 250 µm) from Barro et al. (2019). We corrected the CANDELS photometry for Galactic extinction using the Fitzpatrick (1999) curve, with A V = 0.025, as retrieved from the IRSA DUST tool 13 . To retrieve the X-ray products needed for our fits, we queried the Chandra Source Catalog (CSC) by performing a cone search in 1 around the source position, finding a unique match, 2CXO J033226.4−274035, within 0.41 . We utilized the level 3 CSC spectrum and response files for the single deepest (≈ 163 ks) observation of the source (ObsID 5019). To produce the X-ray photometry used for our fits, we subtracted the scaled background from the spectrum and grouped the X-ray counts into 15 log-spaced bins spanning 0.5-6.0 keV using Sherpa v4.13. We chose this energy range and binning such that the SNR in each bin is >2.
We modeled the resulting SED for J033226.49−274035.5 both with and without an X-ray model. In both cases, we used a stellar population with solar metallicity (i.e., Z = 0.02) and stellar age bins spanning 0-10 Myr, 10-100 Myr, 0.1-1 Gyr, 1-5 Gyr, and 5-5.6 Gyr (the age of the universe at z = 1.03). Both models also included the SKIRTOR UV-to-IR AGN model. Simultaneously constraining the viewing angle and optical depth of the torus is difficult, so we simplified the model by setting τ 9.7 = 7, the middle of the SKIRTOR model's allowed range. Additionally, since J033226.49−274035.5 is classified as a type 1 AGN, we implemented a prior on the cosine of the viewing angle to limit our models to a type 1 AGN. For the fit without an X-ray component, we allow the log of the integrated luminosity of the AGN model to vary between 11.0 and 13.0, with a uniform prior. To attenuate the UV-to-NIR emission (both the stellar and AGN emission that escapes the torus; see Section 2.2.1), we used the base Calzetti et al. (2000) curve. The dust attenuation was set to be in energy balance with the Draine & Li (2007) dust emission model, with U max = 3 × 10 5 , α = 2, and q PAH = 0.0047. We fix q PAH to the minimum allowed value for this example, since high-redshift galaxies like J033226.49−274035.5 are not expected to have strong PAH emission.
For the X-ray model in this example, Lightning automatically includes a stellar component when using an X-ray model, and we additionally use the qsosed X-ray model for the AGN X-ray emission component. X-ray absorption was modeled using the tbabs model, with Wilms et al. (2000) abundances and a Galactic HI column density fixed at N H = 9.19×10 19 cm −2 , as retrieved using the prop_colden tool in CIAO. We fit both models using the affine-invariant MCMC sampler, with an ensemble of 75 walkers running for 4 × 10 4 steps, assuming 10% model uncertainty. We adjusted the proposal distribution width parameter a to 1.8 to achieve accep-13 https://irsa.ipac.caltech.edu/applications/DUST/  tance fractions > 20%. The free parameters and associated priors for these fits are summarized in Table 4. We set Lightning to automatically generate the final chain portion of the posterior distributions from the MCMC chains and keep the last 1000 posterior samples, with the autocorrelation times indicating convergence of the runs. The photometry and resulting best-fitting models are shown in Figure 12. In Figure 13, we show a corner plot of the posterior distributions on the AGN parameters. Since the X-ray AGN model directly normalizes the UVto-IR AGN model, we calculated the equivalent L AGN from its UV-to-IR model to compare with the L AGN estimated from the model without X-ray data. In the bottom right corner of the corner plot, it can be seen that the addition of the X-ray emission results in highly consistent AGN viewing angle (cos i AGN ) and luminosity (L AGN ) estimates when including and excluding X-ray emission. Additionally, the incorporation of X-ray data also gives us estimates of M SMBH andṁ when using the qsosed model, and a separate independent estimate of the obscuration of the source in N H , which indicates  Figure 12. The X-ray-to-IR SED fit for J033226.49−274035.5. In the left panels, we show the instrumental X-ray spectrum (in terms of count rate density) with its best-fit model and residuals. In the right panels, we show the observed UV-to-IR SED (in terms of luminosity), its best-fitting model, and residuals. The best-fit model minimizes the total X-ray and UV-to-IR − log(Ppost). The Lightning X-ray model implementation can provide rudimentary X-ray spectral fits, and directly connect them to the UV-to-IR SED fit.
that this source is unabsorbed, consistent with its spectroscopic identification in the literature. In cases where MIR data is minimal or unavailable, adding the X-ray data can also give independent constraints on obscuration and add reliability to AGN IR luminosity estimates. Therefore, including X-ray observations when fitting the SED of an AGN can add valuable insights on the properties giving rise to the X-ray emission itself and place powerful constraints on the derived properties of the UV-to-IR component of the AGN.

Stellar X-ray Emission in an Inclined Galaxy
To demonstrate how X-rays emitted from the stellar binary population can be used to help constrain the SFR for an inclined galaxy (since X-rays are less sensitive to dust attenuation), we fit the global broadband photometry of the edge-on nearby galaxy, NGC 4631. For the UV-to-submillimeter photometry, we utilized the 30band SINGS/KINGFISH data presented in Table 2 of Dale et al. (2017). We then corrected the data for Galac-tic extinction before fitting, using the E(B−V ) values in Table 1 and A V -normalized extinction values in Table 2 of Dale et al. (2017).
For the X-ray photometry, we made use of the Chandra ACIS-I data for a single ≈58 ks observation (Ob-sID 797). These data were reduced and point-source catalogs were produced following the procedures detailed in Section 3.2 of Lehmer et al. (2019) using CIAO v.4.13. To obtain spectral constraints, we utilized the specextract routine to extract cumulative point-source and background spectral data. For the point-source spectrum, we chose to utilize circular apertures with radii that were four times the 90% encircled energy fraction and centered on the 22 X-ray detected point sources within the K s -band footprint of the galaxy (see Jarrett et al. 2003, for details on this region). We extracted the background spectrum from four large regions (circles with radii spanning 1-1.5 arcmin) outside the galactic footprint that were chosen to be free of bright X-ray detected point-sources. Using Sherpa we fit the background- with qsosed X-ray model without X-ray model J033226.49-274035.5 Figure 13. A corner plot of the AGN model parameters for J033226.49−274035.5, fit with and without an X-ray component as the blue and orange lines, respectively. The X-ray emission is modeled with the theoretical qsosed model and their addition to the SED results in highly consistent UV-to-IR AGN bolometric luminosity and inclination parameters compared to the fit without X-rays. In addition, the theoretical qsosed model also provides indirect estimates of the SMBH mass MSMBH and Eddington ratioṁ =Ṁ /Ṁ edd . subtracted point-source spectrum using a model that consisted of both thermal and absorbed power-law components (i.e., apec + tbabs × pow) to account for diffuse gas and the X-ray binary emission, respectively. Using the best-fit model, we calculated X-ray fluxes in the energy bands of 0.5-1, 1-2, 2-4, and 4-7 keV, with 10% uncertainty on the fluxes. When fitting this X-ray data with Lightning, we accounted for Galactic absorption by assuming a Galactic HI column den-sity of N H = 1.29 × 10 20 cm −2 , as derived from the prop colden tool in CIAO.
To show the effects of including X-rays and inclination-dependence, we modeled and fit the SED using four different permutations that include or exclude X-ray data with either the Calzetti et al. (2000) or inclination-dependent attenuation curves. For all four fits, we used a stellar population with solar metallicity (i.e, Z = 0.02) and SFH age bins of 0-10 Myr, 10-100 Myr, 0.1-1 Gyr, 1-5 Gyr, and 5-13.6 Gyr. Additionally, all fits included the Draine & Li (2007) dust emission model with the fixed values of U max = 3 × 10 5 and α = 2. As for the dust attenuation, which was set to be in energy balance with the dust emission, the models with the Calzetti et al. (2000) curve utilized the base curve extrapolated to the Lyman limit. For the fits with the inclination-dependent curve, we assumed the galaxy to be disk dominated and have a minimal contribution from the bulge (i.e., B/D = 0), a choice motivated by visual inspection. We further assume the youngest three age bins to be part of the young stellar population (i.e., r 0,old = 0 for 0 -1 Gyr and r 0,old = 1 otherwise). Finally, X-ray absorption was included using the tbabs model with Wilms et al. (2000) abundances for the fits that included X-rays.
To fit the four models to the SED, we utilized the affine-invariant MCMC algorithm with 5% model uncertainty, which we ran with 75 walkers for 10 4 and 5×10 4 trials for the Calzetti et al. (2000) and inclinationdependent models, respectively. The drastic increase in trials for the inclination-dependent models is required to reach convergence (i.e., autocorrelation times ≥ 50) of the cos i and τ f B parameters, since they are generally highly correlated (Doore et al. 2021). For all free parameters, we implemented the priors and limited initialization ranges as listed in Table 5. Since we know that NGC 4631 is an edge-on galaxy, we set a tabulated prior on cos i generated from the Monte Carlo method described in Section 3 of Doore et al. (2021), which converts an axis ratio into a distribution of inclination. The axis ratio and its uncertainty were retrieved from Hyper-Leda, which provides the axis ratio calclated from the 25 mag arcsec −2 B-band isophote.
With the described models and algorithm, we used Lightning to fit each model to the SED, assuming a luminosity distance to NGC 4631 of 7.62 Mpc as given in Dale et al. (2017). For each model, we set Lightning to automatically generate the final post-processed chain portion from the autocorrelation times and keep the final 2000 posterior samples. After confirming convergence of the fits from the autocorrelation time, we compared how the inclusion of X-rays influenced the derived properties.
In the right panels of Figure 14, we show the histograms of the resulting posterior distributions of the recent SFR of the last 100 Myr. From these distributions, each of the four models can be seen to have general agreement. However, the Calzetti et al. (2000) models have a stronger variation when including the X-rays compared to the inclination-dependent models, since the Calzetti et al. (2000) attenuation model, which assumes a uniform, spherical distribution of stars and dust, is too  simplistic for edge-on galaxies. Including the inclination dependence allows for a more accurate estimate of the SFR, with the inclusion of the X-rays increasing the precision of the estimate as would be expected when adding additional data. Further, the X-ray data rules out some higher SFR solutions (i.e., SFR > 8 M yr −1 ), as they become more unlikely with the X-ray data constraint.

Bayesian Sampling Code Comparison
To show how Lightning compares with other Bayesian sampling SED fitting codes, we fit the global broadband photometry of the face-on nearby spiral The total best-fit model spectrum to the observed SED for each of the four models. Inlaid in the SED plot is the posterior distributions of inclination (in terms of 1 − cos i, where 1 − cos i = 1 is edge-on) for the inclination-dependent models. This inlay shows that the models are correctly predicting a close to edge-on view. In all plots, the inclination-dependent model with and without X-rays and the Calzetti et al. (2000) model with and without X-rays are shown as the blue, green, orange, and red lines, respectively. While the best-fit model spectra are practically identical for each model, the resulting SFR distributions vary depending on the attenuation model and inclusion of X-ray emission.  (2007)  fitting codes currently available. However, we have chosen to compare directly with Prospector and BAGPIPES due to their inclusion of similar models and algorithms (e.g., non-parametric SFHs and Bayesian sampling algorithms) as Lightning in the interest of a "fair" compar- Prospector and BAGPIPES SFH M j e U (0, 10 12 ) · · · [10 7 , 10 7 ] · · · Note-U (a, b) indicates a uniform distribution from a to b. Fixed parameters have their value listed in the initialization range columns.
a The initialization range specified for Lightning.
b The initialization parameters specified for Prospector. The first value is the median starting point, and the second is the dispersion scale around that point.
c The nested sampling algorithm in BAGPIPES does not require initialization. However, we list this column to show the fixed parameter values.
d BAGPIPES normalizes the Calzetti et al. (2000) curve with A V , which we convert to τ V via τ V = 0.4 ln(10)A V .
e The stellar mass parameter in M . Prospector and BAGPIPES normalize their SFH bins to unit stellar mass, where Lightning normalizes to unit SFR. Therefore, the priors on the parameters are different.
ison. 16 While the codes can be set to have many matching components, there remain some differences between them that will generate differences in their results. We list the components of each code in Table 6 for ease in comparison and note that the differences in IMFs are likely to cause the most variation in results (Kennicutt & Evans 2012;Conroy 2013).
The UV-to-submillimeter photometry of NGC 628 used for this comparison was taken from the 30-band SINGS/KINGFISH data presented in Table 2 of Dale et al. (2017), which we corrected for Galactic extinction using the extinction values given in Table 1 and 2 of Dale et al. (2017). Since Prospector and BAGPIPES do not have a built in model uncertainty method like Lightning, we added in quadrature an additional 10% uncertainty to the quoted uncertainties in Dale et al. 16 For a more complete comparison of the most established SED fitting codes, we recommend Pacifici et al. (2022) and note that Lightning is not included in their comparison due to its more recent development.
(2017) to act as model uncertainties in our fits. We note that this is the typically utilized method in SED fitting to account for model uncertainties, and therefore it is a reasonable method for accounting for additional uncertainty not contained within the data. For the models in each code, we used a stellar population, as given in Table 6, with constant solar metallicity (i.e, Z = 0.02) and SFH age bins of 0-10 Myr, 10-100 Myr, 0.1-1 Gyr, 1-5 Gyr, and 5-13.4 Gyr. We attenuated the stellar emission in all codes using the original Calzetti et al. (2000) curve extrapolated to the Lyman limit. Finally, the dust attenuation was set to be in energy balance with the Draine & Li (2007) dust emission model, where α = 2, U max = 3×10 5 for Lightning, and U max = 10 6 for Prospector and BAGPIPES (the differences in U max have minimal effects on the models, see Section 2.4).
To fit the SED with each code, we utilized their Bayesian sampling algorithms to generate posterior distributions containing 2000 samples. For Lightning and Prospector, we ran their affine-invariant MCMC algo-rithm for 10 4 trials with 75 walkers, which was sufficient for each code to reach convergence using the autocorrelation times. The resulting chains were then postprocessed (i.e., burn-in removed and thinned) to the 2000 samples using the longest autocorrelation time of any parameter. 17 For BAGPIPES, we used the available MultiNest (Feroz et al. 2009(Feroz et al. , 2019 nested sampling algorithm, which we ran using 1000 live points. We note that we increased the number of live points from the default of 400 to 1000 after testing showed that we could get significant variation in the posteriors between runs when using fewer than 1000 live points for our chosen model. In Table 7, we list the free parameters in each model along with their utilized prior function. We note that we used the same priors across each code, except for the SFH parameters. Unlike Lightning, Prospector and BAGPIPES normalize their SFHs by stellar mass rather than SFR. Therefore, they require different priors to accommodate the change in normalization, which is expected to cause minimal to no effects on the fits due to the utilization of the same SFH bins. With the results from the fits, we first compared the computational performance of each code. All codes were run sequentially on the author's 2016 MacBook, which contains a 2-core, 1.2 GHz CPU. Lightning, Prospector, and BAGPIPES took 1279.7 s, 4864.2 s, and 920.0 s, respectively, to complete their fitting. While BAGPIPES can be seen to be almost 1.4 times faster than Lightning (which is 3.8 times faster than Prospector), it is important to note that this is due to the fewer likelihood evaluations needed by the nested sampling algorithm, which was designed (in part) to reduce the number of model evaluations required to produce a full sampling of the posterior (see, e.g., Feroz et al. 2009Feroz et al. , 2019. Where Lightning and Prospector each took 1.125 × 10 6 likelihood evaluations, BAGPIPES only performed an order of magnitude fewer (115,921 evaluations) to fit. Therefore, the difference in algorithms allowed for an overall similar fitting time as Lightning. However, in terms of likelihood evaluations per second (which is a better comparison of the practical speed of different SED fitting codes), Lightning is almost seven times faster than BAGPIPES, which is a result of its designed computational efficiency.
To show how the derived results compare at a base level, we show the observed photometry and best-fitting model spectra in Figure 15. In the lower residual panels for each fit, we also show the derived p-value for the  Figure 15. The total best-fit model spectra to the observed SED of NGC 628 and the associated residuals as generated by Lightning, Prospector, and BAGPIPES as the blue, green, and orange lines, respectively. For each code, the p-value from a PPC is shown in the upper left of each residual plot. In general, it can be seen that all three codes model the data well, both from the PPC and the residuals. fits calculated from a PPC. From the p-values and residuals, it can be seen that all codes appropriately model the data. It is interesting to note that the best-fit model spectra and resulting residuals from the Lightning and Prospector fits are highly similar, since their models only differ by the SSPs. However, BAGPIPES identifies a comparatively unique best-fit spectrum and residuals in the UV-to-NIR. This variation is expected, since the different IMF in BAGPIPES can create a substantial variation in the UV-to-NIR stellar emission models.
Finally, in Figure 16, we show the derived SFHs and posterior distributions for five commonly derived parameters from SED fitting: the (surviving) stellar mass (M ), the SFR of the last 100 Myr, the sSFR of the last 100 Myr, the V -band attenuation (A V ), and L TIR (i.e., the bolometric luminosity from 8-1000 µm). From these distributions in the lower left, it can be seen that The median SFH with the 16th and 84th percentile uncertainty range given as the offset vertical lines for each SED fitting code. For all plots, the results from Lightning, Prospector, and BAGPIPES are given as the blue, green, and orange lines, respectively. From these plots, it can be seen that the results from Lightning and Prospector are highly consistent, while the BAGPIPES results vary in consistency, especially for the SFH.
all three codes have derived parameters, except stellar mass, that are in excellent agreement. The variation in stellar mass between codes is expected and entirely due to the differences in SSPs and IMFs, which dictate the surviving stellar mass of the populations over time. As for the other parameters and SFHs, Lightning and Prospector have near identical results, which is a quality indicator given that both codes were run using almost identical models and were independently developed. As for the differences with BAGPIPES, these are a combination of the differences in fitting algorithms and IMFs, especially for the SFH. However, overall for both Lightning and SED fitting as a whole, it is reassuring that codes with different models and algorithms generally return derived parameters that are in statistical agreement.  Calzetti et al. (2000) Color excess, E(B − V ) 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.7, 0.9 mag Dust em.; Dale et al. (2014) α slope in dM dust ∝ U −α dU 1.

X-ray AGN Code Comparison
With the new additions to Lightning making it capable of fitting the X-ray-to-IR SEDs of AGNs, we examined how its derived AGN properties compared to the only other SED fitting code publicly available (at the time of publication) 18 also capable of modeling Xray-to-IR AGN emission, CIGALE 19 (Boquien et al. 2019;Yang et al. 2020Yang et al. , 2022. Unlike the comparison in Section 4.4.1, Lightning and CIGALE have notably different models and algorithms. For example, commonalities for both include the modified Calzetti et al. (2000) attenuation curve, Draine & Li (2007) dust emission models, SKIRTOR AGN templates, and power-law X-ray AGN models. However, Lightning implements nonparametric SFHs, while CIGALE mainly relies on parametric SFHs. Additionally, when fitting the models to data, Lightning can utilize either Bayesian sampling or maximum-likelihood inferencing, while CIGALE only utilizes a gridded Bayesian statistical inferencing method. Therefore, rather than attempting to use similar models and algorithms like we did in Section 4.4.1, in this example, we fit using the recommended settings of each code when modeling sources that include an X-ray AGN to see how the derived results compare.
For this comparison, we utilized the same X-ray AGN presented in Section 4.2, J033226.49−274035.5. The fit we performed in Section 4.2 already utilized our recommended models and the vast majority of our recommended model settings (see Table 4). The only differ-ence from the recommendation is fixing q PAH and limiting the prior range of cos i AGN . We retain these setting for this example, as q PAH would be a nuisance parameter if not fixed (as can be seen from the negligible contribution to the MIR by the dust emission in Figure 12) and the cos i AGN prior restrics modeling to type 1 AGN. Therefore, we utilized our results from the X-ray model fit in Section 4.2 for this comparison.
For the CIGALE settings, we utilized their recommended settings for an X-ray AGN using the example in their online documentation 20 . This model consisted of a delayed exponential SFH with a Chabrier (2003) IMF and metallicity of Z = 0.02, UV-to-IR AGN emission using the SKIRTOR templates with the Schartmann et al. (2005) intrinsic-disk type, X-ray AGN emission from a power-law model assuming Γ = 1.8, Calzetti et al. (2000) dust attenuation, and Dale et al. (2014) dust emission. The list of all gridded or non-default parameters used for each model are given in Table 8. We note that we deviated slightly from their recommended settings for two parameters in the model: θ (the AGN viewing angle, equivalent to i AGN in Lightning) and δ AGN (the deviation from the default UV/optical slope of the AGN model). We limit the grid points for θ to all possible options for type 1 AGN views and allow δ AGN to vary as is recommended when modeling type 1 AGN (Yang et al. 2022). Additionally, we set "lambda fracAGN" to define the AGN IR fraction as the integrated luminosity between observed-frame 5-1000 µm for ease in comparison with Lightning.
Utilizing the same observational data of J033226.49−274035.5 given in Section 4.2, we then fit the data using CIGALE. We do note that while we used the same UV-to-IR data, CIGALE requires X-ray data be input as absorption-corrected X-ray fluxes rather than the instrumental counts input into Lightning. This conversion to absorption-corrected fluxes requires one to either assume or fit a spectral model to the X-ray spectrum prior to performing the full SED fit, in order to include the X-ray data. For this example, we chose to fit the X-ray spectrum with an absorbed power-law (i.e. tbabs × pow) using Sherpa v4.13, finding a bestfit Γ = 1.96 and N H < 1 × 10 20 cm −2 . From this X-ray spectrum fit, we then calculated the absorptioncorrected fluxes for input into CIGALE in 15 log-spaced bands from 0.5-6.0 keV (the same energy bands used in Section 4.2).
Both codes were then run sequentially on the author's 2016 MacBook, which contains a 2-core, 1.2 GHz CPU. Lightning and CIGALE took ≈1 hr and ≈10 hr, respectively, to complete their fitting. While Lightning fit an order of magnitude faster in absolute terms, like in Section 4.4.1, it is better to compare the likelihood evaluations per time interval for an equivalent speed comparison. For this fit, Lightning performed 3 × 10 6 likelihood evaluations, while CIGALE performed a factor of four more (12,830,400 evaluations) to fit. Therefore, Lightning is approximately two times faster in terms of likelihood evaluations, which again is a result of its designed computational efficiency.
With the results from both fits, we compared the derived physical properties from Lightning to those derived by CIGALE. In Figure 17, we show the best-fit model SED for both codes, with the gray and magenta colored data points in the X-rays correspond to each codes input data, since Lightning had inputs in units of observed counts (which are converted to model dependent luminosities for plotting) and CIGALE had inputs in units of absorption-corrected fluxes. Additionally, since Lightning uses an MCMC algorithm, we generated the range of the top 64% of the best-fit models, in terms of the posterior, to create an uncertainty range on the model SEDs. With the uncertainty range, it can be seen that the best-fit model SED for both codes are highly consistent in the UV-to-IR, but differ in the X-rays. This difference in the X-ray is mainly caused by the difference in models (Lightning uses the more flexible qsosed model, while CIGALE uses a more rigid power-law) and input data (Lightning uses the observed counts, while CIGALE uses absorption-corrected  Figure 17. The total best-fit model spectra to the observed SED of J033226.49−274035.5 and the associated residuals as generated by Lightning and CIGALE as the gray and magenta lines, respectively. The gray and magenta colored data points in the X-rays correspond to each codes input data, since Lightning had inputs in units of observed counts (which are converted to model dependent luminosities for plotting) and CIGALE had inputs in units of absorption-corrected fluxes.
For the Lightning fit, we include each model component as in Figure 12. Additionally, we include the range of the top 64% of best-fit models from the MCMC posterior as the correspondingly colored shaded regions to show the uncertainty range on the component SEDs. For each code, we also show the best-fit χ 2 value in the center of each residual plot. In general, it can be seen that both codes model their data well from the best-fit χ 2 and the residuals.
fluxes). Therefore, differences are to be expected for the models in the X-rays. Finally, we list the commonly quoted parameters (mainly focusing on X-ray and AGN parameters) derived from each code in Table 9. These parameters include: the recent SFR over the last 100 Myr, the surviving stellar mass (M ), the total integrated IR luminosity of the dust emission model (L TIR ), the AGN torus inclination (i AGN ), the UV-to-IR AGN bolometric luminosity (L AGN ), the fraction of AGN IR luminosity to the total IR luminosity from 5-1000 µm (frac AGN ), and the slope of theL 2 keV −L 2500 relationship (α OX ). Since Lightning does not calculate frac AGN or α OX by default, we derived these values from the postpro-  (6): UV-to-IR AGN bolometric luminosity. Column (7): Fraction of the AGN IR luminosity to the total IR luminosity between observed-frame 5-1000 µm. Column (8): Slope of theL 2 keV −L 2500 relationship (α OX = −0.3838 log 10 (L 2500Å /L 2keV )). Column (9): SMBH mass (unique to Lightning). Column (10): SMBH accretion rate, Eddington rate normalized (unique to Lightning). Column (11): Modified optical slope power-law index of the SKIRTOR templates.
cessed results. We calculated frac AGN as the fraction of the model AGN IR luminosity to the total IR luminosity between observed-frame 5-1000 µm. For α OX , we first derived monochromatic luminosities at restframe 2 keV and 2500Åand computed it as α OX = −0.3838 log 10 (L 2500Å /L 2keV ) (Equation 5 in Yang et al. 2020).
From the data in Table 9, it can be seen that these commonly used properties are in excellent statistical agreement between the two codes. Of the derived properties, the SFRs and stellar masses have the largest uncertainty. This is caused by the uncertainty in the stellar models, since the type 1 AGN contributes a significant fraction of the emission at optical wavelengths. As for the AGN parameters, we find that only L AGN disagrees between codes at ≥1σ level, but their values are similar. The slight disagreement arises from the small uncertainties for both fits that are caused by the highly constraining 15 band X-ray data, which strongly limits the allowable AGN luminosity in the UV-to-IR. Comparatively, frac AGN is highly consistent between codes, indicating common contributions from the AGN to the total IR model. Overall, just like the comparison in Section 4.4.1, it is reassuring, for these results and SED fitting in general, that both Lightning and CIGALE return derived parameters that are in excellent statistical agreement despite using different models and algorithms.

SUMMARY AND PLANNED ADDITIONS
In this paper, we have presented the most recent version of the SED fitting code Lightning. The new version of Lightning contains a variety of models and algorithms that can be used to account for any combination of stellar, dust, and AGN emission in an observed X-ray to submillimeter SED. A brief review of each of these models and algorithms is as follows: 1. Stellar emission can be modeled using the SSPs from PÉGASE integrated over the age bins given by the user-defined non-parametric SFH. Stellar Xray emission from the XRBs is linked to the SFH using a power-law spectral model and the empirical parameterizations of L X /M with stellar age given in Gilbertson et al. (2022).
2. AGN emission can be modeled in the UV-to-IR using a subset of the SKIRTOR models. X-ray AGN emission can be modeled as (1) a simple rigid power-law spectra, which is tied to the UV-to-IR AGN emission using the Lusso & Risaliti (2017) L 2 keV −L 2500 relationship or (2) the physicallymotivated qsosed models from Kubota & Done (2018), which directly scale the UV-to-IR emission as a function of the mass and Eddington ratio of the SMBH.
3. Dust attenuation of the UV-to-NIR emission can be modeled using either a modified form of the Calzetti et al. (2000) curve or the inclinationdependent curve described in Doore et al. (2021).
Absorption of X-ray emission, when included, is modeled using either the tbabs or the Sherpa atten models.
4. Dust emission can be modeled using the Draine & Li (2007) model. When included, dust emission can be set to be in energy balance with the dust attenuation, which requires the bolometric luminosity of the dust emission to be equal to the bolometric luminosity of the attenuated light.
5. Algorithms for fitting the models to the data include both maximum likelihood and Bayesian methods. For the maximum likelihood method, Lightning uses the MPFIT implementation of the gradient-descent Levenberg-Marquardt algorithm. For the Bayesian methods, Lightning includes two MCMC algorithms, an adaptive Metropolis-Hastings algorithm from Andrieu & Thoms (2008) and an implementation of the Goodman & Weare (2010) affine-invariant algorithm.
With these models and algorithms, we presented different example applications of Lightning. These examples included (1) deriving spatially resolved stellar properties of M81 using an SED map, (2) demonstrating how the SMBH properties of an AGN can be derived by including X-ray emission, (3) exploring how X-ray emission and inclination-dependent attenuation can be used to constrain the SFR of an edge-on galaxy, (4) comparing the performance of Lightning to similar Bayesian sampling SED fitting codes (Prospector and BAGPIPES), and (5) comparing the X-ray-to-IR AGN properties derived from Lightning and CIGALE. From these examples, we clearly demonstrate the capabilities of Lightning and some of its potential uses.
In future updates to Lightning, we plan to expand our current PÉGASE stellar models. This will include adding additional IMF choices and allowing for a constant but continuous metallicity, which, like other SED fitting codes, we plan to have as an optional free parameter. Additionally, Lightning is currently restricted to using an exact redshift value if it is used as the distance indicator. We plan to allow for increased flexibility in redshift, where it can have an associated prior distribu-tion when fitting using a Bayesian method. This would allow for better propagation of uncertainty when using photometric redshifts, which can have large associated uncertainties.
Further, since Lightning was originally developed to be used in XRB population studies, we plan to add new SSPs like BPASS 21 (Eldridge et al. 2017;Stanway & Eldridge 2018) and/or POSYDON 22 (Fragos et al. 2022) that include binary population evolution. This would allow for more accurate stellar emission models, as binary populations can significantly influence the stellar emission (Eldridge & Stanway 2020). Additionally, binary stars are the progenitors of compact object binaries. Future versions of binary stellar population models may provide predictions for the observed X-ray binary luminosity function and its evolution with age and metallicity, based on the same prescriptions that govern the age and metallicity evolution of the stellar population. By adopting such models, we can self-consistently produce L X /M with our stellar population models rather than relying on empirical relations.
We acknowledge and thank the anonymous referee for their valuable comments that helped improve the quality of this paper. We gratefully acknowledge support from the NASA Astrophysics Data Analysis Program (ADAP) grant 80NSSC20K0444 (KD, EBM, RTE, BDL). KG is supported by an appointment to the NASA Postdoctoral Program at NASA Goddard Space Flight Center, administered by Oak Ridge Associated Universities under contract with NASA. ABZ supported by NASA award number 80GSFC21M0002. This work has made use of the Arkansas High Performance Computing Center, which is funded through multiple National Science Foundation grants and the Arkansas Economic Development Commission. We acknowledge the usage of the HyperLeda database (http://leda.univ-lyon1.fr).