ABSTRACT

Hot dust in the proximity of active galactic nuclei (AGNs) strongly emits in the near-infrared producing a red excess that, in type 2 sources, can be modelled to measure its temperature. In the era of high spatial resolution multiwavelength data, mapping the hot dust around supermassive black holes is important for the efforts to achieve a complete picture of the dust’s role and distribution around these compact objects. In this work, we propose a methodology to detect the hot dust emission in the proximity of type 2 AGNs and measure its temperature using K-band spectra (λc = 2.2 µm). To achieve this, we have developed nirdust, a python package for modelling K-band spectra, estimating the dust temperature, and characterizing the involved uncertainties. We tested synthetic and real spectra in order to check the performance and suitability of the physical model over different types of data. Our tests on synthetic spectra demonstrated that the obtained results are influenced by the signal-to-noise ratio (S/N) of the input spectra. However, we accurately characterized the uncertainties, which remained below ∼150 K for an average S/N per pixel exceeding 20. Applying nirdust to NGC 5128 (Centaurus A), observed with the Gemini South Telescope, we estimated a dust temperature of 662 and 667 K from Flamingos-2 spectra and 697 and 607  K from Gemini Near-Infrared Spectrograph (GNIRS) spectra using two different approaches.

1 INTRODUCTION

The Unified Model for Active Galactic Nuclei (Antonucci 1993; Netzer 2015) proposes orientation with respect to the observer as the fundamental difference between active galactic nucleus (AGN) types. The key to the unification is a nuclear structure of molecular gas and dust that absorbs the radiation in some lines of sight but not in others and was unresolved in the observations that drove the model (e.g. Mason 2015). As modern facilities become available, this structure can be studied at higher spatial resolutions, via infrared spectroscopy and interferometry and sub-mm observations, bringing clues to its true morphology. This structure, called ‘the torus’ due to its theoretical origin, not only absorbs radiation in order to display the differences between AGN types but also could play an important role in the accretion mechanism of the supermassive black hole (SMBH) as it is part of the inflow and outflow processes in the nuclear regions of galaxies.

Over the last decades, the model for the dusty absorber has suffered profound changes. The idea of a ‘clumpy torus’ was adopted after the homogeneity of the torus was ruled out on the basis of the high temperatures the dust reaches under such conditions (Krolik & Begelman 1988; Nenkova, Ivezić & Elitzur 2002; Nenkova et al. 2008). Later, the discovery of a polar component observed with mid-infrared (MIR) interferometry (e.g. Hönig et al. 2012, 2013; López-Gonzaga et al. 2014; Asmus, Hönig & Gandhi 2016; Asmus 2019; Stalevski, Tristram & Asmus 2019) resulted in a new model that incorporates a polar component as an outflowing wind and a thin annulus of high-temperature dust (∼1500 K) observed in near-infrared (NIR) interferometry data at the dust sublimation radius (Hönig 2019).Parallel to this, several authors have found massive and large equatorial components observing in the sub-mm regime, establishing that the dusty absorber can reach sizes up to 100 pc (e.g. Alonso-Herrero et al. 2018; Izumi et al. 2018; Combes et al. 2019; García-Burillo et al. 2021). Regarding the chemical composition of the dust, temperature plays a key role, as different grain compositions and sizes sublimate with increasing temperature. For instance, silicates are destroyed in regions with temperatures above 1200 K while graphite grains can survive up to 1900 K (Hönig & Kishimoto 2017). In this context, the temperature of the hot dust in the torus is a fundamental parameter that is related to the determination of the density, the composition, and the geometry of the torus.

Dust radiates at different wavelengths depending on its temperature. In the case of dust thermally heated by the accretion disc of an SMBH part of the emission is detected in the near-infrared and, in particular, is prominent in the K band (2–2.4 µm). Several authors have reported measurements of the temperature of hot dust in AGNs using NIR colours (e.g. Glass & Moorwood 1985; Alonso-Herrero et al. 1998), NIR interferometry (e.g. Kishimoto et al. 2011; GRAVITY Collaboration 2020; Leftley et al. 2021; Gámez Rosas et al. 2022), or NIR spectroscopy (e.g. Burtscher et al. 2015; Durré & Mould 2018; Gaspar et al. 2019, 2022; Riffel et al. 2022) and found temperatures ranging from 700 to 1500 K.Burtscher et al. (2015) measured the temperature of the hot dust component for 51 near AGNs in 1 arcsec apertures and found a mean of 887 ± 68 K for type 2 AGNs and 1292 ± 46 K for type 1 AGNs. Regarding the extension and temperature distribution of this hot dust, Gaspar et al. (2019, 2022) found T ∼ 1000 K dust in scales of tens of pc and almost constant in temperature in two type 2 Seyfert galaxy nuclei. Gámez Rosas et al. (2022) measured dust temperatures using MIR interferometry with Multi AperTure mid-Infrared SpectroScopic Experiment (MATISSE) at Very Large Telescope Interferometer (VLTI) instrument in several sub-pc regions around the SMBH of NGC 1068, the archetypical type 2 Seyfert galaxy. They found temperatures ranging between 800 and 1500 K in the different regions that ‘do not decrease systematically with distance from the inferred black hole position’. Finally, Riffel et al. (2022) found that hot dust emission extended outside the unresolved core and stated that the emission inside ∼40 pc is dominated by hot dust in a sample of 18 nearby Seyfert galaxies. The variety of these dust detections evidence that more temperature measurements in high angular resolution infrared data are necessary to improve our understanding of the extension, distribution, and origin of the hot dust in the proximities of AGNs.

The resolved NIR spectroscopic observations currently available with ground-based 8–10 m class telescopes are suitable for determining the extent and temperature distribution of the hot dust component in AGNs, as they can resolve the central regions of these objects at resolutions of tens of pc or less in nearby galaxies. Furthermore, the Mid-Infrared Instrument (MIRI) (5–28 µm) and Near Infrared Spectrograph (NIRSpec) (0.6–5 µm) instruments onboard the JWST will enable spectroscopic mapping of dust and molecular gas in the central regions of AGNs with unprecedented sensitivity and spatial resolution.

This work introduces nirdust, a python-based tool designed for detecting and measuring the temperature of hot dust near type 2 AGNs. nirdust utilizes K-band rest-frame spectra to fit the temperature of a blackbody component emitted by the hot dust in the NIR range. The package offers a range of functionalities, including spectrum storage, pre-processing, and fitting. nirdust’s application programming interface facilitates seamless integration with astropy (Astropy Collaboration 2013), particularly with specutils, enabling users to incorporate nirdust into their broader code base when working with NIR spectra.

