# myGym: Modular Toolkit for Visuomotor Robotic Tasks

Michal Vavrecka 1, Nikita Sokovnin 1
Megi Mejdrechova1 , Gabriela Sejnova1 , Marek Otahal1

###### Abstract

We introduce a novel virtual robotic toolkit myGym, developed for
reinforcement learning (RL), intrinsic motivation, and imitation learning
tasks trained in a 3D simulator. The trained tasks can then be easily
transferred to real world robotic scenarios. The modular structure of the
simulator enables users to train and validate their algorithms on a large
number of scenarios with various robots, environments and tasks. Compared to
existing toolkits (e.g. OpenAI Gym, Roboschool) which are suitable for
classical RL, myGym is also prepared for visuomotor (combining vision &
movement) unsupervised tasks that require intrinsic motivation, i.e. the
robots are able to generate their own goals. There are also collaborative
scenarios intended for human-robot interaction. The toolkit provides
pretrained visual modules for visuomotoric tasks allowing rapid prototyping,
and, moreover, users can customize the visual submodules and retrain with
their own set of objects. In practice, the user selects the desired
environment, robot, objects, task and type of reward as simulation parameters,
and the training, visualization and testing themselves are handled
automatically. The user can thus fully focus on development of the neural
network architecture while controlling the behavior of the environment using
predefined parameters.

## Introduction

Figure 1: Examples of robotic workspaces. Top: Pepper robot by a desk (left)
and in a maze (right). Middle: A vertical plane workspace (left), and a
collaborative workspace with two KUKA arms and a humanoid (right). Bottom:
Example of image augmentation during training, where random textures from a
database are applied to the objects (left). A pretrained neural network Yolact
can be used to preprocess the images (right).

Figure 2: Overview of myGym structure. The modular system helps to prepare
desired environments and tasks to combine them compositionaly.The upper row
stands for types of modules that adds specific features to the simulation. The
bottom row standsfor the visualization of the environment after application of
specific module.

With the growing complexity of recent neural network architectures, we need to
increase the amount of training data to provide good balance between neural
network parameters and training instances. The most advanced neural networks
nowadays consists of 175 billions of parameters (Brown et al. 2020) and there
is about 300 billions of instances provided during training. In case of
textual data, there are big internet datasets available. In the area of visual
datasets, the annotation of the data is manual and very time demanding.

In the area of robotics, we require training data that are composed of visual
inputs, actual positions of the objects in the environment, the current state
of the robot actuators and also interactions within the environment. These
type of datasets are the most expensive in terms of preparation and
annotation. On the other hand, it is possible to generate an infinite amount
of such training data in a simulator without any manual effort. There are
virtual simulators that provide fully customizable environments and automatic
annotation of all objects in the scene, serving as a robust source of training
data. Since the recent simulators are also equipped with reliable physics
engines and photorealistic rendering (see bellow), the neural architectures
can be adopted in the real world tasks after successful training in the
simulator.

We introduce novel toolkit called myGym, which serves as a modular platform
for development and evaluation of cognitive architectures. Its main goal is to
let the user fully focus on development and improvement of learning
algorithms, while the simulation, task setup and visualization are
automatically taken care of.

## Related work

The possibility to train agents in a virtual environment has rapidly increased
the progress in the area of reinforcement learning. This has naturally brought
attention to further improvement of the simulation tools and, currently, one
can select among many ready-to-use libraries and toolkits. Some of them are
oriented on the gaming area (Beattie et al. 2016) (Bellemare et al. 2013),
some on navigation tasks (Savva et al. 2017)(Savva et al. 2019) or multi-agent
systems (Song et al. 2019). There is also large focus on a robot control and
robotic manipulation and one of the widely used libraries in this area is
OpenAI Gym (Brockman et al. 2016). It is a collection of environments that
includes algorithmic tasks, Atari and board games, several toy tasks, but also
2D and 3D robotic environments. There is a number of recent frameworks that
extend OpenAI Gym (or are compatible with it) and here we provide their brief
overview. For a detailed comparison, please refer to (Ivaldi et al.
2014)(Ferigo et al. 2020).

An official extension of OpenAI Gym is their Robotic Environments toolkit
(Roboschool), introduced along with a set of Baselines implementations of
state-of-the-art RL algorithms like PPO or HER (Dhariwal et al. 2017).
However, it is based on the Mujoco physics engine (Todorov, Erez, and Tassa
2012) which is not open source, similarly as the photorealistic virtual
simulator from Nvidia ISAAC. An open source alternative is PyBullet (Coumans
and Bai 2016), based on the widely used Bullet physics engine. PyBullet comes
with its own set of environments (PyBullet-gym) with implementations of some
basic agent-task scenarios. Gibson Env is also based on Bullet physics engine
and provides a highly photorealistic environment for development of real-world
perception in active agents (Xia et al. 2018).

The two frameworks that are more closely related to myGym toolkit are
PyRoboLearn (Delhaisse, Rozo, and Caldwell 2020) and Gym-Ignition (Ferigo et
al. 2020). Gym-Ignition is designed for reproducible results and enables
parallel or headless simulation, PyRoboLearn is highly modular and features a
large number of robots and tasks including sports, locomotion, manipulation or
control. However, neither of them has support for unsupervised tasks such as
visuomotor manipulation or includes pretrained tools for visual processing.
Moreover, most of the mentioned toolkits are discontinued (e.g. Roboschool) or
their modular functions are not yet fully implemented (e.g PyRoboLearn).

We would like to overcome the weak parts of the previous toolkits and extend
the usecases to the area of intrinsic motivation and imitation learning by
incorporating pretrained modules for vision. The user can then use them for
both supervised and unsupervised learning. We consider the following features
of our toolkit as the most advantageous:

  * •

interchangeable parts of the simulation

  * •

parametric definition of the environment

  * •

predefined intrinsic reward functions

  * •

predefined task types

  * •

pretrained neural networks for image preprocessing

  * •

domain randomization

  * •

visual augmentation

