How to use ARIADNE

Stellar information setup

To use ARIADNE start by setting up the stellar information, this is done by importing the Star module.

from astroARIADNE.star import Star

After importing, a star has to be defined.

Stars are defined in ARIADNE by their RA and DEC in degrees, a name, and optionally the Gaia DR3 source id, for example:

starname = 'WASP-19'
ra = 148.41676021592826
dec = -45.65910531582427
gaia_id = 5411736896952029568

s = Star(starname, ra, dec, g_id=gaia_id)

The starname is used purely for user identification later on, and the gaia_id is provided to make sure the automatic photometry retrieval collects the correct magnitudes, otherwise ARIADNE will try and get the gaia_id by itself using a cone search centered around the RA and DEC.

Executing the previous block will start the photometry and stellar parameter retrieval routine. ARIADNE will query Gaia DR2 for an estimate on the temperature, radius, parallax and luminosity for display as preliminar information, as it’s not used during the fit, and prints them along with its TIC, KIC IDs if any of those exist, its Gaia DR3 ID, and maximum line-of-sight extinction Av:

                        Gaia DR2 ID : 5411736896952029568
			Gaia Effective temperature : 5458.333 +/- 109.667
			Gaia Stellar radius : 1.001 +/- 0.195
			Gaia Stellar Luminosity : 0.802 +/- 0.009
			Gaia Parallax : 3.751 +/- 0.024
			Bailer-Jones distance : 265.144 +/- 0.620
			Maximum Av : 0.581

If you already know any of those values, you can override the search for them by providing them in the Star constructor with their respective uncertainties. Likewise if you already have the magnitudes and wish to override the on-line search, you can provide a dictionary where the keys are the filters and values are the mag, mag_err tuples.

If you want to check the retrieved magnitudes you can call the print_mags method from Star:

s.print_mags()

This will print the filters used, magnitudes and uncertainties. For WASP-19 this would look like this:

		       Filter       	Magnitude	Uncertainty
		--------------------	---------	-----------
		    SkyMapper_u     	 14.1050 	  0.0070
		    SkyMapper_v     	 13.7370 	  0.0060
		  GROUND_JOHNSON_B  	 16.7920 	  0.1820
		    SkyMapper_g     	 12.4320 	  0.0030
		    GaiaDR2v2_BP    	 12.5227 	  0.0017
		  GROUND_JOHNSON_V  	 16.0100 	  0.0000
		    SkyMapper_r     	 12.0630 	  0.0050
		    GaiaDR2v2_G     	 12.1088 	  0.0005
		    SkyMapper_i     	 11.8740 	  0.0050
		    GaiaDR2v2_RP    	 11.5532 	  0.0014
		    SkyMapper_z     	 11.8610 	  0.0080
		      2MASS_J       	 10.9110 	  0.0260
		      2MASS_H       	 10.6020 	  0.0220
		      2MASS_Ks      	 10.4810 	  0.0230
		    WISE_RSR_W1     	 10.4360 	  0.0230
		    WISE_RSR_W2     	 10.4940 	  0.0200

Note: ARIADNE automatically prints and saves the used magnitudes and filters to a file.

The way the photometry retrieval works is that Gaia DR2 crossmatch catalogs are queried for the Gaia ID, these crossmatch catalogs exist for ALL-WISE, APASS, Pan-STARRS1, SDSS, 2MASS and Tycho-2, so finding photometry relies on these crossmatches. For example, if we were to analyze NGTS-6, there are Pan-STARRS1 photometry which ARIADNE couldn’t find due to the Pan-STARRS1 source not being identified in the Gaia DR2 crossmatch, in this case if you wanted to add that photometry manually, you can do so by using the add_mag method from Star, for example, if you wanted to add the PS1_r mag to our Star object you would do:

s.add_mag(13.751, 0.032, 'PS1_r')

If for whatever reason ARIADNE found a bad photometry point, and you needed to remove it, you can invoke the remove_mag method. For example, you wanted to remove the TESS magnitude due to it being from a blended source, you can just run

s.remove_mag('NGTS')

In the specific example of WASP-19, we see that GROUND_JOHNSON_B and GROUND_JOHNSON_V are likely not the correct photometry. Instead the correct ones are 13.054 and 12.248, respectively. We can correct this mistake:

s.remove_mag('GROUND_JOHNSON_B')
s.remove_mag('GROUND_JOHNSON_V')
s.add_mag(13.054, 0.048, 'GROUND_JOHNSON_B')
s.add_mag(12.248, 0.069, 'GROUND_JOHNSON_V')