In Section 2, we present the fundamental physics basis for the code. In Section 3, the functionalities and quality of the package are described along with an application example, and in Section 4 the results of tests to evaluate the performance of nirdust are presented. The summary and future work are in Section 5.

2 MODELLING THE HOT DUST EMISSION

For type 2 AGNs, in the absence of a non-thermal featureless continuum, the nuclear spectrum continuum in the K band is conformed by two elements: the stellar population spectrum and the emission by hot dust (Thatte et al. 1997; Gratadour et al. 2003, and references therein). This is easily appreciated when observing the slope of the nuclear spectrum continuum that remains constant or even rising when a pure stellar spectrum is expected to fall (e.g. Ferruit et al. 2004; GRAVITY Collaboration 2020), i.e. a red excess is present. This means that if the stellar spectrum can be computed, then the dust emission can be isolated. This hot dust component not only is detected in the nuclear spectrum but also can be found in more external regions in high angular resolution spectra (e.g. Gaspar et al. 2022; Riffel et al. 2022).

At moderate spectral resolution (R = 1000–3000), the dust emission can be modelled by a blackbody function of single temperature. In the real case, the emission could be composed by a continuum of temperatures, but for this kind of spectra experience has shown that using a single temperature is sufficient (e.g. Durré & Mould 2018; Gaspar et al. 2019).

The blackbody radiation can be described by Planck’s law as a function of frequency ν, and parametrized by its temperature T and an amplitude A as

(1)

where h, k, and c correspond to the Planck constant, Boltzmann constant, and speed of light, respectively. The scale factor A is used to model the intensity amplitude of the blackbody radiation for a given temperature T. It is worth noticing that the introduction of this factor does not affect the location of the radiation peak or the general shape of the blackbody function. The dust emission can then be modelled by the free parameters A and T.

To model the target spectrum, where the hot dust emission is expected to be found, we use a linear combination of the stellar component and a blackbody function as in Durré & Mould (2018):

(2)

Here, St is the target spectrum, Ss is the reference spectrum that accounts for the stellar emission, BT is the blackbody function of temperature T, and γ′ is a constant emission that represents the sky background (it is recommended that the sky background is previously subtracted to make γ′ as closer to zero as possible). α and β′ are scale factors of the stellar and blackbody components, respectively. The reformulation of β′ and γ′ is introduced for the technical purpose of exploring the parameter space with a logarithmic distribution (see Section 3 for more details on the technical discussion). Thus, the parameter space is defined by |${\boldsymbol{\theta } \equiv \lbrace T, \alpha , \beta , \gamma \rbrace}$|⁠.

In order to model the hot dust emission with equation (2) spectra do not necessarily have to be flux calibrated. This is important because the uncertainties introduced in the flux-calibration process can reach |$30\!\!-\!\!50{{\ \rm per\, cent}}$| in the NIR due to the fast changes in atmospheric conditions that occur at these wavebands. The telluric correction (essential when studying spectral continuum) eliminates the shape of the atmosphere transmission and the instrumental response from the spectrum and hence the flux calibration is reduced to a multiplicative operation. Given the linear combination that models the target spectrum, this multiplicative factor is absorbed by the β parameter. The temperature information does not change after the multiplicative operation.

2.1 The reference spectrum

As stated above, in order to isolate the hot dust spectrum present in the nuclear region of a type 2 AGN, it is necessary to provide the stellar spectrum of the same region. The hypothesis working here is that Ss is the same as the target spectrum stellar component, in other words, that the stellar population that emitted Ss is the same stellar population that inhabits the region that emitted St.

Ss is, of course, not possible to isolate from the target spectrum without knowing the dust emission, but there are some hypotheses under which this stellar spectrum can be represented by different proxies. The main hypothesis is that in the NIR the late-type population emission is dominant in comparison with the emission by young stars. This allows, for instance, for the nuclear stellar population to be represented by a spectrum or a mix of spectra of late-type stars. Moreover, the late-type population is more homogeneously distributed than the young population that is usually located in clusters or ‘clumps’ as can be seen, for example, in fig. 1 of Lin et al. (2013), where the same galaxy is shown in different wavebands. Taking this into account, it could be considered that a spectrum of the stellar population extracted at a certain distance from the nucleus can serve as reference for the nuclear stellar population.

One crucial consideration of this approximation is the frequent occurrence of starburst activity intertwined with AGN emissions in galactic nuclei. Depending on the intensity and duration of the starburst, deviations from the external stellar population used as a reference can emerge. In this context, the choice of width and location for extracting the reference spectrum becomes pivotal in accurately emulating the nuclear stellar population. Notably, in cases where high spatial resolution data are available, analysing the slopes of the continuum at various radii can offer insights into where the stellar population is most representative of the overall population and where regions of hot dust are prevalent. If a multitude of spectra are extracted from distinct radii, a meticulous examination of the spectral slopes can reveal a radius at which the slope consistently remains unchanged as one moves to larger radii. Once this phenomenon is observed, it becomes reasonable to infer that these particular spectra are indeed representative of the underlying spectral population. This approach provides a valuable tool for pinpointing radii where the stellar population’s characteristics remain relatively stable and unaffected by extreme local conditions such as the presence of hot dust. An illustrative example of this concept can be found in fig. 9 of Gaspar et al. (2022), where the variations in the slope of the spectral continuum are evident between the nucleus and two distinct radii. Interestingly, in these two radii, the spectral slope remains consistent. For a comprehensive set of spectra illustrating these findings, the interested reader is directed to the appendix of the aforementioned publication. In Section 4.2 of this work, we elaborate on the procedure employed to select the reference spectrum for a specific example presented herein, corresponding to the galaxy NGC 5128.

In summary and under the aforementioned hypothesis, Ss could be represented by a reference spectrum in three different ways:

  • A spectrum of the same galaxy extracted in the nuclear region but far enough from the nucleus to avoid the hot dust component. (moderate to high spatial resolution data)

  • A spectrum of a late-type star.

  • A nuclear stellar template of a non-active galaxy of similar type to the host galaxy under study.

It is important to acknowledge that nirdust has been exclusively tested using the first and second scenarios. The first one, which involves analysing a spectrum extracted near the nucleus of the same galaxy, has been previously employed in other studies (e.g. Durré & Mould 2018; Gaspar et al. 2019, 2022) and has also been successfully implemented as an alternative in situations where the population synthesis technique does not yield satisfactory results (Dumont et al. 2020). The second scenario, which was applied by Burtscher et al. (2015), and the third scenario, which remains hypothetical, require further testing to determine their suitability for modelling the stellar population of AGNs using nirdust.