The highly modular structure of myGym and its support for multiple physics
engines (Mujoco, PyBullet) enables the user to easily set up a custom
environment tailored to his needs, while the adjustable randomization and
visual augmentation provide robust results that can be easily implemented in
real world robotic scenarios.

Table 1: Comparison of features provided by different frameworks that are all based on / compatible with OpenAI Gym. Toolkit | Baselines | Pretrained Vision | Augmentation | Parametric env | Parametric rewards
---|---|---|---|---|---
Roboschool | $\bullet$ | $\circ$ | $\circ$ | $\circ$ | $\circ$
PyBulletGym | $\bullet$ | $\circ$ | $\circ$ | $\circ$ | $\circ$
Gibson | $\bullet$ | $\circ$ | $\circ$ | $\circ$ | $\circ$
PyRoboLearn | $\bullet$ | $\circ$ | $\bullet$ | $\bullet$ | $\circ$
Gym-Ignition | $\circ$ | $\circ$ | $\circ$ | $\circ$ | $\circ$
myGym | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$

## Toolkit structure

The novel simulator toolkit - myGym is intended for scientists interested in
development of artificial neural network algorithms with application in
robotic tasks. The user can choose among several workspaces, robots,
manipulation tasks and levels of supervision. See Fig. 2 for details.

### Environment

Compared to previous toolkits that stem from OpenAI Gym structure, we
decomposed the environment into a set of specific modules that can be combined
together depending on user requirements. There are basic types of environments
with basic arrangement that are suitable for manipulation tasks (worskpace,
kitchen table), navigation tasks (maze, office room), curriculum learning
(mini-golf, basketball or football playgrounds with increasing complexity) or
collaborative and imitation learning (collaborative workspace or a kitchen
unit). There are also basic methods to reset and control the environment and
the connections to the physical engines that controls the simulation (Mujoco
or PyBullet). The rest of the environment is modular based on the parameters
defined by the user. Each environment in myGym can be constructed from a few
modules:

  * •

Robot - represents a simulated robot and gives an interface to control it and
get observation data

  * •

Objects - represents objects that the robot can manipulate

  * •

Camera - represents a camera with which images are taken

  * •

Task - represents a specific task to learn (reach, throw, pick-and-place,
etc.)

  * •

Reward - represents the type of reward

Each environment also includes the interactive test mode, where the user can
manipulate the objects and test the neural network effectiveness in a specific
configuration. In the following paragraphs, we describe each module in greater
detail.

### Robot and camera configuration

When the user chooses an environment with static objects, he can add the robot
to the environment just by adding the parameter with specified name, position,
orientation and initial configuration of the robot. The toolkit provides a set
of industrial robots (Kuka IIWA, Franka Emika, Universal robots), research
robots (Reachy arm, Gummi arm, Jaco arm), dual arm robot Yumi and also
humanoid robot for navigation tasks (Pepper robot). The user can also connect
a gripper that is suitable for the selected task as there are non-grasping
grippers (magnetic gripper, Festo universal gripper), non-tactile grippers
(Brunel Hand, Schunk two finger gripper) and also tactile grippers (RightHand
Takktile and Barret Hand). You can see the examples of robots and grippers in
Fig. 3. There are many combinations of robots and grippers, but the toolkit
will not allow user to create an impossible configuration (e.g. to connect a
tactile gripper to a humanoid robot). After a successful initialization of the
the robot in the environment, the action space is automatically set according
to the number of joints and tactile sensors. The initialization of the cameras
is similar to the initialization of the robot. In case that the robot is not
equipped with visual sensors, the user can specify the amount of cameras,
their positions and other properties in the configuration file. Otherwise, the
cameras are fixed to the robot body (e.g. Pepper or Realsense at Franka
Emika). The user then specify the active cameras that will be used during
training and testing.

Figure 3: Examples of robotic arms and exchangeable grippers included in
myGym. The user can easily add any robot to the environment just by specifying
its name and optionally position, orientation and initial configuration.

### Manipulable objects configuration

As most of the visuomotor tasks require object manipulation, the toolkit comes
with several sets of objects such as household items, geometric primitives or
industrial tools. These objects have been chosen as suitable for manipulation
tasks and the available visual modules were pretrained on them. The user can
also easily add custom objects to the toolkit and retrain the visual modules
if needed.

The user defines a desired subset of all available objects during the
configuration stage simply by stating their names. The user can further
specify position in the scene, where each object will be initialized at the
beginning of a training episode, by providing its spatial coordinates. It is
also possible to randomize the initialization of objects by specifying an area
instead of a position, varying the number of initialized objects and
performing random choice of objects from the given subset.

### Task specification

Compared to previous toolkits (e.g. Roboschool, PyRoboLearn, PyBullet-gym), we
simplified the process of task specification and reward calculation. There is
reward and task specified in each environment separately in the previous
toolkits and the user has to manually modify the environment to a novel task.
We improved this part of the simulator as there are generic types presented
that should be easily adopted to any environment. We realize that most of the
manipulation and navigation tasks share a common assumption about the reward
calculation. At the beginning, the user has to identify objects in the
environment that are crucial for reward calculation just by assigning them a
target and a goal identifier. As most of the tasks (reach, push, throw, pick,
place, navigate) requires just to calculate the distance between target and
goal object (the object should be also part of robot body or position on the
objects surface), this process can be easily set just by few parameters. As
the user defines the type of the task and crucial objects, it is also
necessary to define what source of information is used to calculate the
position of objects. This process is described in the next section.

### Reward calculation

The specific training reward depends both on the selected type of task (i.e.
push, pick and place, reach) and type of learning (supervised RL or
unsupervised/intrinsic learning). The common objective for all rewards in our
toolkit is minimization of distance between two instances - usually between
two objects (pick and place, push tasks) or between the agent and an object
(reaching task).

