Emission line analysis#

In this tutorial we will analyze the Pa\(\beta\) emission line in the \(J\) band spectrum of GQ Lup B. This hydrogen emission line is related to the accretion processes in the circumsubstellar disk around this young directly-imaged brown dwarf. The line parameters and uncertainties will be estimated from the continuum-subtracted line profile. For hydrogen lines, the line luminosity will be converted into an accretion luminosity with the coefficients from Aoyama et al. (2021) and Marleau & Aoyama (in prep.).

Getting started#

We start by importing the required Python modules.

[1]:
import urllib.request
import species

Next, we initiate the species workflow by running the SpeciesInit class.

[2]:
species.SpeciesInit()
Initiating species v0.5.4... [DONE]
Creating species_config.ini... [DONE]
Database: /Users/tomasstolker/applications/species/docs/tutorials/species_database.hdf5
Data folder: /Users/tomasstolker/applications/species/docs/tutorials/data
Working folder: /Users/tomasstolker/applications/species/docs/tutorials
Grid interpolation method: linear
Creating species_database.hdf5... [DONE]
Creating data folder... [DONE]
[2]:
<species.core.init.SpeciesInit at 0x151d8d2a0>

We will also download the VLT/SINFONI \(J\) band spectrum of GQ Lup B that has been published by Seifahrt et al. (2007).

[3]:
urllib.request.urlretrieve('https://home.strw.leidenuniv.nl/~stolker/species/gqlupb_sinfoni_j.dat',
                           'gqlupb_sinfoni_j.dat')
[3]:
('gqlupb_sinfoni_j.dat', <http.client.HTTPMessage at 0x151d8d690>)

Adding data of GQ Lup B to the database#

To add data to the HDF5 database, we first create an instance of Database. This object provides read and write access to the HDF5 file.

[4]:
database = species.Database()

Next, we will add the available magnitudes and parallaxes of GQ Lup B to the database (only the parallax is required later on for calculating the line luminosity) by using the add_companion method of Database. This method will add the magnitudes and filter profiles, but it also downloads a flux-calibrated spectrum of Vega and converts the magnitudes into fluxes.

[5]:
database.add_companion('GQ Lup B')
Adding filter: HST/WFPC2-PC.F606W... [DONE]
Adding filter: HST/WFPC2-PC.F814W... [DONE]
Adding filter: HST/NICMOS2.F171M... [DONE]
Adding filter: HST/NICMOS2.F190N... [DONE]
Adding filter: HST/NICMOS2.F215N... [DONE]
Adding filter: Magellan/VisAO.ip... [DONE]
Adding filter: Magellan/VisAO.zp... [DONE]
Adding filter: Magellan/VisAO.Ys... [DONE]
Adding filter: MKO/NSFCam.H... [DONE]
Adding filter: Paranal/NACO.NB405... [DONE]
Adding filter: Paranal/NACO.Mp... [DONE]
Adding filter: Paranal/NACO.Ks... [DONE]
Adding filter: Subaru/CIAO.CH4s... [DONE]
Adding filter: Subaru/CIAO.K... [DONE]
Adding filter: Subaru/CIAO.Lp... [DONE]
Downloading Vega spectrum (270 kB)... [DONE]
Adding Vega spectrum... [DONE]
Adding object: GQ Lup B
   - Parallax (mas) = 6.49 +/- 0.03
   - HST/WFPC2-PC.F606W:
      - Apparent magnitude = 19.19 +/- 0.07
      - Flux (W m-2 um-1) = 5.94e-16 +/- 3.83e-17
   - HST/WFPC2-PC.F814W:
      - Apparent magnitude = 17.67 +/- 0.05
      - Flux (W m-2 um-1) = 1.02e-15 +/- 4.68e-17
   - HST/NICMOS2.F171M:
      - Apparent magnitude = 13.84 +/- 0.13
      - Flux (W m-2 um-1) = 2.91e-15 +/- 3.49e-16
   - HST/NICMOS2.F190N:
      - Apparent magnitude = 14.08 +/- 0.20
      - Flux (W m-2 um-1) = 1.65e-15 +/- 3.06e-16
   - HST/NICMOS2.F215N:
      - Apparent magnitude = 13.40 +/- 0.15
      - Flux (W m-2 um-1) = 1.94e-15 +/- 2.70e-16
   - Magellan/VisAO.ip:
      - Apparent magnitude = 18.89 +/- 0.24
      - Flux (W m-2 um-1) = 3.74e-16 +/- 8.33e-17
   - Magellan/VisAO.zp:
      - Apparent magnitude = 16.40 +/- 0.10
      - Flux (W m-2 um-1) = 2.31e-15 +/- 2.13e-16
   - Magellan/VisAO.Ys:
      - Apparent magnitude = 15.88 +/- 0.10
      - Flux (W m-2 um-1) = 3.08e-15 +/- 2.84e-16
   - MKO/NSFCam.H:
      - Apparent magnitude = 14.01 +/- 0.13
      - Flux (W m-2 um-1) = 3.06e-15 +/- 3.67e-16
   - Paranal/NACO.NB405:
      - Apparent magnitude = 12.29 +/- 0.06
      - Flux (W m-2 um-1) = 4.82e-16 +/- 2.66e-17
   - Paranal/NACO.Mp:
      - Apparent magnitude = 11.97 +/- 0.08
      - Flux (W m-2 um-1) = 3.53e-16 +/- 2.60e-17
   - Paranal/NACO.Ks (4 values):
      - Apparent magnitude = 13.47 +/- 0.03
      - Flux (W m-2 um-1) = 1.88e-15 +/- 5.36e-17
      - Apparent magnitude = 13.39 +/- 0.03
      - Flux (W m-2 um-1) = 2.03e-15 +/- 6.00e-17
      - Apparent magnitude = 13.50 +/- 0.05
      - Flux (W m-2 um-1) = 1.84e-15 +/- 8.47e-17
      - Apparent magnitude = 13.50 +/- 0.03
      - Flux (W m-2 um-1) = 1.83e-15 +/- 4.72e-17
   - Subaru/CIAO.CH4s:
      - Apparent magnitude = 13.76 +/- 0.26
      - Flux (W m-2 um-1) = 3.96e-15 +/- 9.58e-16
   - Subaru/CIAO.K:
      - Apparent magnitude = 13.37 +/- 0.12
      - Flux (W m-2 um-1) = 1.87e-15 +/- 2.07e-16
   - Subaru/CIAO.Lp:
      - Apparent magnitude = 12.44 +/- 0.22
      - Flux (W m-2 um-1) = 5.78e-16 +/- 1.18e-16