3 NIRDUST

The entire experiment of this work is supported by nirdust. nirdust is an object-oriented package (Ram et al. 2003) that is fundamentally designed around the factory-method pattern (Gamma et al. 1995), so the easiest way to create the available objects on the package is usually some function or method. The core features of the project are centred around the NirdustSpectrum class, serving as a comprehensive model of a spectrum. This class is capable of segmenting the spectrum within designated boundaries, transforming the spectral axis into frequencies, standardizing intensity levels, deriving noise from a region defined by the user, and implementing selective masking on distinct portions of the spectrum. In addition, nirdust offers additional pre-processing capabilities for line spectrum construction and resampling of two spectra to match the same resolution. Finally, the main functionality of the package is also available, which consists of estimating and characterizing the dust temperature, this algorithm is explained in detail in Section 3.1.

To guarantee accurate results, nirdust maintains a strict adherence to the Python Enhancement Proposal 8 (PEP-8) coding standard (Van Rossum, Warsaw & Coghlan 2001) and employs the Flake-8 tool for seamless code consistency. With a comprehensive suite of 93 unit tests (Jazayeri 2007), it meticulously validates software components, accommodating python versions spanning from 3.8 to 3.11, and achieving an impressive 99 per cent code coverage (Miller & Maloney 1963).

A detailed description of the use of nirdust, the installation process, the bug report procedure, and a more in-depth discussion of the project quality can be found in Appendix  B.

3.1 Fitting procedure

The model given by equation (2) is fitted using both spectra, target and reference, using a basin-hopping algorithm (Wales & Doye 1997). This is a two-phase method that combines a global minimization algorithm with local minimization at each step. When the local optimization routine finds a local minimum, the basin-hopping algorithm randomly perturbs this local solution and attempts a new local optimization. This technique has provided useful results in many problems where the likelihood surface is hard to explore and simpler minimization algorithms tend to get stuck in local minima.

The parameter space of the problem is given by four free parameters |$\boldsymbol{\theta } \equiv \lbrace T, \alpha , \beta , \gamma \rbrace$|⁠. We consider the case where observations, D, are independent and identically distributed, sampled from the same Gaussian probability function |$f(D| \boldsymbol{\theta })$|⁠, where each parameter is described by a normal distribution |$\mathcal {N}(\mu , \sigma)$|⁠. Then, the likelihood function can be expressed as

(3)

where the index i runs through all observed points in the spectrum. The explicit expression takes the form

(4)

where σ corresponds to the uncertainty value of the observations. Ideally, one should use the covariance matrix, which is not always available in spectroscopic studies. We estimate the uncertainty σ using the standard deviation of the observed data to the fitted continuum.

4 APPLICATION ON SIMULATED AND REAL SPECTRA

In this section, we present two scenarios designed to test the accuracy of the model described in Section 2 and the capability of nirdust’s implementation to recover the true values of its parameters.

We consider the following cases: the use of controlled synthetic spectra, and the real case of NGC 5128 (Centaurus A or Cen A), which also serves as an independent comparison with the results reported by Burtscher et al. (2015) for this galaxy.

The goal of these tests is to evaluate the impact of the characteristics of the input spectra and the masking procedure on the estimated parameters, particularly the temperature. Moreover, the tests over synthetic spectra intend to serve as a guide to estimate the uncertainties involved when using nirdust on real spectra.

4.1 Synthetic spectra: accuracy in the presence of noise

For this test, we built multiple synthetic models that represent the target and reference spectra, where we vary the following conditions of interest to assess their impact.

  • The signal-to-noise ratio (S/N) of both spectra.

  • The percentage and location of points removed when masking the spectral features.

We simultaneously consider a wide range of temperature values, ranging between 500 and 1600 K. For simplicity, we modelled the stellar population of the reference spectrum as a decreasing linear function. The target spectrum was then modelled according to equation (2), where we adopted fixed parameter values of α = 3.5, β = 8.3, and γ = −3.3. However, we have studied the behaviour of the results presented here when these parameters are varied along the default bounds range for fit_blackbody() and found that the results remain unaltered. In all cases, a different Gaussian noise (different random seed) is added to each spectrum to obtain the desired S/N ratio. The spectral resolution of both spectra is set at 3.5 Å pix−1 and they consist of 541 spectral points. Throughout all this work, we have consistently measured the S/N ratio using nirdust’s compute_noise class method. This method calculates the noise by determining the standard deviation of the continuum-subtracted spectrum within a user-defined wavelength interval. This interval must be carefully chosen to represent the noise characteristics of the entire spectrum. Subsequently, the mean signal within this interval is divided by the calculated noise value yielding an average S/N ratio per pixel.

In Fig. 1, the results for the first test are presented. In this case, the S/N of both synthetic spectra are equal and varied from 20 to 500 for each value of temperature, with a higher sampling for S/N < 100. The fitting procedure was run with higher number of iterations (niter = 700) instead of the default value (niter = 200) of iterations to achieve convergence in the low S/N points. The deviation of the fitted temperature from the real value, ΔT, for different S/N is shown in the top left panel. In the top right, bottom left, and bottom right panels, the variations for the α, β, and γ parameters are, respectively, presented. At high S/N (≳150), the temperature is accurately estimated for temperatures below 900 K, but for higher temperatures a systematic and constant overestimation of ≲50 K is obtained. However, with decreasing S/N (20–50), lower temperatures show an increase in the overestimation reaching ΔT ∼ 100 K. As can be seen in the other three panels, α and β follow the same trend, remaining constant for S/N ≳ 150 and deviating for lower S/N in the cases of the lower temperatures (≲900 K). The γ parameter remains constant along the entire T and S/N ranges and represents a minor contribution to the total flux (⁠|$\lt\!\! 1{{\ \rm per\, cent}}$|⁠).

Results of the fitting procedure when varying the S/N of both the target and the reference spectra for different temperatures. The legend shows the temperatures of the blackbody used to construct the synthetic target spectrum in each case. The vertical axis in all panels shows the respective difference of the estimated value of a parameter minus its true value. The horizontal axis corresponds to the S/N ratio of the target spectrum that is equal to the S/N ratio of the reference spectrum in this case.
Figure 1.

Results of the fitting procedure when varying the S/N of both the target and the reference spectra for different temperatures. The legend shows the temperatures of the blackbody used to construct the synthetic target spectrum in each case. The vertical axis in all panels shows the respective difference of the estimated value of a parameter minus its true value. The horizontal axis corresponds to the S/N ratio of the target spectrum that is equal to the S/N ratio of the reference spectrum in this case.