In case of fully supervised learning, the reward is calculated based on ground
truth positions of the agent and objects, as retrieved from the simulator. In
case of the intrinsic motivation, the reward is calculated based on image
input only - in such case, we provide pretrained neural submodules that can
convert the image into a representation that can serve for reward calculation.
For example, the distance between an initial and goal position of a single
object can be calculated from latent representations encoded by a pretrained
variational autoencoder. In case of two or more objects in the scene, the
distance can be retrieved from centroids inferred by a pretrained YOLACT
network.

We also considered that the distance itself can be calculated in many
different ways and thus influence the learning results. Therefore, we leave it
up to the user to experiment and choose among several most commonly used
distance functions (e.g. Euclidean, Manhattan or Mahalanobis distance).

In practice, the user only needs to specify the task, the type of distance to
minimize and whether to calculate the reward from ground truth or pre-
processed image. It is also possible to exchange the pretrained visual
networks for a custom architecture or to design an own, shaped reward to train
on.

### Augmentation

Data augmentation is a key component for simulation-to-reality (Sim2Real)
transfer. In order to successfully apply reinforcement learning to real
robotic tasks, the systems based on computer vision should be sufficiently
accurate. To handle gaps between synthetic data and real-world images with an
endless variety of material appearance, colors and lighting conditions, we use
the simple and fast technique - domain randomization (Tobin et al. 2017). The
main purpose of this approach is to provide enough variability of simulated
data at training time such that the model will be able to successfully
generalize to real-world data during the test time.

The whole environment can be randomized at several levels: random
colors/textures applied to all objects, random camera position and
orientation, variable light position, type and intensity and several
postprocessing filters. There is example of the virtual environment
augmentation can be seen in Fig. 1.

Due to modularity, each environment based on PyBullet in myGym can be easily
randomized using a special wrapper. Randomization parameters can be managed by
a configuration file. Each parameter is controlled by a corresponding
randomizer component.

  * •

Light Randomizer - randomizes the parameters of simulation lighting such as
direction, color, distance, ambient, diffuse and specular. Desirable elements
of light configuration and their ranges can be also specified.

  * •

Camera Randomizer - randomizes the position and orientation of the camera
using for images rendering.

  * •

Texture Randomizer - randomizes objects texture using a specified path to
folder with images.

  * •

Color Randomizer - randomizes color attributes of objects in the scene such as
RGBA color and specular color of the object. It can be applied at the same
time as Texture Randomizer.

  * •

Joint Randomizer - add random position to the robotic actuators

  * •

Postprocessing Randomizer - application of special filters to the camera image

### Pretrained modules

The toolkit allows users to start building their applications in the shortest
possible time by offering three already pre-trained modules helpful in
visuomotor tasks.

YOLACT (Bolya et al. 2019) network performs fundamental image segmentation and
object recognition. It is able to identify all robots and objects used in the
toolkit and to determine their position in the scene from the camera inputs in
real-time. DOPE (Tremblay et al. 2018) extends YOLACT capabilities allowing to
obtain 6DoF information about objects - position and orientation in 3D space.
VAE (Kingma and Welling 2013) extracts important information about the scene
by encoding the camera inputs from pixel to latent space. All three vision
modules can be used for intrinsic motivation reward calculation, the VAE is
even trained an unsupervised manner.

In case the user wants to train a customized version of any of the vision
modules, the toolkit provides an option to generate a training dataset. The
user configures the environment, robots and objects in a very similar way as
before training, and chooses what type and size of dataset to generate. At
this stage, it is very helpful to use randomization that is also easily
configurable. The training itself is further module-specific and the toolkit
provides guidance on how to perform it.

### Training and testing

The purpose of the RL training environment is to apply your model and
continuously get the rewards from your defined utility function. A common
entry point is a main script where you can specify the environment, robots, RL
algorithm (e.g. baselines) or optional pretrained vision module. After
training, it is possible to evaluate your model by running main script with
test parameters, which will give you average final score. There is also an
option to record your (visualized) run, which is useful for further sharing
and analysis. There is a code for baseline training provided in form of
oneliners with simulation parameters.

The Gym comes with pretrained submodules for vision, allowing to quickly start
developing your RL model in a real life tasks. Without the actual need to
train your models for vision, which is demanding on HW resources and a lengthy
procedure. There are 2 main modes for vision - for unsupervised learning
(using variational autoencoders), and supervised learning (YOLACT, DOPE).

There are also functions for the visualization of the results after the
training. As the researchers are interested in comparison of their own
algorithms with respect to the state of the art, we provide evaluation
protocol to directly compare the performance. You can see the example in Fig.
4

Figure 4: Evaluation protocol for direct comparison of the neural networks
performance in the specific environment. The user can compare the performance
of own neural network with the state of the art methods.

### Sim2real - Transferring to physical sensors and robots with ROS

As for many practical use-cases, interfacing with various hardware (sensors,
cameras, robotic arms, etc.) is needed. The industry golden-standard is Robot
Operating System (ROS), which is a set of drivers and various libraries for
designing robotic applications. Therefore our toolkit also has ROS2

The most common use-cases are twofolds:

  * •

cameras (sensors): You can seamlessly switch (or combine) virtual cameras
provided by the simulator environment, and any physical cameras.

  * •

robots (actuators): Since the common pitfalls of real HW equipment (expensive,
slow, wear), users can experiment and train their architectures in the
simulator. The provided models of robots are very precise, so the main idea is
that it is possible to switch the entire pipe-line from simulator to your real
environment.

#### Connecting pretrained network to real cameras and robots using ROS

The concept of ROS is a so called Node (a physical HW or an algorithm), and a
running application consists of a network of nodes that communicate (sending
messages).

As a part of the builtin functionality in our Gym, we provide a ROS2 node for
the Intel RealSense camera

We have extended our Vision submodule to be used as a ROS node. This allows
the user to train the vision system in simulator and then switch to ROS and
use the trained model

## Conclusion

We present a new virtual simulator toolkit myGym, which serves as a modular
extension of the widely used OpenAI Gym simulator. Its main strengths lie in
simple parametric definition of a custom composite environment, predefined
task types, domain randomization and full visual augmentation. Furthermore, it
is the first simulator toolkit with support for intrinsic motivation - this is
achieved by predefined intrinsic reward functions and pretrained neural
networks that can be used for image preprocessing and reward calculation
instead of the ground truth.

