11institutetext: Department of Physics and Astronomy ‘Galileo Galilei’,
University of Padova; Vicolo dell’Osservatorio 3, I-35122 Padova, Italy
11email<EMAIL_ADDRESS>22institutetext: INAF – Osservatorio
Astronomico di Padova; Vicolo dell’Osservatorio 5, I-35122 Padova, Italy
33institutetext: School of Physical Sciences, Faculty of Science, Technology,
Engineering and Mathematics, The Open University; Walton Hall, Milton Keynes,
MK7 6AA, UK

# Introducing MADYS: the Manifold Age Determination for Young Stars

I. Isochronal age estimates and model comparison

V. Squicciarini 1122 M. Bonavita 1133

(Received –; accepted –)

###### Abstract

Context. The unrivalled astrometric and photometric capabilities of the Gaia
mission have given new impetus to the study of young stars: both from an
environmental perspective, as members of comoving star-forming regions, and
from an individual perspective, as targets amenable to planet-hunting direct-
imaging observations.

Aims. In view of the large availability of theoretical evolutionary models,
both fields would benefit from a unified framework that allows a
straightforward comparison of physical parameters obtained by different
stellar and substellar models.

Methods. To this aim, we developed madys, a flexible Python tool for age and
mass determination of young stellar and substellar objects. In this first
release, madys automatically retrieves and cross-matches photometry from
several catalogs, estimates interstellar extinction, and derives age and mass
estimates for individual objects through isochronal fitting.

Results. Harmonizing the heterogeneity of publicly-available isochrone grids,
the tool allows to choose amongst 16 models, many of which with customizable
astrophysical parameters, for a total of $\sim 110$ isochrone grids. Several
dedicated plotting function are provided to allow an intuitive visual
perception of the numerical output.