Summarizing, for all values of temperature we find that the uncertainty remains below 150 K with a tendency to decrease for higher temperatures, as one would intuitively expect. To achieve an uncertainty below 50 K, a target spectrum with S/N ≳ 150 is required. It is important to remark that in all cases with S/N ≥ 20 the uncertainty corresponds to a systematic overestimation. The user can consider this as a ‘correction factor’ that can be applied to the obtained temperature. An example of this correction is presented in Section 4.2.

Three additional cases where the S/N of the reference spectrum is lower than the S/N of the target spectrum are presented in Appendix  A.

In Fig. 2, we present results regarding the application of masks used to remove spectral features or high-noise regions. The test consisted of computing the difference between the fitted temperature and the real temperature of the synthetic spectrum for different masking scenarios. Through this analysis, our objective was to discern the presence of spectral regions and/or a minimum of available spectral points crucial for accurate dust temperature determination, emphasizing the necessity for meticulous consideration when implementing masking strategies. The synthetic spectra are the same as used in the previous test but the S/N is fixed at a high value of 800 to minimize the uncertainties introduced by the Gaussian noise.

Uncertainty in the estimated temperature when varying the percentage of points removed from both the target and the reference spectra for different temperatures using mask_spectrum(). Four masks were implemented: four equidistant and equi-sized intervals distributed along the spectral range (top left panel); one interval in the blue part of the spectrum at 2.07 µm (top right panel); one interval in the red part of the spectrum at 2.24 µm (bottom left panel); and one interval in the central part of the spectral range at 2.15 µm (bottom right panel). In all cases, the percentage of removed points is the sum of the masked intervals. The temperatures shown in the legend are the temperatures of the blackbody model of the synthetic target spectrum in each case.
Figure 2.

Uncertainty in the estimated temperature when varying the percentage of points removed from both the target and the reference spectra for different temperatures using mask_spectrum(). Four masks were implemented: four equidistant and equi-sized intervals distributed along the spectral range (top left panel); one interval in the blue part of the spectrum at 2.07 µm (top right panel); one interval in the red part of the spectrum at 2.24 µm (bottom left panel); and one interval in the central part of the spectral range at 2.15 µm (bottom right panel). In all cases, the percentage of removed points is the sum of the masked intervals. The temperatures shown in the legend are the temperatures of the blackbody model of the synthetic target spectrum in each case.

Four cases were considered:

  • ‘Four intervals’ mask: a mask consisting of four equidistant and equi-sized intervals distributed along the spectral range (top left panel).

  • ‘Blue interval’ mask: a mask consisting of one interval in the blue part of the spectrum at 2.07 µm (top right panel).

  • ‘Red interval’ mask: a mask consisting of one interval in the red part of the spectrum at 2.24 µm (bottom left panel).

  • ‘Central interval’ mask: a mask consisting of one interval in the central part of the spectral range at 2.15 µm (bottom right panel). In all cases, the percentage of removed points is the sum of the masked intervals and increases between 0 and 60 per cent.

In all cases, the uncertainty in the estimated temperature is mainly constant until a 30 per cent of mask percentage is reached. Taking into account that the S/N of the spectra involved in this test is 800, it can be considered that this initial overestimation of the temperature is attributable to the S/N effect as analysed in the previous test. As the masking percentage increases, however, the uncertainty in the obtained temperature changes in all the cases except the four intervals mask (top left panel). In the other three cases, the uncertainty increases and even changes sign. This is expected since the four intervals mask is the one that keeps higher amounts of points distributed along the spectral range and hence is expected to retain the information of the continuum slope. Therefore, if the mask involved consists of one large interval placed anywhere along the spectral range, it is not recommended to mask beyond 30 per cent of the points.

4.2 Real spectra: Centaurus A

We have used three different sets of spectra of NGC 5128 (Cen A), obtained with different instruments and focal plane sampling techniques, to measure the hot dust temperature and compare it to the temperature reported by Burtscher et al. (2015) of (796 ± 1) K for the hot dust in Cen A nucleus measured in a 1 arcsec diameter aperture. Their spectra were selected from the European Southern Observatory's Spectrograph for INtegral Field Observations in the Near Infrared (ESO SINFONI, Eisenhauer et al. 2003) data base and correspond to Integral Field Units (IFU) observations in K band with seeing of 0.6 arcsec taken in 2005 March–April (Neumayer et al. 2007).NGC 5128 was chosen as the primary subject of our comparative study due to the availability of two extra suitable data sets for direct analysis obtained with Flamingos-2 and GNIRS both at the Gemini Observatory. Additionally, Burtscher et al. kindly provided us the SINFONI cube of NGC 5128 they used in their work and the spectrum of the star they used as reference spectrum for the stellar population HD 176617 (private communication). We adopted a distance of 3.84 Mpc for NGC 5128 (Rejkuba 2004) to calculate the linear radii in pc.

For this comparative analysis, we used two different approaches for the reference spectrum in order to compare them: on the one hand an off-nuclear spectrum of the same galaxy as described in Section 2.1 (Reference 1), and on the other hand, the spectrum of the star used in Burtscher et al. (2015) (Reference 2).

We extracted target and reference spectra from our two data sets:

  • Klong long-slit spectra taken with Flamingos-2 (Eikenberry et al. 2008; Gomez et al. 2012) at Gemini South with PA = 151° (semimajor axis direction), on 2021 February 21. The percentage of masking implemented is of 30 per cent and the mask consists of eight narrow intervals distributed along the band with bigger presence in the blue part of the spectra. Taking this into account and considering the behaviours of similar masks presented in Fig. 2 (the cases of masks composed of four intervals and one blue interval), the uncertainty is expected to be dominated by the S/N of the spectra (see Section 4.1).

  • A Klong data cube taken with GNIRS (Elias et al. 2006) at Gemini South in 2005 (Díaz et al. 2021; Program ID GS-2005A-Q-38). The percentage of masking is of 10 per cent and the mask consists of one interval in the red part of the spectrum. For this kind of mask, the uncertainty is expected to be only attributed to the S/N.