In our following work, we would like to further develop full support for all
the mentioned features in both of the physics engines - currently, the closed
nature of Mujoco physics engine does not enable some methods like visual
augmentation or full randomization.

We believe that wide variety of users will find the presented tool useful for
their research. We have developed the toolkit so that the beginners can train
neural networks in various environments while all necessary files are
provided, and experts can modify it according to their needs.

## References

  * Beattie et al. (2016) Beattie, C.; Leibo, J. Z.; Teplyashin, D.; Ward, T.; Wainwright, M.; Küttler, H.; Lefrancq, A.; Green, S.; Valdés, V.; Sadik, A.; et al. 2016\.  Deepmind lab.  _arXiv preprint arXiv:1612.03801_ .
  * Bellemare et al. (2013) Bellemare, M. G.; Naddaf, Y.; Veness, J.; and Bowling, M. 2013.  The arcade learning environment: An evaluation platform for general agents.  _Journal of Artificial Intelligence Research_ 47: 253–279.
  * Bolya et al. (2019) Bolya, D.; Zhou, C.; Xiao, F.; and Lee, Y. J. 2019.  YOLACT: real-time instance segmentation.  In _Proceedings of the IEEE International Conference on Computer Vision_ , 9157–9166.
  * Brockman et al. (2016) Brockman, G.; Cheung, V.; Pettersson, L.; Schneider, J.; Schulman, J.; Tang, J.; and Zaremba, W. 2016.  Openai gym.  _arXiv preprint arXiv:1606.01540_ .
  * Brown et al. (2020) Brown, T. B.; Mann, B.; Ryder, N.; Subbiah, M.; Kaplan, J.; Dhariwal, P.; Neelakantan, A.; Shyam, P.; Sastry, G.; Askell, A.; et al. 2020.  Language models are few-shot learners.  _arXiv preprint arXiv:2005.14165_ .
  * Coumans and Bai (2016) Coumans, E.; and Bai, Y. 2016.  Pybullet, a python module for physics simulation for games, robotics and machine learning .
  * Delhaisse, Rozo, and Caldwell (2020) Delhaisse, B.; Rozo, L.; and Caldwell, D. G. 2020.  Pyrobolearn: A Python framework for robot learning practitioners.  In _Conference on Robot Learning_ , 1348–1358.
  * Dhariwal et al. (2017) Dhariwal, P.; Hesse, C.; Klimov, O.; Nichol, A.; Plappert, M.; Radford, A.; Schulman, J.; Sidor, S.; Wu, Y.; and Zhokhov, P. 2017.  OpenAI Baselines.  https://github.com/openai/baselines.
  * Ferigo et al. (2020) Ferigo, D.; Traversaro, S.; Metta, G.; and Pucci, D. 2020.  Gym-Ignition: Reproducible Robotic Simulations for Reinforcement Learning.  In _2020 IEEE/SICE International Symposium on System Integration (SII)_ , 885–890. IEEE.
  * Ivaldi et al. (2014) Ivaldi, S.; Peters, J.; Padois, V.; and Nori, F. 2014.  Tools for simulating humanoid robot dynamics: a survey based on user feedback.  In _2014 IEEE-RAS International Conference on Humanoid Robots_ , 842–849. IEEE.
  * Kingma and Welling (2013) Kingma, D. P.; and Welling, M. 2013.  Auto-encoding variational bayes.  _arXiv preprint arXiv:1312.6114_ .
  * Savva et al. (2017) Savva, M.; Chang, A. X.; Dosovitskiy, A.; Funkhouser, T.; and Koltun, V. 2017.  MINOS: Multimodal indoor simulator for navigation in complex environments.  _arXiv preprint arXiv:1712.03931_ .
  * Savva et al. (2019) Savva, M.; Kadian, A.; Maksymets, O.; Zhao, Y.; Wijmans, E.; Jain, B.; Straub, J.; Liu, J.; Koltun, V.; Malik, J.; et al. 2019.  Habitat: A platform for embodied ai research.  In _Proceedings of the IEEE International Conference on Computer Vision_ , 9339–9347.
  * Song et al. (2019) Song, Y.; Wang, J.; Lukasiewicz, T.; Xu, Z.; Xu, M.; Ding, Z.; and Wu, L. 2019.  Arena: A general evaluation platform and building toolkit for multi-agent intelligence.  _arXiv preprint arXiv:1905.08085_ .
  * Tobin et al. (2017) Tobin, J.; Fong, R.; Ray, A.; Schneider, J.; Zaremba, W.; and Abbeel, P. 2017.  Domain Randomization for Transferring Deep Neural Networks from Simulation to the Real World.
  * Todorov, Erez, and Tassa (2012) Todorov, E.; Erez, T.; and Tassa, Y. 2012.  Mujoco: A physics engine for model-based control.  In _2012 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 5026–5033. IEEE.
  * Tremblay et al. (2018) Tremblay, J.; To, T.; Sundaralingam, B.; Xiang, Y.; Fox, D.; and Birchfield, S. 2018\.  Deep object pose estimation for semantic robotic grasping of household objects.  _arXiv preprint arXiv:1809.10790_ .
  * Xia et al. (2018) Xia, F.; Zamir, A. R.; He, Z.; Sax, A.; Malik, J.; and Savarese, S. 2018.  Gibson env: Real-world perception for embodied agents.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 9068–9079.


# NIRDust: Probing Hot Dust Emission around Type 2 AGN using K-band Spectra