And a new call to s.print_mags() would yield:

		       Filter       	Magnitude	Uncertainty
		--------------------	---------	-----------
		    SkyMapper_u     	 14.1050 	  0.0070
		    SkyMapper_v     	 13.7370 	  0.0060
		  GROUND_JOHNSON_B  	 13.0540 	  0.0480
		    SkyMapper_g     	 12.4320 	  0.0030
		    GaiaDR2v2_BP    	 12.5227 	  0.0017
		  GROUND_JOHNSON_V  	 12.2480 	  0.0690
		    SkyMapper_r     	 12.0630 	  0.0050
		    GaiaDR2v2_G     	 12.1088 	  0.0005
		    SkyMapper_i     	 11.8740 	  0.0050
		    GaiaDR2v2_RP    	 11.5532 	  0.0014
		    SkyMapper_z     	 11.8610 	  0.0080
		      2MASS_J       	 10.9110 	  0.0260
		      2MASS_H       	 10.6020 	  0.0220
		      2MASS_Ks      	 10.4810 	  0.0230
		    WISE_RSR_W1     	 10.4360 	  0.0230
		    WISE_RSR_W2     	 10.4940 	  0.0200

A list of allowed filters can be found here

Interstellar extinction

ARIADNE has an incorporated prior for the interstellar extinction in the Visual band, \(A_{\rm V}\) which consists of a uniform prior between 0 and the maximum line-of-sight value provided by the SFD dust maps. This, however, can be changed either by a custom prior (see Fitter setup) or by changing the dustmap used. We provide following dustmaps:

• SFD (2011)

• Planck Collaboration (2013)

• Planck Collaboration (2016; GNILC)

• Lenz, Hensley & Doré (2017)

• Bayestar (2019)

These maps are all implemented through the dustmaps package and need to be downloaded. Instructions to download the dustmaps can be found in its documentation.

To change the dustmap you need to provide the dustmap parameter to the Star constructor, for example:

ra = 75.795
dec = -30.399
starname = 'NGTS-6'
gaia_id = 4875693023844840448

s = Star(starname, ra, dec, g_id=gaia_id, dustmap='Bayestar')

This concludes the stellar setup and now we’re ready to set up the parameters for the fitting routine.

Fitter setup

In this section we’ll detail how to set up the fitter for the Bayesian Model Averaging (BMA) mode of ARIADNE. For single models the procedure is very similar.

First, import the fitter from ARIADNE

from astroARIADNE.fitter import Fitter

There are several configuration parameters we have to setup, the first one is the output folder where we want ARIADNE to output the fitting files and results, next we have to select the fitting engine ('dynesty'), number of live points to use, evidence tolerance threshold, bounding method, sampling method, threads, and whether to use the dynamic nested sampler. After selecting all of those, we need to select the models we want to use and finally, we feed them all to the fitter:

Deprecation: 'dynesty' is the engine. The legacy 'multinest' engine is deprecated and will be removed entirely in v2.0 — it only ever supported single-grid fits (not BMA), and dynesty is now faster besides. Use 'dynesty'.

out_folder = 'your folder here'

engine = 'dynesty'
nlive = 500
dlogz = 0.5
bound = 'multi'
sample = 'rwalk'
threads = 4
dynamic = False

setup = [engine, nlive, dlogz, bound, sample, threads, dynamic]

# Feel free to uncomment any unneeded/unwanted models
models = [
	'phoenix',
	'btsettl',
	'btnextgen',
	'btcond',
	'kurucz',
	'ck04',
	'bosz'      # MARCS+ATLAS9/Synspec, 2800-16000 K (BMA-capable)
]

f = Fitter()
f.star = s
f.setup = setup
f.av_law = 'fitzpatrick'
f.out_folder = out_folder
f.bma = True
f.models = models
f.n_samples = 100000

Note: While you can always select all 6 models, ARIADNE has an internal filter put in place in order to avoid having the user unintentionally bias the results. For stars with Teff > 4000 K BT-Settl, BT-NextGen and BT-Cond are identical and thus only BT-Settl is used, even if the three are selected. On the other hand, Kurucz and Castelli & Kurucz are known to work poorly on stars with Teff < 4000 K, thus they aren’t used in that regime.

Additional grids. Three more grids ship with ARIADNE:

grid	regime	Teff (K)	use
bosz	FGK + warm	2800–16000	BMA member (independent MARCS+ATLAS9/Synspec physics)
sphinx	M dwarfs	2000–4000	standalone
tlusty	hot O/B (NLTE)	15000–55000	standalone

bosz can be added to the models list like any other (it is temperature-gated to its validity range). sphinx and tlusty are meant to be run standalone (see below); tlusty uses a wide uniform Teff prior because the FGK population prior cannot reach the hot-star regime.

Running a single model (standalone)

To fit a single grid instead of doing BMA, set f.bma = False and pick the grid with f.grid, then call f.fit():

f = Fitter()
f.star = s
f.setup = ['dynesty', 500, 0.1, 'multi', 'rwalk', 4, False]
f.av_law = 'fitzpatrick'
f.out_folder = out_folder
f.bma = False
f.grid = 'phoenix'      # or 'bosz', 'sphinx', 'tlusty', 'coelho', ...
f.prior_setup = {
	'teff': ('default'), 'logg': ('default'), 'z': ('default'),
	'dist': ('default'), 'rad': ('default'), 'Av': ('default'),
}
f.initialize()
f.fit()                 # single-grid fit (not fit_bma)