In the top panel of Fig. 3, the Flamingos-2 spectra are displayed for different radii, spanning the range from 0 to 43 pc, on both sides of the slit. The nuclear spectrum, designated as the target spectrum, is highlighted in red, while the reference spectrum is depicted in black. Notably, the blue spectra extracted from various radii exhibit a consistent slope of the spectral continuum across all spectra except the nuclear one. This means that the excess punitively produced by the hot dust emission is present only in the nuclear spectrum and hence all the other blue spectra displayed in the figure have the potential to serve as a reliable reference spectrum. Our choice, however, was directed towards one of the more distant spectra, thereby ensuring a substantial separation from the nucleus. For the GNIRS data set, shown in the bottom panel of Fig. 3, a distinct approach was necessary due to the instrument’s limited field of view. Instead of employing annular apertures, circular apertures were extracted from different regions surrounding the nucleus at a fixed distance of 2.25 arcsec. Integrating these circular aperture extractions yielded the reference spectrum. As a result of this methodology, Fig. 3 exclusively showcases the target spectrum and the reference spectrum. It is imperative to note that the validation of the reference spectrum’s integrity was supported by the Flamingos-2 spectra observations, ensuring the absence of red excess beyond the nucleus.

Top panel: Flamingos-2 spectra extracted at different radii (r) at both sides of the slit (SE and NW directions). The target and reference spectra are highlighted. Bottom panel: GNIRS spectra extracted at the nucleus (target) and at 42 pc of radius (reference); in this case, the reference spectrum is a combination of spectra extracted at 42 pc in several regions of the IFU field, for more details see the text. The dotted line is a guide to compare the target and reference spectra slopes.
Figure 3.

Top panel: Flamingos-2 spectra extracted at different radii (r) at both sides of the slit (SE and NW directions). The target and reference spectra are highlighted. Bottom panel: GNIRS spectra extracted at the nucleus (target) and at 42 pc of radius (reference); in this case, the reference spectrum is a combination of spectra extracted at 42 pc in several regions of the IFU field, for more details see the text. The dotted line is a guide to compare the target and reference spectra slopes.

The characteristics of the four extracted spectra, along with the SINFONI spectra from Burtscher et al. (2013), are summarized in Table 1. The table includes information on the masking percentage, the temperature obtained using nirdust, and the expected uncertainty calculated from the S/N of the two spectra involved in each fit. The table is divided into three parts corresponding to the instrument used to capture the target spectrum. For the Flamingos-2 and GNIRS data, there are two reference spectra: one extracted from the same galaxy (Reference 1) and the SINFONI spectrum of HD 176617 (Reference 2). A good fit for the SINFONI target spectrum could not be obtained using a reference spectrum extracted from the same data cube. Therefore, this result is not included in the table.

Table 1.

Description of the three data sets used to measure the dust temperature with nirdust. The obtained results are presented in column 7. The masking percentage is only indicated for the reference spectrum as it may not be the same for all the references. The target spectrum presents the same masking percentage as the reference in each case.

SpectrumData typeS/NApertureRadius ofMaskingTemperatureExpected
diameter (arcsec)extraction (arcsec) percentage(K)ΔT
Flamingos-2
TargetLongslit9510
Reference 1Longslit721.083.343066225–85
Reference 2 (HD 176617)IFU40Point source3066790–110
GNIRS
TargetIFU871.05
Reference 1IFU741.052.25106975–75
Reference 2 (HD 176617)IFU40Point source060790–110
SINFONI
TargetIFU20510
Reference (HD 176617)IFU40Point source30649100–110
SpectrumData typeS/NApertureRadius ofMaskingTemperatureExpected
diameter (arcsec)extraction (arcsec) percentage(K)ΔT
Flamingos-2
TargetLongslit9510
Reference 1Longslit721.083.343066225–85
Reference 2 (HD 176617)IFU40Point source3066790–110
GNIRS
TargetIFU871.05
Reference 1IFU741.052.25106975–75
Reference 2 (HD 176617)IFU40Point source060790–110
SINFONI
TargetIFU20510
Reference (HD 176617)IFU40Point source30649100–110
Table 1.

Description of the three data sets used to measure the dust temperature with nirdust. The obtained results are presented in column 7. The masking percentage is only indicated for the reference spectrum as it may not be the same for all the references. The target spectrum presents the same masking percentage as the reference in each case.

SpectrumData typeS/NApertureRadius ofMaskingTemperatureExpected
diameter (arcsec)extraction (arcsec) percentage(K)ΔT
Flamingos-2
TargetLongslit9510
Reference 1Longslit721.083.343066225–85
Reference 2 (HD 176617)IFU40Point source3066790–110
GNIRS
TargetIFU871.05
Reference 1IFU741.052.25106975–75
Reference 2 (HD 176617)IFU40Point source060790–110
SINFONI
TargetIFU20510
Reference (HD 176617)IFU40Point source30649100–110
SpectrumData typeS/NApertureRadius ofMaskingTemperatureExpected
diameter (arcsec)extraction (arcsec) percentage(K)ΔT
Flamingos-2
TargetLongslit9510
Reference 1Longslit721.083.343066225–85
Reference 2 (HD 176617)IFU40Point source3066790–110
GNIRS
TargetIFU871.05
Reference 1IFU741.052.25106975–75
Reference 2 (HD 176617)IFU40Point source060790–110
SINFONI
TargetIFU20510
Reference (HD 176617)IFU40Point source30649100–110

In Fig. 4, we present the nirdust fitting results for the five combinations of spectra described in Table 1. We obtained temperatures of 662 K (Reference 1) and 667 K (Reference 2) from the Flamingos-2 spectrum and of 697 K (Reference 1) and 607 K (Reference 2) from the GNIRS data cube. As previously mentioned, the synthetic spectra tests of Section 4.1 showed that the S/N of the spectra causes an overestimation of the temperature. We present the overestimation of the temperature estimated for each pair of spectra in column 8 of Table 1. For the SINFONI spectrum, the obtained dust temperature is 649 K, which is 147 K lower than the temperature of 796 K reported by Burtscher et al. (2015). Nevertheless, all temperatures obtained with nirdust are consistent with each other within a typical uncertainty of ∼100 K. Since both the target and reference spectra from SINFONI are the same as in Burtscher et al. (2015), this difference in temperature may arise from differences in the methodology, specifically the equation for the linear combination used to perform the fitting. Despite this difference, it is worth noting that all dust temperatures measured by nirdust and by Burtscher et al. (2015) fall in the range of relatively low temperatures that have been previously reported for type 2 Seyfert galaxy nuclei such as the one hosted by NGC 5128 (Riffel et al. 2009).

nirdust fitting results for the five combinations of spectra described in Table 1. The first row corresponds to the Flamingos-2 target spectra, the second row to the GNIRS target spectra, and the third row to the SINFONI target spectra. For the first two rows, the left column corresponds to the reference spectrum extracted from the same galaxy and the right column corresponds to the cases where the spectrum of HD 176617 was used as reference. In the five cases, the plots are the actual output of nirdust, to see the obtained temperature reference to Table 1.
Figure 4.