Gaia Gaspar 1,5 Martín Chalela 1,2,5 Juan Cabral 2,3,5 José Alacoria 4 Damián
Mast 1,5 and Rubén Díaz 1,6
1Observatorio Astronómico de Córdoba, Universidad Nacional de Córdoba, Laprida
854, X5000BGR, Córdoba, Argentina (OAC–UNC).
2Instituto De Astronomía Teórica y Experimental, Córdoba, Argentina
(IATE–CONICET).
3Gerencia De Vinculación Tecnológica Comisión Nacional de Actividades
Espaciales (CONAE), Falda del Cañete, Córdoba, Argentina.
4Instituto de Ciencias Astrómicas, de la Tierra y el Espacio, Av. España Sur
1512, San Juan, Argentina (ICATE–CONICET–UNSJ).
5Consejo de Investigaciones Científicas y Técnicas de la República Argentina,
Buenos Aires, Argentina (CONICET).
6Gemini Observatory, NSF’s NOIRLab, USA E-mail<EMAIL_ADDRESS>0000-0001-9293-4449 0000-0002-9143-7974 0000-0002-7351-0680
0000-0001-9726-8991 0000-0003-0469-3193 0000-0001-9716-5335

(Accepted 21 Dec 2023)

###### Abstract

Hot dust in the proximity of AGNs strongly emits in the Near Infrared
producing a red excess that, in type 2 sources, can be modeled to measure its
temperature. In the era of high spatial-resolution multi-wavelength data,
mapping the hot dust around Supermassive Black Holes is important for the
efforts to achieve a complete picture of the dust 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 ($\lambda_{c}$ = 2.2 $\mu$m). To achieve this, we have
developed NIRDust, a Python package for modeling K-band spectra, estimate the
dust temperature and characterize 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 $\sim$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 GNIRS spectra
using two different approaches.

###### keywords:

ISM: dust – Infrared: ISM – Galaxies: nuclei – Galaxies: active – Techniques:
spectroscopic – Methods: data analysis –

††pubyear: 2023††pagerange: NIRDust: Probing Hot Dust Emission around Type 2
AGN using K-band Spectra–B.1