This is the right mode for the standalone grids: sphinx (M dwarfs), tlusty (hot stars), and coelho (optical-only). The output goes to the same files as a BMA run minus the model-averaging products.

Grid wavelength coverage. Some grids do not span the full set of filters ARIADNE supports. The Coelho 2014 grid, in particular, only covers the optical (≲ 0.9 µm); its infrared photometry (2MASS, WISE, Spitzer, Herschel) was extrapolated when the grid was built and is unreliable, especially for cool stars. For this reason Coelho is not part of the default BMA set. If you use it as a standalone grid (f.bma = False; f.grid = 'coelho'), ARIADNE automatically excludes any photometry redder than the grid’s coverage from the fit (a warning lists the dropped bands). Coverage limits are declared in grid_wave_coverage in config.py. The exclusion is not applied in BMA mode, where every grid must be fit to the same bands for the evidences to be comparable — so Coelho should only be used standalone.

Performance & tuning knobs. A few optional attributes control runtime:

• f.walks (default 25): number of rwalk MCMC steps per proposal. Can also be passed as an 8th element of setup. Lower is faster but mixes less.

• f.n_grid_jobs (default 1): number of BMA model grids to fit concurrently, one core each. Set to len(models) to run them all at once. CAUTION: this spawns that many processes — pick a value your machine’s cores and memory can handle.

• f.isochrone_dlogz (default 0.01): evidence tolerance for the MIST isochrone age/mass fit (the dominant cost of a BMA run, parallelised across threads). Raising it trades age/mass precision for speed.

We allow the use of four different extinction laws:

• fitzpatrick

• cardelli

• odonnell

• calzetti

The next step is setting up the priors to use:

f.prior_setup = {
	'teff': ('default'),
	'logg': ('default'),
	'z': ('default'),
	'dist': ('default'),
	'rad': ('default'),
	'Av': ('default')
}

A quick explanation on the priors:

The default prior for Teff is an empirical prior drawn from the RAVE survey temperatures distribution, the distance prior is drawn from the Bailer-Jones distance estimate from Gaia EDR3, and the radius has a flat prior ranging from 0.5 to 20 R\(_\odot\). The default prior for the metallicity z and log g are also their respective distributions from the RAVE survey, the default prior for Av is a flat prior that ranges from 0 to the maximum of line-of-sight as per the SFD map, finally the excess noise parameters all have gaussian priors centered around their respective uncertainties.

We offer customization on the priors as well, those are listed in the following table.

Prior	Hyperparameters
Fixed	value
Normal	mean, std
TruncNorm	mean, std, lower_lim, uppern_lim
Uniform	ini, end
RAVE (log g only)	—
Default	—

So if you knew (from a spectroscopic analysis, for example) that the effective temperature is 5600 +/- 100 and the metallicity is [Fe/H] = 0.09 +/- 0.05 and you wanted to use them as priors, and the star is nearby (< 70 pc), so you wanted to fix Av to 0, your prior dictionary should look like this:

f.prior_setup = {
	'teff': ('normal', 5600, 100),
	'logg': ('default'),
	'z': ('normal', 0.09, 0.05),
	'dist': ('default'),
	'rad': ('default'),
	'Av': ('fixed', 0)
}

Though leaving everything at default usually works well enough.

f.prior_setup = {
    'teff': ('default'),
    'logg': ('default'),
    'z': ('default'),
    'dist': ('default'),
    'rad': ('default'),
    'Av': ('default')
}

After having set up everything we can finally initialize the fitter and start fitting

f.initialize()
f.fit_bma()

Now we wait for our results!

Visualization

After the fitting has finished, we need to visualize our results. ARIADNE includes a plotter object to do just that! We first star by importing the plotter:

from astroARIADNE.plotter import SEDPlotter

The setup for the plotter is already made for you, but if you really want to change them, instructions on how to change it can be found here

If you want SED model plots, make sure you have either the spectra cache (recommended) or the full model grids installed — see Model Spectra for SED Plotting above. The spectra cache is used automatically if present; the full model grids are used via the ARIADNE_MODELS environment variable as a fallback.

Now we only need to specify the results file location and the output folder for the plots!

in_file = out_folder + 'BMA_out.pkl'
plots_out_folder = 'your plots folder here'

Now we instantiate the plotter and call the desired plotting methods! We offer 5 different plots:

• A RAW SED plot

• A SED plot with the model and synthetic photometry

• A corner plot

• An HR diagram taken from MIST isochrones

• Histograms showing the parameter distributions for each model.

artist = SEDPlotter(in_file, plots_out_folder)
artist.plot_SED_no_model()
artist.plot_SED()
artist.plot_bma_hist()
artist.plot_bma_HR(10)
artist.plot_corner()

The number given to plot_bma_HR is the number of extra tracks you want to plot, drawn randomly from the posterior distribution.

If you’re iterating through lots of stars you can call the SEDPlotter clean method to clear opened figures with artist.clean()

If you don’t have either the spectra cache or the full model grids, then the plot_SED method will be skipped.

An example usage file is provided in the repository called test_bma.py demonstrating the recommended BMA (Bayesian Model Averaging) approach.