nirdust fitting results for the five combinations of spectra described in Table 1. The first row corresponds to the Flamingos-2 target spectra, the second row to the GNIRS target spectra, and the third row to the SINFONI target spectra. For the first two rows, the left column corresponds to the reference spectrum extracted from the same galaxy and the right column corresponds to the cases where the spectrum of HD 176617 was used as reference. In the five cases, the plots are the actual output of nirdust, to see the obtained temperature reference to Table 1.

More extensive experiments on the impact of the reference spectrum using different options and data sets are necessary in order to further constrain the uncertainties involved in the estimation of the dust temperature using the method described in Section 2.

5 SUMMARY AND FUTURE WORK

In this work, we have presented an analysis of the proposed model to describe the NIR emission of hot dust in type 2 AGNs. We designed a set of multiple synthetic spectra that allowed a thorough analysis of the behaviour of the mathematical model that describes the physics of the emitting dust. This study provided useful insights at the moment of analysing real spectra from Flamingos-2 (Gemini South), GNIRS (Gemini North), and SINFONI (VLT).

The results presented here can be summarized in the following key points:

  • The parametrization choice of the model in terms of |$\boldsymbol{\theta } = \lbrace T, \alpha , \beta , \gamma \rbrace$| and a basin-hopping optimization algorithm are suitable to detect the hot dust component in real spectra.

  • The tests performed over synthetic spectra show that estimated temperature presents a systematic error whose value depends primarily on the S/N of the target spectrum and can be accurately characterized for each data set. In most cases, the error of the temperatures is below 100 K for S/N ≥ 20; for the extreme cases of low T and low S/N, the error can rise up to 150 K.

  • The five dust temperatures obtained for NGC 5128 (Cen A) using three independent data sets from Flamingos-2 (662 and 667 K), GNIRS (697 and 607 K), and SINFONI (653 K) and two different reference spectra for the first two data sets are consistent within a typical uncertainty of ∼100 K. The obtained temperatures for the three data sets fall in the low-temperature regime for type 2 AGN NIR emitting dust.

  • The comparison of the dust temperatures obtained with nirdust for the SINFONI spectrum used by Burtscher et al. (2015) for NGC 5128 yields a difference of 147 K with their reported temperature of 796 K.

To achieve these results, we developed nirdust, an open-source python package that uses K-band rest-frame spectra to detect hot dust in type 2 AGNs and measure its temperature. This package is open to the astronomical community and is an effort to provide a standardized and reproducible procedure to measure the temperature of hot dust present in these objects. The performance of nirdust in the presence of Gaussian noise is satisfactory and the uncertainties involved can be estimated.

As a continuation of this work, we intend to characterize the impact of different reference spectrum choices on the estimated temperature. New features will be added, like an optional power-law function in the model to account for any scattered disc accretion emission present in the target spectrum.

ACKNOWLEDGEMENTS

The authors would like to thank Leonard Burtscher and Orban De Xivry Gilles for generously providing the NGC 5128 data they used in their work as these data significantly contributed to our analysis. The authors also thank their families and friends, OAC and IATE astronomers for useful comments and suggestions, and the anonymous referee for a deep reading and useful suggestions. This work was partially supported by the Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET, Argentina), the Secretaría de Ciencia y Tecnología de la Universidad Nacional de Córdoba (SeCyT-UNC, Argentina), and grant PICT 2017-3301 awarded by Fondo para la Investigación Científica y Tecnológica (FonCyT). Gaia Gaspar, José Alacoria, Juan Cabral, and Martín Chalela were supported by a fellowship from CONICET. Damián Mast is a member of the Carrera de Investigador Científico of CONICET. This research has made use of the http://adsabs.harvard.edu/, Cornell University arxiv.org repository, adstex (https://github.com/yymao/adstex), astropy, and the python programming language.

DATA AVAILABILITY

The code used within this article is fully available at https://github.com/Gaiana/nirdust along with some example data. The raw GNIRS data are available from the Gemini Observatory Archive (https://archive.gemini.edu) with Program ID GS-2005A-Q-38.

Footnotes

References

Allen
A.
,
Schmidt
J.
,
2015
,
J. Open Res. Softw.
,
3
,
E15

Alonso-Herrero
A.
,
Simpson
C.
,
Ward
M. J.
,
Wilson
A. S.
,
1998
,
ApJ
,
495
,
196

Alonso-Herrero
A.
et al. ,
2018
,
ApJ
,
859
,
144

Antonucci
R.
,
1993
,
ARA&A
,
31
,
473

Asmus
D.
,
2019
,
MNRAS
,
489
,
2177

Asmus
D.
,
Hönig
S. F.
,
Gandhi
P.
,
2016
,
ApJ
,
822
,
109

Astropy Collaboration
,
2013
,
A&A
,
558
,
A33

Burtscher
L.
et al. ,
2013
,
A&A
,
558
,
A149

Burtscher
L.
et al. ,
2015
,
A&A
,
578
,
A47

Combes
F.
et al. ,
2019
,
A&A
,
623
,
A79

Díaz
R. J.
et al. ,
2021
,
Bol. Asoc. Argent. Astron. La Plata Argent.
,
62
,
219

Dumont
A.
,
Seth
A. C.
,
Strader
J.
,
Greene
J. E.
,
Burtscher
L.
,
Neumayer
N.
,
2020
,
ApJ
,
888
,
19

Durré
M.
,
Mould
J.
,
2018
,
ApJ
,
867
,
149

Eikenberry
S.
et al. ,
2008
, in
McLean
I. S.
,
Casali
M. M.
, eds,
Proc. SPIE Conf. Ser. Vol. 7014, Ground-based and Airborne Instrumentation for Astronomy II
.
SPIE
,
Bellingham
, p.
70140V

Eisenhauer
F.
et al. ,
2003
, in
Iye
M.
,
Moorwood
A. F. M.
, eds,
Proc. SPIE Conf. Ser. Vol. 4841, Instrument Design and Performance for Optical/Infrared Ground-based Telescopes
.
SPIE
,
Bellingham
, p.
1548

Elias
J. H.
,
Joyce
R. R.
,
Liang
M.
,
Muller
G. P.
,
Hileman
E. A.
,
George
J. R.
,
2006
, in
McLean
I. S.
,
Iye
M.
, eds,
Proc. SPIE Conf. Ser. Vol. 6269, Ground-based and Airborne Instrumentation for Astronomy
.
SPIE
,
Bellingham
, p.
62694C

Feigenbaum
A.
,
1983
, in
Robinson
J.
,
Worth
T. L.
, eds,
Total Quality Control
.
McGraw Hill Education
,
New York, NY