##  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 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 it could play an important role in the accretion mechanism of
the 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 et al., 2002, 2008).
Later, the discovery of a polar component observed with MIR interferometry
(e.g. Hönig et al., 2012, 2013; López-Gonzaga et al., 2014; Asmus et al.,
2016; Stalevski et al., 2019; Asmus, 2019) resulted in a new model that
incorporates a polar component as an outflowing wind and a thin annulus of
high temperature dust ($\sim$ 1500 K) observed in 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. Izumi et al., 2018; Alonso-Herrero 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 composition 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 with 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 a SMBH part of the
emission is detected in the near-infrared and, in particular, is prominent in
the K-band (2-2.4$\mu$m). Several authors have reported measurements of the
temperature of hot dust in AGNs using JHK colors (e.g. Glass & Moorwood, 1985;
Alonso-Herrero et al., 1998), NIR interferometry (e.g. Kishimoto et al., 2011;
Leftley et al., 2021; Gravity Collaboration et al., 2020; 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” apertures and found a mean of 887
$\pm$ 68 K for type 2 AGNs and 1292 $\pm$ 46 K for type 1 AGNs. Regarding the
extension and temperature distribution of this hot dust, Gaspar et al. (2019,
2022) found T $\sim 1000$ K dust in scales of tens of pc and almost constant
in temperature in two Seyfert 2 nuclei. Gámez Rosas et al. (2022) measured
dust temperatures using MIR interferometry with MATISSE/ESO/VLTI instrument in
several sub-pc regions around the SMBH of NGC 1068, the archetypical Seyfert 2
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 hot dust
emission extended outside the unresolved core and stated that the emission
inside $\sim$ 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 on 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 MIRI (5-28 $\mu$m) and NIRSpec (0.6
- 5 $\mu$m) instruments on board the James Webb Space Telescope 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 (API) facilitates seamless
integration with Astropy (Astropy Collaboration et al., 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 et al., 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 is not only detected in
the nuclear spectrum but it 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
modeled 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 $\nu$, and parameterized by its temperature, $T$, and an amplitude,
$A$, as:

$B_{T}(\nu)=A\frac{2h\nu^{3}}{c^{2}}\frac{1}{e^{h\nu/kT}-1}$ (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 modeled 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):

$\displaystyle S_{t}$ $\displaystyle=\alpha
S_{s}+\beta^{\prime}B_{T}+\gamma^{\prime}$ (2) $\displaystyle\beta^{\prime}$
$\displaystyle=10^{\,\beta}$ $\displaystyle\gamma^{\prime}$
$\displaystyle=10^{\,\gamma}$

Here, $S_{t}$ is the Target Spectrum, $S_{s}$ is the Reference Spectrum that
accounts for the stellar emission, $B_{T}$ is the blackbody function of
temperature T, and $\gamma^{\prime}$ is a constant emission that represents
the sky background (it is recommended that the sky background is previously
subtracted to make $\gamma^{\prime}$ the closer to zero as possible). $\alpha$
and $\beta^{\prime}$ are scale factors of the stellar and blackbody components
respectively. The reformulation of $\beta^{\prime}$ and $\gamma^{\prime}$ 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\\{T,\alpha,\beta,\gamma\\}$.

In order to model the hot dust emission with Eq. 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\%$ in the NIR due
to the fast changes in atmospheric conditions that occur at this 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 $\beta$ 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 $S_{s}$ is
the same as the Target Spectrum stellar component, in other words, that the
stellar population that emitted $S_{s}$ is the same stellar population which
inhabits the region that emitted $S_{t}$.

$S_{s}$ is , of course, not possible to isolate from the Target Spectrum
without knowing the dust emission, but there are some hypothesis 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 emitted 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) were 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 is available, analyzing 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 Sec. 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, $S_{s}$ 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 analyzing 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, that was
applied by Burtscher et al. (2015), and the third scenario, which remains
hypothetical, require further testing to determine their suitability for
modeling the stellar population of AGNs using NIRDust.

##  3 NIRDust

The entire experiment of this work are 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 centered 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 preprocessing 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 in estimating
and characterizing the dust temperature, this algorithm is explained in detail
in the Subsection 3.1.

To guarantee accurate results, NIRDust maintains a strict adherence to the
PEP-8 coding standard (Van Rossum et al., 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%
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\\{T,\alpha,\beta,\gamma\\}$. 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:

$\mathcal{L}(D|\boldsymbol{\theta})=\prod^{n}_{i=1}f(\nu_{i}|\boldsymbol{\theta})$
(3)

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

$\mathcal{L}(D|\boldsymbol{\theta})=\frac{1}{(2\pi\sigma^{2})^{n/2}}\exp\left(-\frac{1}{2}\sum^{n}_{i=1}\frac{(D(\nu_{i})-S_{T}(\nu_{i}|\boldsymbol{\theta}))^{2}}{\sigma^{2}}\right)$
(4)

where $\sigma$ 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 $\sigma$ 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 (CenA), 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 asses
their impact.

  * •

The signal-to-noise (S/N) ratio 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 modeled the stellar population of the
Reference Spectrum as a decreasing linear function. The Target spectrum was
then modeled according to equation 2 where we adopted fixed parameter values
of $\alpha=3.5$, $\beta=8.3$, and $\gamma=-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_blabckbody() and found that the results
remain unaltered. In all cases a different gaussian noise (different random
seed) is added to each spectra to obtain the desired S/N ratio. The spectral
resolution of both spectra is set at 3.5 Å/pix and they consist of 541
spectral points. Throughout all this work, we’ve 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, $\Delta 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 $\alpha$, $\beta$, and $\gamma$
parameters are respectively presented. At high S/N ($\gtrsim 150$) the
temperature is accurately estimated for temperatures below 900 K, but for
higher temperatures a systematic and constant overestimation of $\lesssim$ 50
K is obtained. However, with decreasing S/N (20-50), lower temperatures show
an increase in the overestimation reaching $\Delta T\sim$ 100 K. As can be
seen in the other three panels, $\alpha$ and $\beta$ follow the same trend,
remaining constant for S/N $\gtrsim$ 150 and deviating for lower S/N in the
cases of the lower temperatures ($\lesssim$ 900 K). The $\gamma$ parameter
remains constant along the entire T and S/N ranges and represents a minor
contribution to the total flux ($<1\%$).

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$\gtrsim$ 150 is required. It is important to remark that in
all cases with S/N $\geq$ 20 the uncertainty correspond 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 Subsec. 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.

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.

In Fig. 2 we present results regarding the application of masks used to remove
spectral features or high noise regions. The test consisted on 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.

Four cases where considered:

  * •

"4 intervals" mask: a mask consisting of 4 equi-distant 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 $\mu$m (top right panel).

  * •

"Red interval" mask: a mask consisting of one interval in the red part of the
spectrum at 2.24 $\mu$m (bottom left panel).

  * •

"Central interval" mask: and a mask consisting in one interval in the central
part of the spectral range at 2.15 $\mu$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%.

In all cases the uncertainty in the estimated temperature is mainly constant
until a 30% of mask percentage is reached. Taking into account that the S/N of
the spectra involved in this test is of 800 it can be considered that this
initial overestimation of the temperature is attributable to the S/N effect as
analyzed in the previous test. As the masking percentage increases, however,
the uncertainty in the obtained temperature changes in all the cases except
the 4 intervals mask (top left panel). In the other three cases, the
uncertainty increases and even changes sign. This is expected since the 4
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% of the points.

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: 4 equi-
distant and equi-sized intervals distributed along the spectral range (top
left panel); one interval in the blue part of the spectrum at 2.07 $\mu$m (top
right panel); one interval in the red part of the spectrum at 2.24 $\mu$m
(bottom left panel); and one interval in the central part of the spectral
range at 2.15 $\mu$m (bottom right panel). In all cases the percentage of
removed points is the sum of the masked intervals. The temperatures showed in
the legend are the temperatures of the blackbody model of the synthetic Target
Spectrum in each case.

###  4.2 Real spectra: Centaurus A

We have used three different sets of spectra of NGC 5128 (Centaurus 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\pm 1)$ K for the hot dust in Cen A nucleus
measured in a 1" diameter aperture. Their spectra were selected from the ESO
SINFONI (Eisenhauer et al., 2003) database and correspond to IFU observations
in K-band with seeing of 0.6" taken in March-April, 2005 (Neumayer et al.,
2007). NGC 5128 was chosen as the primary subject of our comparative study due
to the availability of two extra suitable datasets 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: in one hand an off-nuclear
spectrum of the same galaxy as described in Sec.2.1 (Reference 1), and in 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 datasets:

  * •

$K_{long}$ longslit spectra taken with Flamingos-2 (Eikenberry et al., 2008;
Gomez et al., 2012) at Gemini South with PA =151° (semimajor-axis direction),
on February 21, 2021. The percentage of masking implemented is of 30% and the
mask consist of 8 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 4 intervals and one blue interval) the uncertainty is
expected to be dominated by the S/N of the spectra (see subsec. 4.1).

  * •

A $K_{long}$ 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% and the mask consist 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 dataset,
shown in the bottom panel of Fig. 3, a distinct approach was necessary due to
the instrument’s limited field of view (FOV). Instead of employing annular
apertures, circular apertures were extracted from different regions
surrounding the nucleus at a fixed distance of 2.25". 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’s 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.

The characteristics of the four extracted spectra, along with the SINFONI
spectra from Burtscher et al. (2013), are summarized in Table
LABEL:table:data. The table includes information on the masking percentage,
the temperature obtained using NIRDust, and the expected uncertainty
calculated from the S/N 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 may not be the same for all the references. The Target Spectrum presents the same masking percentage as the Reference in each case. Spectrum | data type | S/N | aperture | radius of | masking | temperature | expected
---|---|---|---|---|---|---|---
|  |  | diameter ["] | extraction ["] | percentaje | [K] |  $\Delta$T
|  |  | Flamingos-2 |  |  |  |
Target | longslit | 95 | 1 | 0 | - | - | -
Reference 1 | longslit | 72 | 1.08 | 3.34 | 30 | 662 | 25 - 85
Reference 2 (HD 176617) | IFU | 40 | point source | - | 30 | 667 | 90 - 110
|  |  | GNIRS |  |  |  |
Target | IFU | 87 | 1.05 | - | - | - | -
Reference 1 | IFU | 74 | 1.05 | 2.25 | 10 | 697 | 5 - 75
Reference 2 (HD 176617) | IFU | 40 | point source | - | 0 | 607 | 90-110
|  |  | SINFONI |  |  |  |
Target | IFU | 205 | 1. | 0 | - | - | -
Reference (HD 176617) | IFU | 40 | point source | - | 30 | 649 | 100 - 110