Now that the photometric data of GQ Lup B are stored in the database, we will add the \(J\) band spectrum by using the add_object method. Here, it is important that we use the same argument for object_name as was used as name with add_companion (i.e. GQ Lup B). The argument of spectrum contains a dictionary with the spectra. Each value is a tuple with the spectrum, covariance matrix (set to None in this case), and the spectral resolution.

[6]:
database.add_object('GQ Lup B',
                    app_mag=None,
                    flux_density=None,
                    spectrum={'SINFONI': ('gqlupb_sinfoni_j.dat', None, 2500.)},
                    deredden=None)
Adding object: GQ Lup B
   - Spectrum:
      - Database tag: SINFONI
      - Filename: gqlupb_sinfoni_j.dat
      - Data shape: (1014, 3)
      - Wavelength range (um): 1.20 - 1.35
      - Mean flux (W m-2 um-1): 3.27e-15
      - Mean error (W m-2 um-1): 2.79e-17
   - Spectral resolution:
      - SINFONI: 2500.0

Initiating the line analysis#

We will now analyze the emission line by first creating an instance of EmissionLine. The spectrum is read from the database by providing the same object name and spectrum name as were used with add_object. We will select a limited wavelength range around the Pa\(\beta\) line by setting the argument of wavel_range.

We need to provide either a name for a hydrogen line as argument of hydrogen_line or the rest wavelength as argument of lambda_rest. If both arguments are set to None, a list names of the available hydrogen lines will be shown and the input can be provided. For the available lines, the rest wavelength (in vacuum) will be automatically calculated and the accretion relation from Aoyama et al. (2021) will be used to convert the measured line luminosity into an accretion luminosity.

[7]:
line_analysis = species.EmissionLine(object_name='GQ Lup B',
                                     spec_name='SINFONI',
                                     hydrogen_line=None,
                                     lambda_rest=None,
                                     wavel_range=(1.26, 1.32))