Ferruit
P.
,
Mundell
C. G.
,
Nagar
N. M.
,
Emsellem
E.
,
Pécontal
E.
,
Wilson
A. S.
,
Schinnerer
E.
,
2004
,
MNRAS
,
352
,
1180

Gámez Rosas
V.
et al. ,
2022
,
Nature
,
602
,
403

Gamma
E.
,
Helm
R.
,
Johnson
R.
,
Vlissides
J.
,
1995
,
Design Patterns: Elements of Reusable Object-Oriented Software, Addison-Wesley Professional Computing Series
,
Vol. 99
.
Addison-Wesley Reading
,
Massachusetts

García-Burillo
S.
et al. ,
2021
,
A&A
,
652
,
A98

Gaspar
G.
,
Díaz
R. J.
,
Mast
D.
,
D’Ambra
A.
,
Agüero
M. P.
,
Günthardt
G.
,
2019
,
AJ
,
157
,
44

Gaspar
G.
,
Díaz
R. J.
,
Mast
D.
,
Agüero
M. P.
,
Schirmer
M.
,
Günthardt
G.
,
Schmidt
E. O.
,
2022
,
AJ
,
163
,
230

Glass
I. S.
,
Moorwood
A. F. M.
,
1985
,
MNRAS
,
214
,
429

Gomez
P. L.
et al. ,
2012
, in
American Astronomical Society Meeting Abstracts #219
. p.
413.07

Gratadour
D.
,
Clénet
Y.
,
Rouan
D.
,
Lai
O.
,
Forveille
T.
,
2003
,
A&A
,
411
,
335

GRAVITY Collaboration
,
2020
,
A&A
,
634
,
A1

Holscher
E.
,
Leifer
C.
,
Grace
B.
,
2010
, https://docs.readthedocs.org

Hönig
S. F.
,
2019
,
ApJ
,
884
,
171

Hönig
S. F.
,
Kishimoto
M.
,
2017
,
ApJ
,
838
,
L20

Hönig
S. F.
,
Kishimoto
M.
,
Antonucci
R.
,
Marconi
A.
,
Prieto
M. A.
,
Tristram
K.
,
Weigelt
G.
,
2012
,
ApJ
,
755
,
149

Hönig
S. F.
et al. ,
2013
,
ApJ
,
771
,
87

Izumi
T.
,
Wada
K.
,
Fukushige
R.
,
Hamamura
S.
,
Kohno
K.
,
2018
,
ApJ
,
867
,
48

Jazayeri
M.
,
2007
, in
Briand
L. C.
,
Wolf
A. L.
, eds,
Future of Software Engineering (FOSE’07)
.
Institute of Electrical and Electronics Engineers
,
Los Alamitos, California
, p.
199

Kishimoto
M.
,
Hönig
S. F.
,
Antonucci
R.
,
Barvainis
R.
,
Kotani
T.
,
Tristram
K. R. W.
,
Weigelt
G.
,
Levin
K.
,
2011
,
A&A
,
527
,
A121

Krolik
J. H.
,
Begelman
M. C.
,
1988
,
ApJ
,
329
,
702

Leftley
J. H.
,
Tristram
K. R. W.
,
Hönig
S. F.
,
Asmus
D.
,
Kishimoto
M.
,
Gandhi
P.
,
2021
,
ApJ
,
912
,
96

Lin
L.
,
Zou
H.
,
Kong
X.
,
Lin
X.
,
Mao
Y.
,
Cheng
F.
,
Jiang
Z.
,
Zhou
X.
,
2013
,
ApJ
,
769
,
127

López-Gonzaga
N.
,
Jaffe
W.
,
Burtscher
L.
,
Tristram
K. R. W.
,
Meisenheimer
K.
,
2014
,
A&A
,
565
,
A71

Mason
R. E.
,
2015
,
Planet. Space Sci.
,
116
,
97

Miller
J. C.
,
Maloney
C. J.
,
1963
,
Commun. ACM
,
6
,
58

Nenkova
M.
,
Ivezić
Ž.
,
Elitzur
M.
,
2002
,
ApJ
,
570
,
L9

Nenkova
M.
,
Sirocky
M. M.
,
Ivezić
Ž.
,
Elitzur
M.
,
2008
,
ApJ
,
685
,
147

Netzer
H.
,
2015
,
ARA&A
,
53
,
365

Neumayer
N.
,
Cappellari
M.
,
Reunanen
J.
,
Rix
H. W.
,
van der Werf
P. P.
,
de Zeeuw
P. T.
,
Davies
R. I.
,
2007
,
ApJ
,
671
,
1329

Ram
S. L.
et al. ,
2003
,
Dr Alan Kay on the Meaning of Object-Oriented Programming
.

Rejkuba
M.
,
2004
,
A&A
,
413
,
903

Riffel
R.
,
Pastoriza
M. G.
,
Rodríguez-Ardila
A.
,
Bonatto
C.
,
2009
,
MNRAS
,
400
,
273

Riffel
R.
et al. ,
2022
,
MNRAS
,
512
,
3906

Stalevski
M.
,
Tristram
K. R. W.
,
Asmus
D.
,
2019
,
MNRAS
,
484
,
3334

Thatte
N.
,
Quirrenbach
A.
,
Genzel
R.
,
Maiolino
R.
,
Tecza
M.
,
1997
,
ApJ
,
490
,
238

Van Rossum
G.
,
Warsaw
B.
,
Coghlan
N.
,
2001
,
Python. org
,
1565
:

Wales
D. J.
,
Doye
J. P. K.
,
1997
,
J. Phys. Chem. A
,
101
,
5111

APPENDIX A: TESTS

Here, we present three complementary tests to the ones described in Section 4.1 using synthetic spectra. In these cases, we consider different relative proportions of S/N between the target and reference spectra. The construction of the spectra follows the same steps as described before: fixed parameter values of α = 3.5, β = 8.3, and γ = −3.3, and temperature values ranging between 500 and 1600 K. The stellar population of the reference spectrum is modelled as a decreasing linear function.

In Figs A1, A2, and A3, we present the results for the scenarios where the S/N of the reference spectrum are 75, 50, and 25 per cent of the S/N of the target spectrum, respectively.

Uncertainty in the estimated parameters when varying the S/N of both the target and the reference spectra for different temperatures. The temperatures shown in the legend are the temperatures of the synthetic target spectrum in each case, the S/N of the reference spectrum corresponds to 75 per cent of the S/N of the target spectrum.
Figure A1.