Figure 3: Top panel: Flamingos-2 spectra extracted at different radii (r) at
both sides of the slit (SE and NW directions). In red and black 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 black line
is a guide to compare the Target and Reference spectra slopes.

Figure 4: NIRDust fitting results for the 5 combinations of spectra described
in Table LABEL:table:data. 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 5 cases the plots are the actual output of NIRDust, to see
the obtained temperatures reference to Table LABEL:table:data

In Fig. 4 we present the NIRDust fitting results for the 5 combinations of
spectra described in Table LABEL:table:data. 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 Sec.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 LABEL:table:data. Sec.4.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
$\sim$ 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 is worth
noting that all dust temperatures measured by NIRDust and by Burtscher et al.
(2015) fall in the range of relatively low temperature that have been
previously reported for Seyfert 2 nucleii such as the one hosted NGC 5128
(Riffel et al., 2009).

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 Sec. 2.

##  5 Summary and future work

In this work we have presented an analysis on 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 on 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}=\\{T,\alpha,\beta,\gamma\\}$ 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 which value depends primarily on the S/N of the
Target Spectrum and can be accurately characterized for each dataset. In most
cases the error of the temperatures is below 100 K for S/N $\geq$ 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 (Centaurus A) using three
independent data sets from Flamingos-2 (662 K and 667 K), GNIRS (697 K 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 $\sim$ 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 disk 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 it
significantly contributed to our analysis; to 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) and 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 Foncyt. G.G, J.A., J.B.C and M.Ch. are supported by a
fellowship from CONICET. D.M. 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 xxx.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.

## References

  * Allen & Schmidt (2015) Allen A., Schmidt J., 2015, Journal of Open Research Software, 3, E15
  * Alonso-Herrero et al. (1998) Alonso-Herrero A., Simpson C., Ward M. J., Wilson A. S., 1998, ApJ, 495, 196
  * Alonso-Herrero et al. (2018) Alonso-Herrero A., et al., 2018, ApJ, 859, 144
  * Antonucci (1993) Antonucci R., 1993, ARA&A, 31, 473
  * Asmus (2019) Asmus D., 2019, MNRAS, 489, 2177
  * Asmus et al. (2016) Asmus D., Hönig S. F., Gandhi P., 2016, ApJ, 822, 109
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Burtscher et al. (2013) Burtscher L., et al., 2013, A&A, 558, A149
  * Burtscher et al. (2015) Burtscher L., et al., 2015, A&A, 578, A47
  * Combes et al. (2019) Combes F., et al., 2019, A&A, 623, A79
  * Díaz et al. (2021) Díaz R. J., et al., 2021, Boletin de la Asociacion Argentina de Astronomia La Plata Argentina, 62, 219
  * Dumont et al. (2020) Dumont A., Seth A. C., Strader J., Greene J. E., Burtscher L., Neumayer N., 2020, ApJ, 888, 19
  * Durré & Mould (2018) Durré M., Mould J., 2018, ApJ, 867, 149
  * Eikenberry et al. (2008) Eikenberry S., et al., 2008, in McLean I. S., Casali M. M., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 7014, Ground-based and Airborne Instrumentation for Astronomy II. p. 70140V, doi:10.1117/12.788326
  * Eisenhauer et al. (2003) Eisenhauer F., et al., 2003, in Iye M., Moorwood A. F. M., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 4841, Instrument Design and Performance for Optical/Infrared Ground-based Telescopes. pp 1548–1561 (arXiv:astro-ph/0306191), doi:10.1117/12.459468
  * Elias et al. (2006) 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, SPIE Conference Series Vol. 6269, Ground-based and Airborne Instrumentation for Astronomy. p. 62694C, doi:10.1117/12.671817
  * Feigenbaum (1983) Feigenbaum A., 1983, Total quality control
  * Ferruit et al. (2004) Ferruit P., Mundell C. G., Nagar N. M., Emsellem E., Pécontal E., Wilson A. S., Schinnerer E., 2004, MNRAS, 352, 1180
  * Fowler & Foemmel (2006) Fowler M., Foemmel M., 2006, Continuous integration
  * Gámez Rosas et al. (2022) Gámez Rosas V., et al., 2022, Nature, 602, 403
  * Gamma et al. (1995) Gamma E., Helm R., Johnson R., Vlissides J., Patterns D., 1995, Elements of reusable object-oriented software. Addison-Wesley Professional Computing Series Vol. 99, Addison-Wesley Reading, Massachusetts
  * García-Burillo et al. (2021) García-Burillo S., et al., 2021, A&A, 652, A98
  * Gaspar et al. (2019) Gaspar G., Díaz R. J., Mast D., D’Ambra A., Agüero M. P., Günthardt G., 2019, AJ, 157, 44
  * Gaspar et al. (2022) 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 & Moorwood (1985) Glass I. S., Moorwood A. F. M., 1985, MNRAS, 214, 429
  * Gomez et al. (2012) Gomez P. L., et al., 2012, in American Astronomical Society Meeting Abstracts #219. p. 413.07
  * Gratadour et al. (2003) Gratadour D., Clénet Y., Rouan D., Lai O., Forveille T., 2003, A&A, 411, 335
  * Gravity Collaboration et al. (2020) Gravity Collaboration et al., 2020, A&A, 634, A1
  * Holscher et al. (2010) Holscher E., Leifer C., Grace B., 2010, URL< https://docs. readthedocs. org
  * Hönig (2019) Hönig S. F., 2019, ApJ, 884, 171
  * Hönig & Kishimoto (2017) Hönig S. F., Kishimoto M., 2017, ApJ, 838, L20
  * Hönig et al. (2012) Hönig S. F., Kishimoto M., Antonucci R., Marconi A., Prieto M. A., Tristram K., Weigelt G., 2012, ApJ, 755, 149
  * Hönig et al. (2013) Hönig S. F., et al., 2013, ApJ, 771, 87
  * Izumi et al. (2018) Izumi T., Wada K., Fukushige R., Hamamura S., Kohno K., 2018, ApJ, 867, 48
  * Jazayeri (2007) Jazayeri M., 2007, in Future of Software Engineering (FOSE ’07). pp 199–213, doi:10.1109/FOSE.2007.26
  * Kishimoto et al. (2011) 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 & Begelman (1988) Krolik J. H., Begelman M. C., 1988, ApJ, 329, 702
  * Leftley et al. (2021) Leftley J. H., Tristram K. R. W., Hönig S. F., Asmus D., Kishimoto M., Gandhi P., 2021, ApJ, 912, 96
  * Lin et al. (2013) Lin L., Zou H., Kong X., Lin X., Mao Y., Cheng F., Jiang Z., Zhou X., 2013, ApJ, 769, 127
  * López-Gonzaga et al. (2014) López-Gonzaga N., Jaffe W., Burtscher L., Tristram K. R. W., Meisenheimer K., 2014, A&A, 565, A71
  * Mason (2015) Mason R. E., 2015, Planet. Space Sci., 116, 97
  * Miller & Maloney (1963) Miller J. C., Maloney C. J., 1963, Communications of the ACM, 6, 58
  * Nenkova et al. (2002) Nenkova M., Ivezić Ž., Elitzur M., 2002, ApJ, 570, L9
  * Nenkova et al. (2008) Nenkova M., Sirocky M. M., Ivezić Ž., Elitzur M., 2008, ApJ, 685, 147
  * Netzer (2015) Netzer H., 2015, ARA&A, 53, 365
  * Neumayer et al. (2007) 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 et al. (2003) Ram S. L., et al., 2003, Dr. Alan Kay on the Meaning of Object-Oriented Programming
  * Rejkuba (2004) Rejkuba M., 2004, A&A, 413, 903
  * Riffel et al. (2009) Riffel R., Pastoriza M. G., Rodríguez-Ardila A., Bonatto C., 2009, MNRAS, 400, 273
  * Riffel et al. (2022) Riffel R., et al., 2022, MNRAS, 512, 3906
  * Stalevski et al. (2019) Stalevski M., Tristram K. R. W., Asmus D., 2019, MNRAS, 484, 3334
  * Thatte et al. (1997) Thatte N., Quirrenbach A., Genzel R., Maiolino R., Tecza M., 1997, ApJ, 490, 238
  * Van Rossum et al. (2001) Van Rossum G., Warsaw B., Coghlan N., 2001, Python. org, 1565
  * Wales & Doye (1997) Wales D. J., Doye J. P. K., 1997, Journal of Physical Chemistry A, 101, 5111