Conclusions. After extensive testing, we have made the tool publicly available
(https://github.com/vsquicciarini/madys). We demonstrate here the capabilities
of madys, summarising already published results as well providing several new
examples.

###### Key Words.:

planets and satellites: fundamental parameters – stars: fundamental parameters
– methods: data analysis

##  1 Introduction

The advent of Gaia (Gaia Collaboration et al., 2016) has brought our view of
the Galaxy to its grandest level of sharpness, opening up tremendous
possibilities that stretch from the direct detection of exoplanets (e.g.,
Bonavita et al., 2022) to the realm of galactic archaeology (e.g., Helmi et
al., 2018). The availability of precise photometric and distance measurements
for $\sim 2\cdot 10^{9}$ stars enabled by Gaia EDR3 (Gaia Collaboration et
al., 2021) has paved the way to precise large-scale measurements of stellar
luminosity and effective temperature, which in turn allow to discern exquisite
features inside theoretical Hertzsprung-Russell diagrams (HRD) or their
observational counterpart, color-magnitude diagrams (CMD) (Gaia Collaboration
et al., 2018a).

The position of a star in a CMD, corrected for distance and interstellar
extinction, can in principle be compared with theoretical grids derived from
stellar evolution models to derive several astrophysical parameters, including
mass, age and metallicity. However, this widely-used technique, known as
”isochrone fitting”, is weakened by a strong degeneracy between age and
metallicity (Howes et al., 2019). Even when metallicity is independently
assessed, obtaining precise ages is an extremely challenging task for main-
sequence (MS) stars (Soderblom, 2010): the reason lies in the feeble variation
of luminosity and effective temperature over main-sequence lifetime, keeping
the isochrones very close to one another in the HRD.

Substantial efforts have thus been undertaken to explore complementary age-
dating techniques111Each one focusing on specific phases of stellar evolution
and/or mass ranges; see Barrado (2016) for a review. such as gyrochronology
(Barnes, 2007; Epstein & Pinsonneault, 2014), chromospheric activity (Mamajek
& Hillenbrand, 2008; Zhang et al., 2021), lithium depletion (Soderblom et al.,
2014) and, very promisingly, asteroseismology (e.g. Rodrigues et al., 2014;
Martig et al., 2016; Mackereth et al., 2021). Over the last few years, first
steps have been taken in combining infrared photometry from 2MASS (Skrutskie
et al., 2006) and AllWISE (Cutri et al., 2021) with spectroscopic constraints
on effective temperature and surface gravity to consistently derive stellar
parameters for millions of main-sequence and post-main sequence stars (e.g.,
Steinmetz et al., 2020; Mints, 2020; Fouesneau et al., 2022). The conjunction
of isochrone fitting and gyrochronology is the idea behind stardate (Angus et
al., 2019b), which shows that the integration of rotational periods into the
classical isochrone-based framework can lead to a gain of up to a factor 3 in
the accuracy on age for main-sequence FGK stars.

In addition to this, the exquisite precision of Gaia’s proper motion
measurements has been enormously beneficial to the study of stellar moving
groups, associations and clusters, leading to the compilation of large
catalogs of confirmed members for known regions (e.g., Cantat-Gaudin et al.,
2018; Gagné et al., 2018; Luhman & Esplin, 2020) and even to the discovery of
new ones (e.g. Cantat-Gaudin et al., 2019b). A complete census of star-forming
regions is, in turn, the first step toward resolving kinematic substructures
within them and connecting these structures with star formation history (e.g.
Kuhn et al., 2019; Cantat-Gaudin et al., 2019a; Pang et al., 2022). Again, the
need for differential age estimates within star-forming regions implies
precise individual age determinations; the picture is further complicated by
the fact that, at such young ages ($t\lesssim 100$ Myr), a large fraction of
the stars is still in the pre-main sequence (PMS) phase.

In principle, the dependence of luminosity on age during the PMS is steep
enough to allow a simultaneous determination of age and mass; indeed, the
spread between age predictions by different models – due to factors like
initial helium abundance, metallicity, mixing length, convective core
overshooting (Prada Moroni et al., 2016)– is generally acceptable ($\sim
20\%$) for F-G stars (Soderblom, 2010). However, the accuracy rapidly degrades
at later spectral types, so that the inter-model spread can be as high as a
factor of 4-5 for late-K and M stars (Soderblom, 2010) as a consequence of
rotation, activity and magnetic fields (Feiden, 2016b).

Not surprisingly, the problem continues to exacerbate below the hydrogen-
burning limit ($\sim 0.08M_{\odot}$) which separates stars from brown dwarfs
(Spiegel & Burrows, 2012). Thanks to the development of high-contrast imaging
facilities, young luminous brown dwarfs and giant exoplanets are being
increasingly found both in isolation (Miret-Roig et al., 2022) and as
companions to members of moving groups and associations (e.g., Vigan et al.,
2021). In this case, a simultaneous isochronal estimate of age and mass is no
more feasible, and independent constraints are needed to lift the degeneracy
between age and mass. Given the importance of the derived properties of these
young substellar objects, to study the low-mass end of the stellar initial
mass function (IMF) on the one hand (e.g. Kirkpatrick et al., 2019), and
exoplanet demographics on the other hand (e.g., Nielsen et al., 2019), it
becomes a crucial task to compare the different predictions done by different
models. In fortunate cases, mass estimates can be compared to model-
independent dynamical masses, possibly disentangling among formation
mechanisms (e.g. Marleau & Cumming, 2014; Brandt et al., 2021a).

With an eye on the study of star-forming regions and the other on directly-
imaged substellar objects, we developed the Manifold Age Determination of
Young Stars (madys): a tool that allows a robust inference of the properties
of stellar (substellar) objects based on stellar (substellar) evolution
models, and a straightforward comparison between different suites of models.
In this first paper we present the first, fully documented, public release of
madys222A preliminary version of madys has been already used, in in several
publications (Janson et al., 2021; Bonavita et al., 2021; Squicciarini et al.,
2021, 2022), now available at https://github.com/vsquicciarini/madys333The
current version on GitHub is an alpha version, only comprising the isochrones
used in this paper. and the underlying algorithm for isochronal age and mass
determination. We defer to a future paper the implementation of an indirect
method for age determination in star-forming regions, based on empirical
kinematic properties (see Squicciarini et al., 2021).

This paper is structured as follows. Section 2 introduces the key concepts
behind madys, with a special focus on data retrieval and extinction estimates;
the algorithm for age and mass determination is presented in Section 3; a few
applications are then provided in Section 4. In Section 5, we discuss the
strengths and the limitations of the tool, and anticipate its future
developments. Finally, in Section 6 we give a short summary of our results and
the main conclusions.

##  2 Introducing MADYS

The Manifold Age Determination of Young Stars (madys), written in Python, is a
tool that allows to derive the properties of arbitrarily large lists of young
objects by comparison with several published evolution models. madys combines
extensive cross-match capabilities, a careful assessment of photometric data,
the ability to estimate interstellar extinction in hundreds of photometric
bands, and the homogenization of a large collection of stellar and substellar
evolutionary models. The derivation of ages and masses can be performed under
several configurations depending on the science case, allowing for the
presence of multi-modal posterior distributions.

The tool can be used in two different modes, depending on the kind of input
provided:

  * •

a list of object names (mode 1);

  * •

a table containing object names, photometry and, optionally, astrometry (mode
2).

Several examples of code execution are provided in a dedicated Jupyter
Notebook within the GitHub repository
(https://github.com/vsquicciarini/madys/blob/main/examples/MADYS_examples.ipynb).

Generally speaking, the execution of the program triggers the following chain
of operations:

  * •

retrieval of photometric and astrometric data (only in mode 1, see Section
2.1);

  * •

selection of reliable photometry (Section 2.2);

  * •

estimation of interstellar extinction in all the bands of interest (Section
2.3);

  * •

creation of the final database to be used for age and mass determination
(Section 3.2).

The estimation of physical parameters is not done during initialization, but
rather by calling a dedicated method that acts upon the database: in this way,
it is possible to inspect data, to carefully decide the (sets of) theoretical
model(s) suitable to the science case (Section 3.1), and to repeat multiple
times the analysis of the same database.

###  2.1 Data retrieval

####  2.1.1 Mode 1

Building on the capabilities of
astroquery444https://astroquery.readthedocs.io/en/latest/ and
tap555https://github.com/mfouesneau/tap. to handle existing cross-matches
between Gaia and other catalogs (Marrese et al., 2019), madys queries the Gaia
Archive to return a single catalog containing astrometric, kinematic and
photometric information from Gaia DR2 (Gaia Collaboration et al., 2018b), Gaia
EDR3 (Gaia Collaboration et al., 2021) and 2MASS (Skrutskie et al., 2006).
Optionally, AllWISE (Cutri et al., 2021) and Pan-STARRS DR1 (Chambers et al.,
2016) can be added to the query as well666The integration of the forthcoming
Gaia DR3 into the described querying scheme is expected in the full version of
the program..

Although it is recommended to use Gaia DR2 or Gaia EDR3 IDs, it is also
possible to use other naming conventions. In the latter case, input names are
converted into their Gaia DR2 counterparts through a query of the SIMBAD
database (Wenger et al., 2000).

The results of the query are provided as a table whose $i$-th row corresponds
to the $i$-th star of the initial list of objects to allow an unambiguous
pairing of input and output data. In other words, the query algorithm of madys
is able to handle the cases in which more than one Gaia EDR3 source is
associated to the same Gaia DR2 source (or vice versa), selecting as best-mach
the source with Gaia DR2 ID = Gaia EDR3 ID or, if missing, the source in the
VizieR catalog (Ochsenbein et al., 2000) having the closest $G$, $G_{BP}$ and
$G_{RP}$ photometry in a suitable region accounting for the astrometric motion
of the source over the two epochs 777In the latter case, the source $X$ is
considered a cross-match to source $0$ only if $|G_{0}-G_{X}|<0.2$ mag and
$|G_{BP,0}-G_{BP,X}|<0.2$ mag and $|G_{RP,0}-G_{BP,X}|<0.2$ mag.. This
stratagem is able to find a cross-match for some high-proper motion stars
which are not paired by the SIMBAD database. Likewise, a small ($\sim 0.5\%$)
fraction of missing 2MASS matches is recovered by indirectly exploiting the
AllWISE-2MASS cross-match, or –if explicitly required– by directly inspecting
the SIMBAD database.

The typical speed of the query is about $\sim 100-150$ stars s-1, meaning that
a list of 1000 objects is fully recovered within few seconds. Large lists of
objects ($\sim 10^{4}-10^{5}$) are handled efficiently by splitting the query
into smaller chunks and later reassembling the results in the original order.

In any case, the resulting database always comprises data from Gaia DR2 and
Gaia EDR3. Parallaxes from Gaia EDR3888If the parallax for a source is present
in DR2 but not in EDR3, values from Gaia DR2 are used. and extinctions
(Section 2.3) are combined with apparent photometry to get absolute
magnitudes. Quality flags from selected surveys are retained with the aim of
identifying reliable photometric data (Section 2.2).

####  2.1.2 Mode 2

In mode 2, a table containing full information needed for age and mass
determination is provided. This mode is thought for objects that are not
present in Gaia, such as self-luminous direct-imaged exoplanets and brown
dwarfs.

Minimum requirements, in this case, consist of a column of object names and a
column with magnitudes. If parallaxes are provided, input magnitudes are
considered as apparent; otherwise, they are considered as absolute. By
providing two columns with equatorial or galactic coordinates, it is possible
for the program to evaluate interstellar extinction in the direction of the
object(s) and to take it into account (see Section 2.3). Alternatively, values
for the E(B-V) color excess can be manually provided.

More than 250 photometric filters are available in this mode, meaning that
there is at least one theoretical model which they can be compared to. The
full list of filters –including, e.g., the full suite of the James Webb Space
Telescope (JWST, Gardner et al., 2006) filters– can be found in the
documentation.

###  2.2 Selection of appropriate photometry values

This Section describes the conditions for a photometric measurement to be
retained in the final database.

By default, madys’ mode 1 collects photometric measurements from Gaia DR2/EDR3
($G$, $G_{BP}$, $G_{RP}$) and 2MASS ($J$, $H$, $K_{s}$). Gaia EDR3 $G$
magnitudes are corrected by adopting the prescriptions by Riello et al.
(2021). As regards $G_{BP}$ and $G_{RP}$, which are known to be intrinsically
much more sensitive than $G$ to contamination from nearby sources or from the
background (Evans et al., 2018), the phot_bp_rp_excess_factor is used as a
proxy to evaluate the quality of photometric measurements. Following Riello et
al. (2021), a color-independent corrected $BP/RP$ excess factor was defined
for both Gaia DR2 and Gaia EDR3:

$C^{*}=C+k_{0}+k_{1}\Delta G+k_{2}\Delta G^{2}+k_{3}\Delta G^{3}+k_{4}G$ (1)

where $\Delta G=(G_{BP}-G_{RP})$.

The $BP/RP$ excess factor has an expected value of 0 for well-behaved sources
at all magnitudes but, when considering subsamples of stars with similar
brightness, it tends to widen out for fainter $G$; a varying standard
deviation $\sigma(G)$ can be defined (Riello et al., 2021) as:

$\sigma_{C^{*}}(G)=c_{0}+c_{1}\cdot G^{m}$ (2)

We exclude $G_{BP}$ and $G_{RP}$ magnitudes with a corrected excess factor
larger, in absolute value, than 3 $\sigma_{C^{*}}(G)$. Values for the
constants for Eq. 1\- 2 are taken from Riello et al. (2021) for EDR3 and
Squicciarini et al. (2021) for DR2, and are provided in Table 1.

Table 1: Adopted values for Eq. 1-2. DR2 | $\Delta G<0.5$ | $\text{ }0.5\leq\Delta G<3.5$ | $\Delta G\geq 3.5$
---|---|---|---
$k_{0}$ | -1.121221 | -1.1244509 | -0.9288966
$k_{1}$ | +0.0505276 | +0.0288725 | -0.168552
$k_{2}$ | -0.120531 | -0.0682774 | 0
$k_{3}$ | 0 | 0.00795258 | 0
$k_{4}$ | -0.00555279 | -0.00555279 | -0.00555279
$c_{0}$ | 0.004
$c_{1}$ | $8\cdot 10^{-12}$
$m$ | 7.55
EDR3 | $\Delta G<0.5$ | $\text{ }0.5\leq\Delta G<4$ | $\Delta G\geq 4$
$k_{0}$ | 1.154360 | 1.162004 | 1.057572
$k_{1}$ | 0.033772 | 0.011464 | 0.140537
$k_{2}$ | 0.032277 | 0.049255 | 0
$k_{3}$ | 0 | -0.005879 | 0
$k_{4}$ | 0 | 0 | 0
$c_{0}$ | 0.0059898
$c_{1}$ | $8.817481\cdot 10^{-12}$
$m$ | 7.618399

From 2MASS and AllWISE, only sources with photometric flag ph_qual=’A’ are
kept. If needed, a different value for the worst quality flag still considered
reliable can be selected.

###  2.3 Interstellar extinction

The estimate of extinction (reddening) in a given band (color) is performed by
integrating along the line of sight a suitable 3D extinction map. The
integration algorithm draws a line from the position of the Sun toward that of
the star of interest; the value of each pixel crossed by the line is weighted
according to the portion of the total distance spent by the line in the pixel
itself. This method ensures a rapid and precise evaluation of the integral,
allowing 10000 stars to be handled in $\sim 1$ s under typical PC
performances.

Two extinction maps can be selected:

  * •

the STILISM 3D extinction map by Lallement et al. (2019): a Sun-centered
(6000x6000x800) pc grid, with step equal to 5 pc;

  * •

the Galactic extinction catalog by Leike et al. (2020): a Sun-centered
(740x740x540) pc grid with step equal to 1 pc.

Since the file with the selected map must be present in the local path where
madys has been installed, the errors on the derived estimates –which would
require the download of additional files– are currently not returned by the
program.

Coordinate transformations from the equatorial/galactic frame to the right-
handed galactocentric frame (i.e., a Cartesian galactic frame) is performed by
means of the astropy package999Default parameters from the ’pre-v4.0’ are
used: galcen_distance=8.3 kpc, galcen_v_sun=(11.1, 232.24, 7.25) km s-1,
z_sun=27.0 pc, roll=0.0 deg. See
https://docs.astropy.org/en/stable/api/astropy.coordinates.Galactocentric.html
for details. (Astropy Collaboration et al., 2013).

As a test for the accuracy of the algorithm, we provide in Figure 1 a visual
comparison between the integrated absorption in the Upper Scorpius region
(already used in Squicciarini et al., 2021) obtained through the map by Leike
et al. (2020) and the intensity Stokes map returned by the Plank satellite
(Planck Collaboration et al., 2020) in the far infrared ($\nu=857$ GHz
$=350\leavevmode\nobreak\ \mu$m). Given the large galactic latitude of the
region ($l\in[10^{\circ},30^{\circ}]$), we expect the large majority of the
integrated intensity in the Planck image to be produced by the association,
with only a negligible background contribution. Indeed, the agreement between
the two images is excellent.

Figure 1: Integrated extinction toward Upper Scorpius. (A) G-band extinction
map produced by madys by integrating the 3D map by Leike et al. (2020) up to a
distance $d=160$ pc. (B) Intensity Stokes map at $350\leavevmode\nobreak\
\mu$m produced by the Planck satellite.

The conversion between extinction and reddening is mediated by a total-to-
selective absorption ratio $R=3.16$ (Wang & Chen, 2019). The extinction law is
obtained by combining the extinction law by Wang & Chen (2019) in the range
$[0.3,2]$ $\mu$m and the diffuse average extinction by Gordon et al. (2021) in
the range $[6.5,40]$ $\mu$m; a linear combination of the two is used in the
intermediate range $[2,6.5]$ $\mu$m (Figure 2):

$\frac{A_{\lambda}}{A_{V}}=\cases{\displaystyle}\sum_{i=0}^{7}b_{i}\leavevmode\nobreak\
\xi^{i}&\leavevmode\nobreak\ \lambda\in[0.3,1]\leavevmode\nobreak\ \mu m\\\
h_{2}\lambda^{\beta_{2}}\equiv f_{2}(\lambda)\leavevmode\nobreak\
\lambda\in[1,2]\leavevmode\nobreak\ \mu m\\\
f_{2}(\lambda)+q(\lambda)f_{4}(\lambda)\leavevmode\nobreak\
\lambda\in[2,6.5]\leavevmode\nobreak\ \mu m\\\
h_{4}\lambda^{\beta_{4}}+S_{1}D_{1}(\lambda)+S_{2}D_{2}(\lambda)\equiv
f_{4}(\lambda)\leavevmode\nobreak\ \lambda\in[6.5,40]\leavevmode\nobreak\ \mu
m\\\ {}$ (3)

where:

$\xi=\frac{1}{\lambda}-1.82\leavevmode\nobreak\ \mu\text{m}^{-1}$ (4)
$q(\lambda)=\frac{\lambda-2\leavevmode\nobreak\
\mu\text{m}}{6.5\leavevmode\nobreak\ \mu\text{m}-2\leavevmode\nobreak\
\mu\text{m}}$ (5)

$D_{1}$ and $D_{2}$ are two modified Drude profiles, used to model the
silicate absorption features at $\sim 10\leavevmode\nobreak\ \mu\text{m}$ and
$\sim 20\leavevmode\nobreak\ \mu\text{m}$:

$D(\lambda)=\frac{(\gamma(\lambda)/\lambda_{0})^{2}}{((\lambda/\lambda_{0}-\lambda_{0}/\lambda)^{2}+(\gamma(\lambda)/\lambda_{0})^{2})}$
(6)

while $\gamma(\lambda)$ is in turn given by:

$\gamma(\lambda)=\frac{2\gamma_{0}}{1+\exp(a_{0}(\lambda-\lambda_{0}))}$ (7)

(Gordon et al., 2021). We list in Table 2 all the coefficients from Eq. 3-7,
where ($\gamma_{1}$, $\lambda_{1}$, $a_{1}$) and ($\gamma_{2}$, $\lambda_{2}$,
$a_{2}$) indicate the coefficients for $D_{1}(\lambda)$ and $D_{2}(\lambda)$,
respectively.

Table 2: Adopted values for the coefficients in Eq. 3-7. Name | Value | Name | Value
---|---|---|---
$b_{0}$ | 1 | $h_{4}$ | 0.366
$b_{1}$ | 0.7499 | $\beta_{4}$ | -1.48
$b_{2}$ | -0.1086 | $S_{1}$ | 0.06893
$b_{3}$ | -0.08909 | $S_{2}$ | 0.02684
$b_{4}$ | 0.02905 | $\lambda_{1}$ | 9.865
$b_{5}$ | 0.01069 | $\gamma_{1}$ | 2.507
$b_{6}$ | 0.001707 | $a_{1}$ | -0.232
$b_{7}$ | -0.001002 | $\lambda_{2}$ | 19.973
$h_{2}$ | 0.3722 | $\gamma_{2}$ | 16.989
$\beta_{2}$ | -2.07 | $a_{2}$ | -0.273

The adopted extinction law goes farther in the mid-infrared than widely-used
parametrizations, as those offered by the extinction
package101010ttps://extinction.readthedocs.io/en/latest/, delving into
wavelength ranges amenable to forthcoming JWST observations. Individual
extinction coefficients $A_{\lambda}$ are directly taken from Table 3 of Wang
& Chen (2019) whenever possible, or computed through Eq. 3 adopting as
$\lambda$ the mean wavelength indicated by the SVO Filter Profile Service
(Rodrigo et al., 2012; Rodrigo & Solano, 2020).

Figure 2: Adopted extinction law (solid line). Blue dots represent absorption
coefficients in the visible and NIR (Wang & Chen, 2019, Table 3), green dots
in the MIR (Gordon et al., 2021, Table 8).

##  3 Age and mass determination

###  3.1 Loading isochrone tables

As already mentioned, the determination of ages and masses in this first
release of madys is performed via isochronal fitting, comparing the database
obtained in Section 2.2 to the selected set of isochrones. We will
interchangeably refer to these estimates as ”photometric” or ”isochronal”
estimates.

Overcoming the multifarious conventions of isochrone tables found in the
literature, madys employs an extended version of the evolution
routine111111https://github.com/avigan/Python-utils/tree/master/vigan/astro
that currently supports 16 different stellar and/or substellar evolutionary
models. Some of these models offer the additional possibility to customize
astrophysical parameters such as metallicity, rotational velocity, helium
content, alpha enhancement and the spot fraction, for a total of $\sim 110$
different isochrone tracks (Table 3.

Table 3: Isochrone models currently supported by madys with their bibliographic references. Adopted values for solar metallicity ($Z_{\odot}$) and the initial helium abundance ($Y_{0}$) are reported (u=unknown), together with mass ($M$) and age ranges ($t$). Customizable parameters (c.p.): metallicity (m), helium content (Y), alpha enhancement ($\alpha$), rotational velocity (v), spot fraction (s). Name | $Y_{0}$ | $Z_{\odot}$ | $M$ | $t$ | c. p. | reference
---|---|---|---|---|---|---
|  |  | $M_{\odot}$ | Myr |  |
ATMO2020 | 0.275 | 0.0169 | [0.001, 0.075] | $[10^{0},10^{4}]$ | — | Phillips et al. (2020)
B97 | u | u | [$10^{-3}$,0.04] | $[10^{0},10^{4}]$ | — | Burrows et al. (1997)
BHAC15 | 0.271 | 0.0153 | [0.01,1.4] | [0.5, $10^{4}$] | — | Baraffe et al. (2015)
Geneva | 0.266 | 0.014 | [0.8,120] | [0.1, $10^{2}$] | — | Haemmerlé et al. (2019)
MIST | 0.2703 | 0.0142 | [0.1,150] | $[10^{-1},2\cdot 10^{4}]$ | m, $\alpha$, v | Dotter (2016); Choi et al. (2016)
PARSEC | 0.2485 | 0.01524 | [0.09,350] | $[10^{-1},10^{4}]$ | m | Marigo et al. (2017)
(PHOENIX) AMES-COND | u | u | [0.005,1.4] | $[1,1.2\cdot 10^{4}]$ | — | Allard et al. (2001)
(PHOENIX) AMES-Dusty | u | u | [0.005,1.4] | $[1,1.2\cdot 10^{4}]$ | — | Allard et al. (2001)
(PHOENIX) Bt-Settl | 0.2485 | 0.0153 | [0.01,1.4] | $[10^{0},10^{4}]$ | — | Allard (2016)
(PHOENIX) NextGen | 0.28 | 0.02 | [0.01,1.4] | $[10^{0},1.2\cdot 10^{4}]$ | — | Hauschildt et al. (1999)
PM13 | — | — | [0.075,27] | — | — | Pecaut & Mamajek (2013)
SB12 | u | u | [$10^{-3}$,$10^{-2}$] | [1,100] | — | Spiegel & Burrows (2012)
Sonora Bobcat | 0.2735 | 0.0153 | $[5\cdot 10^{-4},10^{-1}]$ | $[10^{-1},10^{4}]$ | — | Marley et al. (2021)
SPOTS | 0.2676 | 0.0165 | [0.1,1.3] | $[10^{0},4\cdot 10^{3}]$ | s | Somers et al. (2020)
STAREVOL | 0.269 | 0.0134 | [0.2,1.5] | $[1,1.26\cdot 10^{4}]$ | m, v | Amard et al. (2019)
YAPSI | 0.2775 | 0.0142 | [0.15,5] | $[0.5,10^{4}]$ | m, Y | Spada et al. (2017)

Mass and age ranges for the grid can be externally imposed; otherwise, the
program computes suitable mass limits that take into account both the
dynamical range of the model and the range of values expected from a rough
preliminary evaluation of the sample’s absolute Gaia EDR3 $G$-band and, if
applicable, 2MASS $K_{s}$-band magnitudes based on the tables by Pecaut &
Mamajek (2013).

After selecting a model, the program creates two vectors of masses and ages,
whose points are equally spaced on a logarithmic scale, and builds the
theoretical photometric grid. Linear interpolation between consecutive points
into the final grid is performed for every quantity, and no extrapolation is
attempted outside the grids. The usage of a fixed grid, rather than real-time
interpolation, was chosen for a twofold reason: to handle arbitrarily large
group of stars while keeping computation times reasonable, and to allow a
homogeneous treatment of statistical uncertainties. The spacing between
consecutive steps can be adjusted, but is by default significantly smaller
than any realistic uncertainty value.

Generally speaking, the choice of the isochrone set for a given sample should
be carefully pondered depending on the expected mass and age ranges of the
sample and on the photometric filters of interest. We notice that a few
isochrone sets do not come equipped with Gaia and 2MASS filters: hence, they
can only be used in mode 2. The program is not halted if it encounters a
filter that is not available for the selected model, but –after throwing a
warning– it neglects all the photometry provided in that filter. For this
reason, it is always recommended to look at the execution log produced by the
program.

###  3.2 Age and mass determination

For each object in the final database (Section 2.2), madys seeks the minimum
of a suitable $\chi^{2}$ function:

$\chi^{2}=\sum_{k}\left(\frac{M^{th}_{k}-M^{obs}_{k}}{\sigma_{M^{obs}_{k}}}\right)^{2}\equiv\sum_{k}s^{2}_{k}$
(8)

which can be thought as a 2D distance matrix with same shape as the age-mass
grid and elements:

$\chi^{2}_{ij}=\sum_{k}\left(\frac{M^{th}_{ijk}-M^{obs}_{k}}{\sigma_{M^{obs}_{k}}}\right)^{2}$
(9)

where $M^{th}_{ijk}$ is the theoretical magnitude in the $k$-th filter
corresponding to the $i$-th mass and $j$-th age of the model grid, $M_{obs,k}$
is the observed magnitude in the same filter and $\sigma_{M_{obs,k}}$ its
associated uncertainty. The sum is done only over the filters $k$ passing the
prescriptions:

  1. 1.

an error on the absolute magnitude smaller than 0.2 mag;

  2. 2.

a best-match $M^{th}_{i_{0}j_{0}k}$ such that
$|M^{th}_{i_{0}j_{0}k}-M^{obs}|<0.2$ mag.

Individual age ranges can be provided for each target, and this is
particularly useful when external constraints are available; the only caveat
is that the kind of input should be the same for every target. In particular,
if the age of each object star is explicitly imposed, or a triplet [optimal
age, minimum age, maximum age] is provided (case 1), a single filter is
sufficient for parameter estimation; conversely, if no age constraint is
given, or just a doublet [minimum age, maximum age] is provided (case 2), the
estimation is performed only if the following conditions are met:

  1. 3.

at least three filters passed the prescriptions 1. and 2.;

  2. 4.

after identifying the minimum $\chi^{2}$, its third smallest associated
$s^{2}_{k}$ ¡ 9, or alternatively its third smallest
$|M^{th}_{k}-M^{obs}_{k}|<0.1$ mag.

In order to embed photometric uncertainty into the final parameter estimate,
the procedure is repeated $q=1000$ times while randomly varying, using a Monte
Carlo approach, the apparent photometry and the parallax according to their
uncertainties (which are assumed to be distributed in a normal fashion).

In case 1, the age is not fitted, and the resulting mass distribution is
assumed to be unimodal: i.e., the masses corresponding to the
$16^{th}$,$50^{th}$ and $84^{th}$ percentile of the sample composed by the $q$
best-fit solutions are returned.

In case 2, the algorithm considers the possibility of a multimodal posterior
distribution for both age and mass. At each iteration $q$, having a minimum
$\chi^{2}=\chi^{2}_{q}$, the set of $(i,j)$ mass and age steps such that:

$\chi^{2}_{i,j}<\chi^{2}_{q}+\Delta\chi^{2}$ (10)

are collected and added to a single array, $\bar{P}$. We decided to adopt
$\Delta\chi^{2}=3.3$ as it defines the 68.3% confidence region around the
best-fit joint parameter estimate for a two-parameter model (see, e.g., Verde,
2010).

The final outcome of this procedure is an array of solutions, $\bar{P}$. The
”hottest points” are the indices recurring more frequently; each occurrence of
a point has an associated $\chi^{2}$, and this should be properly reflected
into the final weighted average. In general, the ensemble of points in
$\bar{P}$ will not be connected, meaning that multiple families of solutions
in the age-mass space can be possible.

An intuitive approach to identify these families consists in identifying
connected regions in the age-mass grid. In order to reduce the strong
dependence of the connection on random realizations of data perturbations, we
decided to define as ”attraction points” the points which appear at least in
the 10% of the interactions in $\bar{P}$. Each isolated attraction point
defines a family of solutions; a group of contiguous attraction points is
treated as a single attraction point located in the group’s center of mass,
hence defining a single family of solutions as well. The remaining points are
then assigned to the family of the closest attraction point.

Each family of solutions $p$ corresponds, from a physical perspective, to a
different physical solution; its associated age and mass estimates
($t_{p}$,$m_{p}$) are defined as the average of the $i$-th mass and the $j$-th
age, weighted by a coefficient $1/\chi^{2}_{ij,q}$:

$\log_{10}{m_{p}}=\frac{\sum_{(i,j)\in
p}\log_{10}{m_{i}}\cdot(\chi^{2}_{ij,q})^{-1}}{\sum_{(i,j)\in
p}(\chi^{2}_{ij,q})^{-1}}$ (11)

$\log_{10}{a_{p}}=\frac{\sum_{(i,j)\in
p}\log_{10}{a_{j}}\cdot(\chi^{2}_{ij,q})^{-1}}{\sum_{(i,j)\in
p}(\chi^{2}_{ij,q})^{-1}}$ (12)

where, of course, points $(i,j)$ repeating in different iterations are summed
each time with a weight corresponding to the $\chi^{2}_{ij,q}$ of the $q$-th
iteration.

The variances associated to $\log_{10}{m_{p}}$ and $\log_{10}{a_{p}}$ are
given by:

$\sigma^{2}_{m_{p}}=\frac{\sum_{(i,j)\in
p}(\log_{10}{m_{i}}-\log_{10}{m_{p}})^{2}\cdot(\chi^{2}_{ij,q})^{-1}}{\sum_{(i,j)\in
p}(\chi^{2}_{ij,q})^{-1}}$ (13) $\sigma^{2}_{a_{p}}=\frac{\sum_{(i,j)\in
p}(\log_{10}{a_{i}}-\log_{10}{a_{p}})^{2}\cdot(\chi^{2}_{ij,q})^{-1}}{\sum_{(i,j)\in
p}(\chi^{2}_{ij,q})^{-1}}$ (14)

Couples of solutions ($p_{1}$,$p_{2}$) that are consistent with representing
the same solution, i.e. with:

$\Delta
d=\frac{(\log_{10}{m_{p_{1}}}-\log_{10}{m_{p_{2}}})^{2}}{\sigma^{2}_{m_{p_{1}}}+\sigma^{2}_{m_{p_{2}}}}+\frac{(\log_{10}{a_{p_{1}}}-\log_{10}{a_{p_{2}}})^{2}}{\sigma^{2}_{a_{p_{1}}}+\sigma^{2}_{a_{p_{2}}}}<8$
(15)

are merged together. The outcome of the process is a set of solutions
$\\{p\\}$, each one bearing a fraction of the total region of the solutions
$\bar{P}$ equal to:

$w_{p}=\frac{\sum_{(i,j)\in
p}(\chi^{2}_{ij,q})^{-1}}{\sum_{(i,j)\in\bar{P}}(\chi^{2}_{ij,q})^{-1}}$ (16)

The solution with the maximum $w_{p}$ is returned as the best-fit solution,
but the other solutions can be inspected as well. Both the $\chi^{2}$ map for
nominal photometric values and the weight map $W$, defined as the 2D matrix
with elements:

$w_{ij}=\frac{1}{\sum_{(i,j)\in\bar{P}}(\chi^{2}_{ij,q})^{-1}}$ (17)

referring instead to the whole fitting process, can be returned and plotted
through a dedicated method.

##  4 Applications

###  4.1 Age substructures

The ability of madys to handle thousands of stars at once makes it
particularly suited to the systematic study of star-forming regions.

As a possible application of the code, we compute here the age of confirmed
members of the Scorpius-Centaurus association. The association, that is the
nearest star-forming region to the Sun, is classically divided into three
subgroups: Upper Scorpius (US), Upper Centaurus-Lupus (UCL) and Lower
Centaurus-Crux (LCC) (de Zeeuw et al., 1999).

We start from the list of bona fide Scorpius-Centaurus members compiled by
Damiani et al. (2019) using Gaia DR2 data. In order to define the sub-regions,
we employ classical coordinate boundaries as in de Zeeuw et al. (1999) and
subsequent works: for US, $l\in[343^{\circ},360^{\circ}]$,
$b\in[10^{\circ},30^{\circ}]$; for UCL, $l\in[313^{\circ},343^{\circ}]$,
$b\in[0^{\circ},25^{\circ}]$; for LCC, $l\in[280^{\circ},313^{\circ}]$,
$b\in[-10^{\circ},23^{\circ}]$.

Starting from Gaia DR2 IDs, madys recovers the photometry and computes
extinction values as described in Section 2.3. The age and mass determination,
initialized with only a modest constraint on age ($t\in[1,300]$ Myr), is done
here with the BHAC15 models (Baraffe et al., 2015).

A visual inspection of the ($G_{BP}-G_{RP},G$) CMD shows that some stars
appear to be too old to be members of the association, and indeed they have
fitted ages $\gtrsim 100$ Myr. Therefore, we exclude the stars with best-fit
ages greater than 60 Myr; for the few stars with multiple possible solutions,
i.e. with an overlap between a pre-MS solution and an evolved MS solution, we
pick the youngest one.

Figure 3 shows the derived ages and masses for the three subgroups. Group
ages, computed as the $16^{th}$,$50^{th}$ and $84^{th}$ percentile of the age
distribution of their members, are:

$US:6.6_{-3.0}^{+5.6}\text{ Myr,}$ $UCL:9.1_{-3.7}^{+4.8}\text{ Myr,}$
$LCC:8.7_{-3.2}^{+5.6}\text{ Myr.}$

Figure 3: Derived age and mass distribution of the Sco-Cen sample, restricted
to stars with $M>0.8M_{\odot}$. The insets represent a kernel-density estimate
of the underlying distributions. The clump of sources at the upper mass end is
an artifact due to the sharp mass cut of BHAC15 isochrones at $1.4M_{\odot}$.

We recover some facts which are already known from the literature (see, e.g.
Pecaut & Mamajek, 2016):

  * •

US is younger than UCL and LCC;

  * •

a positive correlation exists between age and mass; in other words, M stars
appear younger than their F and G counterparts. Although an age spread between
the two cannot be completely ruled out, most of the observed spread is likely
due to a slowed-down contraction of low-mass stars caused by magnetic fields
(Feiden, 2016a). See Squicciarini et al. (2021) for a review;

  * •

indeed, if we restrict to stars with a best-fit $M>0.8M_{\odot}$, the results
become: $US:13.1_{-7.1}^{+8.3}\text{ Myr}$, $UCL:16.0_{-5.0}^{+6.9}\text{
Myr}$, $LCC:15.3_{-6.4}^{+7.9}\text{ Myr}$. The results are similar, both in
the median value and in the associated scatter, to the estimates by Pecaut &
Mamajek (2016).

With these caveats in mind, the possibility of computing individual age
estimates for pre-MS stars with madys opens up important opportunities for the
study of young star-forming regions, whose exquisite substructures are being
more and more connected with their star formation history (e.g., Kerr et al.,
2021; Krolikowski et al., 2021; Squicciarini et al., 2021).

###  4.2 Stellar physical parameters

Although by construction madys is able to return age estimates for main-
sequence stars (MS), we stress that in this case they should be regarded as
purely indicative. Nevertheless, the argument can be reversed: if external age
constraints are available, madys can return precise determination of stellar
parameters such as mass, effective temperature, radius and surface gravity for
large samples of stars.

As an example of this possibility, we recovered from the literature a
collection of stars with interferometrically measured angular diameters. Our
sample combines the six main-sequence stars studied by Huber et al. (2012) and
the full samples by Boyajian et al. (2012a, b), spanning a spectral range that
stretches from A- to M-type. We combined angular measurements with the latest
parallaxes from Gaia EDR3 to have radius estimates with median precision $\sim
1\%$.

Our parameter estimates are based on PARSEC isochrones; we applied only a
modest age constraint ($t\in[500,7000]$ Myr); with respect to metallicity, we
refer to [Fe/H] estimates from the original studies. Under the assumption
[Fe/H]$\approx$[M/H], we used for each star the isochrone grid with the
closest metallicity rather than interpolating 121212The available
metallicities for this example were:
[Fe/H]=[-1.0,-0.75,-0.5,-0.25,0.0,0.13,0.25,0.5,0.75,1.00]..

The results are shown in Figure 4. The mean and standard deviation of the
fractional difference between interferometric radii and those returned by
madys are +1% and 6%, respectively.

Figure 4: Comparison between photometric radius estimates obtained by madys
and interferometric radii from the literature: Huber et al. (2012) (H12),
Boyajian et al. (2012a) (B12a), Boyajian et al. (2012b) (B12b).

###  4.3 Mass of directly-imaged substellar companions

madys has been conceived to connect models spanning from young self-luminous
gas giants to old stars. Thanks to the complete coverage of JWST filters
offered by a handful of models, madys will be a valuable tool to study the
first data coming from the instrument, and to help characterize the properties
of newly-discovered objects.

Indeed, mode 2 is intended for objects that are not found in Gaia, such as
objects discovered in direct imaging studies, either in isolation or as
companions to stellar objects. In the latter case, madys can be used in two
steps of the chain: to obtain the age of the stellar target –either directly
or by exploiting its kinematic neighbors– and to use such information to
derive a mass estimate for the companion. A combination of indirect kinematic
constraints and literature data was recently applied in Squicciarini et al.
(2022) to derive age and mass estimates for both the stellar target and its
companions.

We present here an application of madys to the HR 8799 system (Marois et al.,
2008), one of the cornerstones of direct imaging studies. With four giant
planets detected around a $1.47^{+0.11}_{-0.08}M_{\odot}$ primary (Wang et
al., 2018), the system is a unique laboratory to test the accuracy of
substellar evolutionary models.

Several age estimates have been derived for the system in the recent
literature: $10-23$ Myr from Sepulveda & Bowler (2022), $33^{+7}_{-13}$ Myr
from Baines et al. (2012) or $42^{+6}_{-4}$ Myr from Bell et al. (2015); we
notice that the last estimate, based on the lithium-depletion boundary for the
Colomba association which the star appears to be a member of, is independently
indicated by madys when inspecting the kinematic neighborhood of the star.
Indeed, we identified three stars131313Gaia EDR3 2838213864935858816, Gaia
EDR3 2835796794780262912 and Gaia EDR3 2830197806693622272. with projected
separation ¡ 3 pc and tangential velocity difference ¡ 3 km s-1: all of them
have a best-fit mass $\in[0.7,1]M_{\odot}$ and age $\sim 40$ Myr.

Nevertheless, we conservatively start from two possible age intervals,
$t\in[10,23]$ Myr and $t\in[30,60]$ Myr, to compare our estimates with already
published results. As a consequence of the uncertainty on age, we expect the
model-dependent uncertainty on the derived photometric mass estimates to be
broadened.

Table 4 reports literature estimates for the masses of the four planets,
obtained with photometric or dynamical methods, together with new estimates
obtained by madys. We collected contrasts measurements from Zurlo et al.
(2016) in the SPHERE bands $J$ ($\lambda_{peak}=1.245\leavevmode\nobreak\
\mu$m), $H_{2}$ ($\lambda_{peak}=1.593\leavevmode\nobreak\ \mu$m), $H_{3}$
($\lambda_{peak}=1.667\leavevmode\nobreak\ \mu$m), $K_{1}$
($\lambda_{peak}=2.110\leavevmode\nobreak\ \mu$m) and $K_{2}$
($\lambda_{peak}=2.251\leavevmode\nobreak\ \mu$m), and combined them to 2MASS
magnitudes and Gaia EDR3 parallax for HR 8799 to obtain absolute magnitudes.

Mass estimates were obtained through four models: namely, AMES-Dusty, AMES-
Cond, ATMO2020 (in particular, the chemical equilibrium grid) and Sonora
Bobcat 141414For ATMO2020 and Sonora Bobcat, which currently lack SPHERE
filters, we employed theoretical magnitudes in the closest photometric system
available: the Mauna Kea Observatories photometric system (MKO, Tokunaga et
al., 2002), and 2MASS, respectively. In particular, $J_{\text{SPHERE}}\sim
J_{\text{MKO}}\sim J_{\text{2MASS}}$, $K_{\text{SPHERE}}\sim
K_{\text{MKO}}\sim K_{s,\text{2MASS}}$,
$0.5\cdot(H_{2,\text{SPHERE}}+H_{3,\text{SPHERE}})\sim H_{\text{MKO}}\sim
H_{\text{2MASS}}$..

The results are also summarised in Figure 5. While the results of photometric
estimates can significantly differ from one another even in the same age
window, tighter dynamical constraints151515The small errorbar of the dynamical
mass estimates by Goździewski & Migaszewski (2020) is a consequence of the
assumption that the planets are in an exact 8:4:2:1 mean-motion resonance, and
should therefore be taken with caution. coming from thorough astrometric
follow-up in the next few years will help distinguishing among them, shedding
light into the still poorly-constrained cooling timescale of young self-
luminous Super Jupiters.

Table 4: Mass estimates for the planets in the HR 8799 system. We compare the best-fit results from madys with photometric and dynamical masses taken from the literature. Two age ranges ($t\in 10-23$ Myr, $t\in 30-60$ Myr) and four models are used for the independent estimates. source | age | b | c | d | e | method
---|---|---|---|---|---|---
| Myr | $M_{\mathrm{J}}$ | $M_{\mathrm{J}}$ | $M_{\mathrm{J}}$ | $M_{\mathrm{J}}$ |
Marois et al. (2008) | $30-60$ | $5-7$ | $7-10$ | $7-10$ | $7-10$ | photometric
Wilner et al. (2018) | — | $5.8^{+7.9}_{-3.1}$ | — | — | — | dynamical
Wang et al. (2018) | $42\pm 5$ | $5.8\pm 0.5$ | $7.2_{-0.7}^{+0.6}$ | $7.2_{-0.7}^{+0.6}$ | $7.2_{-0.7}^{+0.6}$ | photometric
Goździewski & Migaszewski (2020) |  | $5.7\pm 0.4$ | $7.8\pm 0.5$ | $9.1\pm 0.2$ | $7.4\pm 0.6$ | dynamical
Brandt et al. (2021b) | $42^{+24}_{-16}$ | — | — | — | $9.6^{+1.8}_{-1.9}$ | dynamical
Sepulveda & Bowler (2022) | $10-23$ | $2.7-4.9$ | $4.1-7.0$ | $4.1-7.0$ | $4.1-7.0$ | photometric
Zurlo et al. (2022) | — | $5.54-6.20$ | $6.84-8.10$ | $9.07-10.05$ | $7.11-10.66$ | dynamical
madys (this work) | $10-23$ | $6.9_{-1.2}^{+0.9}$ | $7.8_{-1.3}^{+1.0}$ | $7.5_{-1.3}^{+1.0}$ | $7.9_{-1.3}^{+1.0}$ | AMES-Dusty
$30-60$ | $10.7_{-1.1}^{+0.4}$ | $11.3_{-0.6}^{+0.3}$ | $11.1_{-0.8}^{+0.4}$ | $11.4_{-0.6}^{+0.3}$ | AMES-Dusty
$10-23$ | $3.8_{-0.7}^{+0.5}$ | $6.9_{-1.2}^{+0.9}$ | $7.5_{-1.4}^{+1.0}$ | $7.2_{-1.3}^{+1.0}$ | ATMO2020
$30-60$ | $6.5_{-1.0}^{+0.9}$ | $10.9_{-1.2}^{+0.8}$ | $11.4_{-1.0}^{+0.5}$ | $11.2_{-1.1}^{+0.6}$ | ATMO2020
$10-23$ | $4.0_{-0.8}^{+0.6}$ | $5.8_{-1.1}^{+0.9}$ | $6.0_{-1.1}^{+1.0}$ | $6.1_{-1.1}^{+1.0}$ | Ames-Cond
$30-60$ | $7.0_{-1.0}^{+0.9}$ | $9.9_{-1.3}^{+0.8}$ | $10.1_{-1.4}^{+0.6}$ | $10.4_{-1.4}^{+0.5}$ | Ames-Cond
$10-23$ | $4.5_{-0.8}^{+0.7}$ | $7.2_{-1.1}^{+1.0}$ | $7.7_{-1.1}^{+1.0}$ | $7.5_{-1.1}^{+1.0}$ | Sonora Bobcat
$30-60$ | $7.5_{-1.1}^{+0.8}$ | $11.0_{-1.0}^{+0.6}$ | $11.3_{-0.8}^{+0.4}$ | $11.2_{-0.9}^{+0.5}$ | Sonora Bobcat
Figure 5: Literature (black) and new (blue) mass estimates for the HR 8799
planets. Each panel refers to the eponymous planet. For the sake of clarity,
dynamical estimates (D) have been placed in the gray region, which visually
separates the estimates based on a younger (10-23 Myr) age (Y) from those
estimates based on an older (30-60 Myr) age (O). Different symbols for madys
estimates refer to different models.

##  5 Discussion

Since the algorithm behind age determination in madys is based on isochronal
fitting, the tool automatically inherits the same drawbacks of this technique,
which have been the subject of extensive discussion in the literature (see,
e.g., Soderblom, 2010; Barrado, 2016).

In particular, madys alone can not be used to understand if a star is young
($t\lesssim 100$ Myr) or not: a degeneracy exists between young pre-MS stars
and evolved stars that have left the MS, and this is naturally reflected into
different families of solutions that arise if the age is left completely
unconstrained. A young solution is to be preferred if independent youth
indicators (e.g., activity indicators such as X-ray, UV, $H_{\alpha}$
emission) are available.

A conceptually different youth indicator is the membership to a young star-
forming region. Indeed, the integration of kinematic information into madys
will be the subject of a second version of the tool. For the moment being,
madys can exploit existing lists of confirmed members of these regions to
unveil correlations between the star formation history and kinematic
substructures (Squicciarini et al., 2021).

A strong asset of madys is the ability to collect and handle photometric data
for thousands of stars. The age determination is rapid, taking about one
second per star under typical conditions. In this sense, our tool constitutes
a step forward in the automation of the process with respect to existing tools
such as PARAM (Rodrigues et al., 2017), ARIADNE (Vines & Jenkins, 2022),
stardate (Angus et al., 2019a) or isochrones (Morton, 2015); on the other
hand, unlike them, it does not currently allow to exploit additional
information coming from, e.g., asteroseismology or spectroscopy, during the
fitting process.

A second strength of madys is the possibility to employ a large collection of
stellar and substellar models, allowing the possibility to evaluate the impact
of different input physics into the final results. This is particularly
important not only for the pre-MS phase of stars, but also in the context of
direct imaging studies of cooling substellar objects, where there is still no
established standard on how photometric mass determinations are to be
obtained.

##  6 Conclusions

We introduced a Python tool, madys, aimed at obtaining photometric age and
mass estimates for arbitrarily large groups of young stellar or substellar
objects. The main strengths of the tool are:

  * •

the ability to query and cross-match different catalogs to yield highly
reliable catalogs of potentially large lists of objects;

  * •

the possibility to take interstellar extinction into account;

  * •

the ability to derive photometric ages and mass estimates by comparing dozens
of filters with a large suite of substellar or stellar evolutionary models;

  * •

the possibility to unambiguously compare the predictions of different models,
and to see the effect of changing astrophysical parameters;

  * •

the large plotting options for efficient data visualization.

These features give madys a large number of possible scientific applications,
such as:

  * •

the study of young star-forming regions, connecting kinematic data with age
determinations;

  * •

direct-imaging studies, including forthcoming JWST observations. Even in the
case of a nondetection, the tool can be useful to turn contrast limit into
mass limits, paving the way to a systematic assessment of the demographics of
direct-imaged exoplanets and brown dwarfs.

###### Acknowledgements.

We are grateful to Dr. Núria Miret-Roig and Dr. Raffaele Gratton for their
precious advice. This research has made use of the VizieR catalogue access
tool, CDS, Strasbourg, France (DOI : 10.26093/cds/vizier). The original
description of the VizieR service was published in 2000, A&AS 143, 23. This
research has made use of the SIMBAD database, operated at CDS, Strasbourg,
France. This work has made use of data from the European Space Agency (ESA)
mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data
Processing and Analysis Consortium (DPAC,
http://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement. This publication makes use
of data products from the Two Micron All Sky Survey, which is a joint project
of the University of Massachusetts and the Infrared Processing and Analysis
Center/California Institute of Technology, funded by the National Aeronautics
and Space Administration and the National Science Foundation. VS acknowledges
the support of PRIN-INAF 2019 Planetary Systems At Early Ages (PLATEA). For
the purpose of open access, the authors have agreed to apply a Creative
Commons Attribution (CC BY) licence to any Author Accepted Manuscript version
arising from this submission.

## References

  * Allard (2016) Allard, F. 2016, in SF2A-2016: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, ed. C. Reylé, J. Richard, L. Cambrésy, M. Deleuil, E. Pécontal, L. Tresse, & I. Vauglin, 223–227
  * Allard et al. (2001) Allard, F., Hauschildt, P. H., Alexander, D. R., Tamanai, A., & Schweitzer, A. 2001, ApJ, 556, 357
  * Amard et al. (2019) Amard, L., Palacios, A., Charbonnel, C., et al. 2019, A&A, 631, A77
  * Angus et al. (2019a) Angus, R., Morton, T., & Foreman-Mackey, D. 2019a, The Journal of Open Source Software, 4, 1469
  * Angus et al. (2019b) Angus, R., Morton, T. D., Foreman-Mackey, D., et al. 2019b, AJ, 158, 173
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Baines et al. (2012) Baines, E. K., White, R. J., Huber, D., et al. 2012, ApJ, 761, 57
  * Baraffe et al. (2015) Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42
  * Barnes (2007) Barnes, S. A. 2007, ApJ, 669, 1167
  * Barrado (2016) Barrado, D. 2016, in EAS Publications Series, Vol. 80-81, EAS Publications Series, 115–175
  * Bell et al. (2015) Bell, C. P. M., Mamajek, E. E., & Naylor, T. 2015, MNRAS, 454, 593
  * Bonavita et al. (2022) Bonavita, M., Fontanive, C., Gratton, R., et al. 2022, arXiv e-prints, arXiv:2205.02213
  * Bonavita et al. (2021) Bonavita, M., Gratton, R., Desidera, S., et al. 2021, arXiv e-prints, arXiv:2103.13706
  * Boyajian et al. (2012a) Boyajian, T. S., McAlister, H. A., van Belle, G., et al. 2012a, ApJ, 746, 101
  * Boyajian et al. (2012b) Boyajian, T. S., von Braun, K., van Belle, G., et al. 2012b, ApJ, 757, 112
  * Brandt et al. (2021a) Brandt, G. M., Brandt, T. D., Dupuy, T. J., Li, Y., & Michalik, D. 2021a, AJ, 161, 179
  * Brandt et al. (2021b) Brandt, G. M., Brandt, T. D., Dupuy, T. J., Michalik, D., & Marleau, G.-D. 2021b, ApJ, 915, L16
  * Burrows et al. (1997) Burrows, A., Marley, M., Hubbard, W. B., et al. 1997, ApJ, 491, 856
  * Cantat-Gaudin et al. (2018) Cantat-Gaudin, T., Jordi, C., Vallenari, A., et al. 2018, A&A, 618, A93
  * Cantat-Gaudin et al. (2019a) Cantat-Gaudin, T., Jordi, C., Wright, N. J., et al. 2019a, A&A, 626, A17
  * Cantat-Gaudin et al. (2019b) Cantat-Gaudin, T., Krone-Martins, A., Sedaghat, N., et al. 2019b, A&A, 624, A126
  * Chambers et al. (2016) Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv e-prints, arXiv:1612.05560
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102
  * Cutri et al. (2021) Cutri, R. M., Wright, E. L., Conrow, T., et al. 2021, VizieR Online Data Catalog, II/328
  * Damiani et al. (2019) Damiani, F., Prisinzano, L., Pillitteri, I., Micela, G., & Sciortino, S. 2019, A&A, 623, A112
  * de Zeeuw et al. (1999) de Zeeuw, P. T., Hoogerwerf, R., de Bruijne, J. H. J., Brown, A. G. A., & Blaauw, A. 1999, AJ, 117, 354
  * Dotter (2016) Dotter, A. 2016, ApJS, 222, 8
  * Epstein & Pinsonneault (2014) Epstein, C. R. & Pinsonneault, M. H. 2014, ApJ, 780, 159
  * Evans et al. (2018) Evans, D. W., Riello, M., De Angeli, F., et al. 2018, A&A, 616, A4
  * Feiden (2016a) Feiden, G. A. 2016a, A&A, 593, A99
  * Feiden (2016b) Feiden, G. A. 2016b, in Young Stars & Planets Near the Sun, ed. J. H. Kastner, B. Stelzer, & S. A. Metchev, Vol. 314, 79–84
  * Fouesneau et al. (2022) Fouesneau, M., Andrae, R., Dharmawardena, T., et al. 2022, arXiv e-prints, arXiv:2201.03252
  * Gagné et al. (2018) Gagné, J., Mamajek, E. E., Malo, L., et al. 2018, ApJ, 856, 23
  * Gaia Collaboration et al. (2018a) Gaia Collaboration, Babusiaux, C., van Leeuwen, F., et al. 2018a, A&A, 616, A10
  * Gaia Collaboration et al. (2018b) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018b, A&A, 616, A1
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1
  * Gardner et al. (2006) Gardner, J. P., Mather, J. C., Clampin, M., et al. 2006, Space Sci. Rev., 123, 485
  * Gordon et al. (2021) Gordon, K. D., Misselt, K. A., Bouwman, J., et al. 2021, ApJ, 916, 33
  * Goździewski & Migaszewski (2020) Goździewski, K. & Migaszewski, C. 2020, ApJ, 902, L40
  * Haemmerlé et al. (2019) Haemmerlé, L., Eggenberger, P., Ekström, S., et al. 2019, A&A, 624, A137
  * Hauschildt et al. (1999) Hauschildt, P. H., Allard, F., Ferguson, J., Baron, E., & Alexander, D. R. 1999, ApJ, 525, 871
  * Helmi et al. (2018) Helmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Nature, 563, 85
  * Howes et al. (2019) Howes, L. M., Lindegren, L., Feltzing, S., Church, R. P., & Bensby, T. 2019, A&A, 622, A27
  * Huber et al. (2012) Huber, D., Ireland, M. J., Bedding, T. R., et al. 2012, ApJ, 760, 32
  * Janson et al. (2021) Janson, M., Squicciarini, V., Delorme, P., et al. 2021, A&A, 646, A164
  * Kerr et al. (2021) Kerr, R. M. P., Rizzuto, A. C., Kraus, A. L., & Offner, S. S. R. 2021, ApJ, 917, 23
  * Kirkpatrick et al. (2019) Kirkpatrick, J. D., Martin, E. C., Smart, R. L., et al. 2019, ApJS, 240, 19
  * Krolikowski et al. (2021) Krolikowski, D. M., Kraus, A. L., & Rizzuto, A. C. 2021, AJ, 162, 110
  * Kuhn et al. (2019) Kuhn, M. A., Hillenbrand, L. A., Sills, A., Feigelson, E. D., & Getman, K. V. 2019, ApJ, 870, 32
  * Lallement et al. (2019) Lallement, R., Babusiaux, C., Vergely, J. L., et al. 2019, A&A, 625, A135
  * Leike et al. (2020) Leike, R. H., Glatzle, M., & Enßlin, T. A. 2020, A&A, 639, A138
  * Luhman & Esplin (2020) Luhman, K. L. & Esplin, T. L. 2020, AJ, 160, 44
  * Mackereth et al. (2021) Mackereth, J. T., Miglio, A., Elsworth, Y., et al. 2021, MNRAS, 502, 1947
  * Mamajek & Hillenbrand (2008) Mamajek, E. E. & Hillenbrand, L. A. 2008, ApJ, 687, 1264
  * Marigo et al. (2017) Marigo, P., Girardi, L., Bressan, A., et al. 2017, ApJ, 835, 77
  * Marleau & Cumming (2014) Marleau, G. D. & Cumming, A. 2014, MNRAS, 437, 1378
  * Marley et al. (2021) Marley, M. S., Saumon, D., Visscher, C., et al. 2021, ApJ, 920, 85
  * Marois et al. (2008) Marois, C., Macintosh, B., Barman, T., et al. 2008, Science, 322, 1348
  * Marrese et al. (2019) Marrese, P. M., Marinoni, S., Fabrizio, M., & Altavilla, G. 2019, A&A, 621, A144
  * Martig et al. (2016) Martig, M., Fouesneau, M., Rix, H.-W., et al. 2016, MNRAS, 456, 3655
  * Mints (2020) Mints, A. 2020, arXiv e-prints, arXiv:2012.09690
  * Miret-Roig et al. (2022) Miret-Roig, N., Bouy, H., Raymond, S. N., et al. 2022, Nature Astronomy, 6, 89
  * Morton (2015) Morton, T. D. 2015, isochrones: Stellar model grid package, Astrophysics Source Code Library, record ascl:1503.010
  * Nielsen et al. (2019) Nielsen, E. L., De Rosa, R. J., Macintosh, B., et al. 2019, AJ, 158, 13
  * Ochsenbein et al. (2000) Ochsenbein, F., Bauer, P., & Marcout, J. 2000, A&AS, 143, 23
  * Pang et al. (2022) Pang, X., Tang, S.-Y., Li, Y., et al. 2022, arXiv e-prints, arXiv:2204.06000
  * Pecaut & Mamajek (2013) Pecaut, M. J. & Mamajek, E. E. 2013, ApJS, 208, 9
  * Pecaut & Mamajek (2016) Pecaut, M. J. & Mamajek, E. E. 2016, MNRAS, 461, 794
  * Phillips et al. (2020) Phillips, M. W., Tremblin, P., Baraffe, I., et al. 2020, A&A, 637, A38
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A1
  * Prada Moroni et al. (2016) Prada Moroni, P. G., Valle, G., Dell’Omodarme, M., & Degl’Innocenti, S. 2016, Astronomische Nachrichten, 337, 819
  * Riello et al. (2021) Riello, M., De Angeli, F., Evans, D. W., et al. 2021, A&A, 649, A3
  * Rodrigo & Solano (2020) Rodrigo, C. & Solano, E. 2020, in XIV.0 Scientific Meeting (virtual) of the Spanish Astronomical Society, 182
  * Rodrigo et al. (2012) Rodrigo, C., Solano, E., & Bayo, A. 2012, SVO Filter Profile Service Version 1.0, IVOA Working Draft 15 October 2012
  * Rodrigues et al. (2017) Rodrigues, T. S., Bossini, D., Miglio, A., et al. 2017, MNRAS, 467, 1433
  * Rodrigues et al. (2014) Rodrigues, T. S., Girardi, L., Miglio, A., et al. 2014, MNRAS, 445, 2758
  * Sepulveda & Bowler (2022) Sepulveda, A. G. & Bowler, B. P. 2022, AJ, 163, 52
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163
  * Soderblom (2010) Soderblom, D. R. 2010, ARA&A, 48, 581
  * Soderblom et al. (2014) Soderblom, D. R., Hillenbrand, L. A., Jeffries, R. D., Mamajek, E. E., & Naylor, T. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 219
  * Somers et al. (2020) Somers, G., Cao, L., & Pinsonneault, M. H. 2020, ApJ, 891, 29
  * Spada et al. (2017) Spada, F., Demarque, P., Kim, Y. C., Boyajian, T. S., & Brewer, J. M. 2017, ApJ, 838, 161
  * Spiegel & Burrows (2012) Spiegel, D. S. & Burrows, A. 2012, ApJ, 745, 174
  * Squicciarini et al. (2021) Squicciarini, V., Gratton, R., Bonavita, M., & Mesa, D. 2021, MNRAS, 507, 1381
  * Squicciarini et al. (2022) Squicciarini, V., Gratton, R., Janson, M., et al. 2022, arXiv e-prints, arXiv:2205.02279
  * Steinmetz et al. (2020) Steinmetz, M., Guiglion, G., McMillan, P. J., et al. 2020, AJ, 160, 83
  * Tokunaga et al. (2002) Tokunaga, A. T., Simons, D. A., & Vacca, W. D. 2002, PASP, 114, 180
  * Verde (2010) Verde, L. 2010, in Lecture Notes in Physics, Berlin Springer Verlag, ed. G. Wolschin, Vol. 800, 147–177
  * Vigan et al. (2021) Vigan, A., Fontanive, C., Meyer, M., et al. 2021, A&A, 651, A72
  * Vines & Jenkins (2022) Vines, J. I. & Jenkins, J. S. 2022, MNRAS[arXiv:2204.03769]
  * Wang et al. (2018) Wang, J. J., Graham, J. R., Dawson, R., et al. 2018, AJ, 156, 192
  * Wang & Chen (2019) Wang, S. & Chen, X. 2019, ApJ, 877, 116
  * Wenger et al. (2000) Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9
  * Wilner et al. (2018) Wilner, D. J., MacGregor, M. A., Andrews, S. M., et al. 2018, ApJ, 855, 56
  * Zhang et al. (2021) Zhang, J., Zhao, J., Oswalt, T. D., et al. 2021, ApJ, 908, 207
  * Zurlo et al. (2022) Zurlo, A., Goździewski, K., & Lazzoni, C. 2022, submitted on A&A
  * Zurlo et al. (2016) Zurlo, A., Vigan, A., Galicher, R., et al. 2016, A&A, 587, A57