Downloading coefficients for accretion relation (2.1 kB)... [DONE]
Adding coefficients for accretion relation (2.1 kB)...
/Users/tomasstolker/applications/species/species/analysis/emission_line.py:130: UserWarning: Provide an argument for either the 'hydrogen_line' or 'lambda_rest' parameter. Only one of the two arguments should be set to None.
  warnings.warn(
 [DONE]
Please cite Aoyama et al. (2021) when using the accretion relation in a publication
Available hydrogen lines:
['Lya', 'Lyb', 'Lyg', 'Ha', 'Hb', 'Hg', 'Hd', 'H7', 'H8', 'H9', 'H10', 'Paa', 'Pab', 'Pag', 'Pad', 'Pa8', 'Pa9', 'Pa10', 'Bra', 'Brb', 'Brg', 'Brd', 'Br9']
Please provide the name of the hydrogen line: Pab
Hydrogen line = Pab
Rest wavelength (um) = 1.2815 um

Let’s have a look at the data of the cropped \(J\) spectrum that is extracted.

[8]:
line_analysis.spectrum
[8]:
array([[1.26000500e+00, 3.28487860e-15, 1.58013906e-17],
       [1.26014996e+00, 3.30313622e-15, 1.06466256e-17],
       [1.26029503e+00, 3.32310600e-15, 1.08635306e-17],
       ...,
       [1.31959999e+00, 3.36931744e-15, 8.78280340e-18],
       [1.31974494e+00, 3.40047958e-15, 6.23790556e-17],
       [1.31989002e+00, 3.39223842e-15, 3.46253392e-17]])

Converting line luminosity into accretion luminosity#

The list_hydrogen_lines method can be used to list the available names of the hydrogen lines for which the line luminosity can be converted into an accretion luminosity. The conversion is done by using the coefficient for the accretion relation from Aoyama et al. (2021).

[9]:
_ = line_analysis.list_hydrogen_lines()
Available hydrogen lines:
['Lya', 'Lyb', 'Lyg', 'Ha', 'Hb', 'Hg', 'Hd', 'H7', 'H8', 'H9', 'H10', 'Paa', 'Pab', 'Pag', 'Pad', 'Pa8', 'Pa9', 'Pa10', 'Bra', 'Brb', 'Brg', 'Brd', 'Br9']

Since we have set the hydrogen_line argument of EmissionLine, we can use the accretion_luminosity method to manually convert a line luminosity (in \(L_\odot\)) into an accretion luminosity (in \(L_\odot\)). Let’s give that a try for the Pa\(\beta\) line with \(L_\mathrm{line} = 10^{-5}\) \(L_\odot\).

[10]:
acc_lum = line_analysis.accretion_luminosity(1e-5)
print(f"Accretion luminosity (Lsun) = {acc_lum:.2e}")
Accretion luminosity (Lsun) = 7.82e-03

Subtracting the continuum flux#

In the \(J\) band, we need to subtract the pseudo-continuum of the atmospheric emission. To do so, we make use of the subtract_continuum method of EmissionLine. The continuum is modeled with a 3rd order polynomial (after smoothing the spectrum) and we use a linear least square fitting algorithm to optimize the parameters. A plot with the best-fit polynomial and the continuum-subtracted spectrum is created, which can be inspected before continuing with the line analysis.

[11]:
line_analysis.subtract_continuum(poly_degree=3,
                        plot_filename=None)
Fitting continuum... [DONE]
Plotting continuum fit... [DONE]
../_images/tutorials_emission_line_28_1.png

Integrating the line flux#

Now that the continuum is subtracted, we can determine the line characteristics in two different ways. Here, we will measure the line properties directly from the spectrum. Later on, we will fit a Gaussian profile and adopt the sampled parameters for inferring the line properties.

We use the integrate_flux method to integrate the line flux across the specified wavelength range of wavel_int after interpolating the spectrum to \(R = 100000\). Several other characteristics are also calculated from the spectral line. The uncertainties are estimated by sampling 1000 random spectra from the errors on the fluxes.

[12]:
line_analysis.integrate_flux(wavel_int=(1.280, 1.285),
                             interp_kind='linear',
                             plot_filename=None)
Plotting integrated line... [DONE]
../_images/tutorials_emission_line_31_1.png
Mean wavelength (nm) = 1282.09 +/- 0.02
FWHM (km s-1) = 203.06 +/- 8.26
Radial velocity (km s-1) = 144.6 +/- 5.6
Line flux (W m-2) = 1.43e-18 +/- 2.68e-20
Line luminosity (Lsun) = 1.07e-06 +/- 1.99e-08
Line luminosity log10(L/Lsun) = -5.97 +/- 0.01
Inflating the uncertainty on the accretion luminosity by 0.3 dex
 to account for the the model uncertainty (see Aoyama et al. 2021)...
Accretion luminosity log10(L/Lsun) = -2.95 +/- 0.30
[12]:
(1.4349177338040668e-18, 2.6772263940095102e-20)

Modeling the Pa\(\beta\) emission line#

We will now use the fit_gaussian method of EmissionLine to fit the spectrum with a Gaussian function while mapping the posterior distributions of the 3 free parameters (amplitude, mean, and standard deviation) with the nested sampling algorithm of UltraNest.

The argument of bounds is a dictionary with the prior boundaries of the three parameters (gauss_amplitude, gauss_mean, and gauss_sigma). If any parameter is missing or bounds=None then conservative boundaries will be estimated from the spectrum. The posterior samples will be stored in the database with the tag name and a plot is created with a comparison of the data and best-fit model (i.e. the median parameter values).

[13]:
line_analysis.fit_gaussian(tag='pa-beta',
                           min_num_live_points=400,
                           bounds={'gauss_amplitude': (0., 2e-15)},
                           output='ultranest',
                           plot_filename=None,
                           show_status=False)
Creating directory for new run ultranest/run1
[ultranest] Sampling 400 live points from prior ...
/Users/tomasstolker/.pyenv/versions/3.10.0/envs/species3.10/lib/python3.10/site-packages/ultranest/store.py:195: DeprecationWarning: `np.float` is a deprecated alias for the builtin `float`. To silence this warning, use `float` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.float64` here.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
  'points', dtype=np.float,
[ultranest] Explored until L=-1e+07
[ultranest] Likelihood function evaluations: 50946
[ultranest] Writing samples and results to disk ...
[ultranest] Writing samples and results to disk ... done
[ultranest]   logZ = -1.325e+07 +- 0.19
[ultranest] Effective samples strategy satisfied (ESS = 1696.1, need >400)
[ultranest] Posterior uncertainty strategy is satisfied (KL: 0.47+-0.07 nat, need <0.50 nat)
[ultranest] Evidency uncertainty strategy is satisfied (dlogz=0.19, need <0.5)
[ultranest]   logZ error budget: single: 0.23 bs:0.19 tail:0.01 total:0.19 required:<0.50
[ultranest] done iterating.
Log-evidence = -13250633.99 +/- 0.36
Best-fit parameters (mean +/- std):
   - gauss_amplitude = 1.23e-15 +/- 1.05e-17
   - gauss_mean = 1.28e+00 +/- 5.67e-06
   - gauss_sigma = 4.29e-04 +/- 5.77e-06
Maximum likelihood sample:
   - Log-likelihood = -1.33e+07
   - gauss_amplitude = 1.23e-15
   - gauss_mean = 1.28e+00
   - gauss_sigma = 4.30e-04
Calculating line fluxes...Inflating the uncertainty on the accretion luminosity by 0.3 dex
 to account for the the model uncertainty (see Aoyama et al. 2021)...
 [DONE]
Log-evidence = -13250633.99 +/- 0.36
Integrated autocorrelation time:
   - gauss_amplitude: fixed
   - gauss_mean: 0.93
   - gauss_sigma: 0.99
   - gauss_fwhm: 0.99
   - line_flux: fixed
   - log_line_lum: 1.02
   - log_acc_lum: 1.07
   - line_eq_width: 1.00
   - line_vrad: 0.93
Plotting best-fit line profile... [DONE]
../_images/tutorials_emission_line_34_3.png

Plotting the posterior distributions#

Now that the posterior samples are available in the database, we can use the plot_posterior function for visualizing the 1D and 2D projected distributions with corner.py. The argument of tag is set to the same database tag that was used with fit_gaussian.

[14]:
species.plot_posterior(tag='pa-beta',
                       burnin=None,
                       title=None,
                       title_fmt=['.2e', '.3f', '.3f', '.1f', '.2e', '.3f', '.2f', '.2f', '.0f'],
                       offset=(-0.4, -0.35),
                       output=None)
Median sample:
   - gauss_amplitude = 1.23e-15
   - gauss_mean = 1.28e+00
   - gauss_sigma = 4.30e-04
   - gauss_fwhm = 2.37e+02
   - line_flux = 1.32e-18
   - log_line_lum = -6.01e+00
   - log_acc_lum = -2.98e+00
   - line_eq_width = -3.72e+00
   - line_vrad = 1.39e+02
   - parallax = 6.49e+00
Plotting the posterior...
../_images/tutorials_emission_line_37_1.png
 [DONE]

The first three parameters (from left to right) are the free parameters while the FWHM velocity, line flux, line luminosity, accretion luminosity, equivalent width, and radial velocity have been calculated from the posterior distributions of these three Gaussian model parameters.