Uncertainty in the estimated parameters when varying the S/N of both the target and the reference spectra for different temperatures. The temperatures shown in the legend are the temperatures of the synthetic target spectrum in each case, the S/N of the reference spectrum corresponds to 75 per cent of the S/N of the target spectrum.

Same as in Fig. A1 but the S/N of the reference spectrum corresponds to 50 per cent of the S/N of the target spectrum.
Figure A2.

Same as in Fig. A1 but the S/N of the reference spectrum corresponds to 50 per cent of the S/N of the target spectrum.

Same as in Fig. A1 but the S/N of the reference spectrum corresponds to 25 per cent of the S/N of the target spectrum.
Figure A3.

Same as in Fig. A1 but the S/N of the reference spectrum corresponds to 25 per cent of the S/N of the target spectrum.

In the three cases, the behaviour of the parameters T, α, β, and γ is similar to the case where the S/N of both spectra are equal (Fig. 1). However, a slight difference is noticeable in the low S/N regime: as the S/N fraction of the reference spectrum decreases, the point where the lowest temperatures raise their uncertainty occurs at higher target spectrum S/N.

APPENDIX B: APPLICATION EXAMPLE

In this section, we present an application example of the simplest case scenario. A more complete example can be found in the Github repository of the project.1 In this example, the fitting is performed over the nuclear spectrum of the galaxy but the reader should keep in mind that nirdust can be applied at any radius where the dust component is present.

The first step is to read the spectrum from an FITS (Flexible Image Transport System) file, providing the name of the file and the redshift of the galaxy. Then, it is necessary to follow some steps in order to ‘clean’ the spectrum of undesirable spectral features such as emission/absorption lines or residuals from the reduction process. The user may need to cut the borders of the spectrum due to higher noise or to remove spectral features such as the CO band usually present in the end of the Klong filter. All these steps can be performed with the functionalities described in the previous section. The impact of these operations on the best-fitting solution is evaluated in Section 4. These procedures must be followed for both spectra: the target and the reference. Once both spectra are pre-processed, the function fit_blackbody() can be executed to obtain the model that best describes the data. The initial values and bounds for the parameters can be provided by the user along with the amount of iterations required. Also, the contribution of the γ term can be limited through the gamma_target_fraction parameter. In this example, the bounds, the amount of iterations, and the limit on the γ term are initially given to the function. In Fig. B1, an example of the output plot is presented.

Plot showing the result of the fitting procedure as given by the method plot() of the class NirdustResults. The fitting was performed over the Flamingos-2 data set described in Section 4.2. The individual components involved in the modelling process of equation (2), labelled as α, β, and γ terms, are shown in the top panel. The residual of the fit is shown in the bottom panel and it is consistent with zero. The masked parts in the target and reference spectra are easily recognized as straight lines.
Figure B1.

Plot showing the result of the fitting procedure as given by the method plot() of the class NirdustResults. The fitting was performed over the Flamingos-2 data set described in Section 4.2. The individual components involved in the modelling process of equation (2), labelled as α, β, and γ terms, are shown in the top panel. The residual of the fit is shown in the bottom panel and it is consistent with zero. The masked parts in the target and reference spectra are easily recognized as straight lines.

>>> import nirdust as nd

>>> import matplotlib.pyplot as plt

# Read the two spectra from the FITS file

# the redshift for this galaxy is 0.00183

>>> target_spectrum = nd.read_fits(... ’nuclear_spectrum.fits’, z = 0.00183)

>>> ref_spectrum = nd.read_fits(... ’external_spectrum.fits’, z = 0.00183)

# Cut the borders of the spectrum in Angstroms

# This is user defined

>>> cut_tgt = target_spectrum.cut_edges(... 20600, 22700)

>>> cut_ref = ref_spectrum.cut_edges(... 20600, 22700)

# Find spectral lines and construct the masks

>>> lines_tgt, intervals_tgt = \ ... nd.line_spectrum(cut_tgt, noise_factor = 5.5)

>>> lines_ref, intervals_ref = \ ... nd.line_spectrum(cut_ref, noise_factor = 3.5)

# Apply the mask to the spectra

>>> mask_tgt = cut_tgt.mask_spectrum(... line_intervals = intervals_tgt)

>>> mask_ref = cut_ref.mask_spectrum(... line_intervals = intervals_ref)

# Match the spectral axis of both spectra

# so they are equally sampled:

>>> tgt, ref = nd.match_spectral_axes(... mask_tgt, mask_ref)

# Fit

>>> bounds = (... (400, 2000.0), (0, 20),... (6, 20), (-10, 10))

>>> result = nd.fit_blackbody(... tgt, ref, niter = 700, ... gamma_target_fraction = 0.01, ... bounds = bounds)

# Obtain the value of the fitted temperature

>>> T  = result.temperature.value

# Visualize the hot dust spectrum and the model

>>> axis = result.plot(show_components = True)

>>> plt.show()

B1 Quality

Quality is a user-determined characteristic and not a technical feature, which is based on the user’s experience with a product with respect to their interests and requirements, which may vary over time (Feigenbaum 1983). From developer’s perspective, software quality can be composed of qualitative and quantitative metrics. For nirdust, we determine our quality threshold by defending a coding standard, code coverage, unit tests, and continuous integration.

We chose to follow the PEP-8 standard (Van Rossum et al. 2001), which defines a coding style convention proposed by the python community, both for the standard language library and in community projects. Following PEP-8 ensures that any python programmer has a reasonable understanding of the code written and makes it easier to find collaborators and maintain the project over time. This standard is enforced using the Flake-8 tool,2 which automatically checks deviations in style and will help minimize the code errors in future versions.

In addition to style and maintainability, it is clearly obvious that we want nirdust to work properly on all computers on which it is deployed, taking into account the heterogeneity of operating systems and python versions. For this reason, we have implemented 93 unit tests to validate that the individual components of the software work correctly (Jazayeri 2007). nirdust is tested for python 3.8, 3.9, and 3.10, reaching 99 per cent of code coverage (Miller & Maloney 1963), or in other words our test-suite executes almost all the code of the project.

The source code is hosted in a public GitHub repository3 with an open-source MIT (Instituto Tecnológico de Massachusetts) License.4The committed versions of code are automatically tested with a continuous-integration workflow implemented with GitHub Actions5 (Fowler & Foemmel 2006), and the documentation is built from the repository and made public in the read-the-docs service6 (Holscher, Leifer & Grace 2010).

The project is available in the Python Package Index (PyPi)7 throught the pip installer.8 Finally, the project is registered in the Astrophysics Source Code Library (ASCL; Allen & Schmidt 2015).9

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.