##  Appendix A Tests

Here we present three complementary tests to the ones described in Subsec. 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 follow the same steps as described before: fixed parameter
values of $\alpha=3.5$, $\beta=8.3$, and $\gamma=-3.3$, and temperatures
values ranging between 500 and 1600 K. The stellar population of the Reference
spectrum is modeled as a decreasing linear function.

In Figures 5, 6 and 7 we present the results for the scenarios where the S/N
of the Reference Spectrum is $75\%$, $50\%$ and $25\%$ of the S/N of the
Target Spectrum, respectively.

In the three cases the behaviour of the parameters T, $\alpha$, $\beta$, and
$\gamma$ is similar to the case where the S/N of both spectra is 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.

Figure 5: Uncertainty in the estimated parameters when varying the S/N of both
the Target and the Reference spectra for different temperatures. The
temperatures showed in the legend are the temperatures of the synthetic Target
Spectrum in each case, the S/N of the Reference Spectrum corresponds to 75% of
the S/N of the Target Spectrum. Figure 6:  Same as in Fig. 5 but the S/N of
the Reference Spectrum corresponds to 50% of the S/N of the Target Spectrum.
Figure 7: Same as in Fig. 5 but the S/N of the Reference Spectrum corresponds
to 25% of the S/N of the Target Spectrum.

##  Appendix B Application example

In this subsection we present an application example of the simplest case
scenario. A more complete example can be found in the Github repository of the
project111https://github.com/Gaiana/nirdust. 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 a FITS 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
$K_{long}$ filter. All this steps can be performed with the functionalities
described in the previous subsection. The impact of these operations on the
best-fit 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. As
well, the contribution of the $\gamma$ term can be limited through the
gamma_target_fraction parameter. In this example the bounds, the amount of
iterations, and the limit on the $\gamma$ term are initially given to the
function. In Fig. 8 an example of the output plot is presented.

⬇

>>> 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()

Figure 8: 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 dataset described in Sec. 4.2. The individual components involved
in the modeling process of equation 2, labeled as $\alpha$, $\beta$ and
$\gamma$ 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.

###  B.1 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 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 are enforced using the Flake-8 tool 222https://flake8.pycqa.org,
that 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-testing 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$\%$ 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 repository
333https://github.com/Gaiana/nirdust with a open-source MIT-License
444https://opensource.org/licenses/MIT The committed versions of code are
automatically tested with a continuous-integration workflow implemented with
GitHub Actions555https://github.com/Gaiana/nirdust/actions (Fowler & Foemmel,
2006), and the documentation is built from the repository and made public in
the read-the-docs service666https://nirdust.readthedocs.io/ (Holscher et al.,
2010).

The project is available in the Python-Package-Index
(PyPi)777https://pypi.org/project/nirdust/ throught the pip installer
888https://pypi.org/project/pip/. Finally, the project is registered in the
Astrophysics Source Code Library (ASCL, Allen & Schmidt,
2015)999https://ascl.net/code/v/3419.