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.

Getting started#

We start by importing the required Python modules.

import urllib.request
import species

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

Initiating species v0.5.1... [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
Creating species_database.hdf5... [DONE]
Creating data folder... [DONE]
<species.core.init.SpeciesInit at 0x14dcdab00>

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

('gqlupb_sinfoni_j.dat', <http.client.HTTPMessage at 0x10e72b160>)

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.

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.

database.add_companion('GQ Lup B')
Downloading Vega spectrum (270 kB)... [DONE]
Adding Vega spectrum... [DONE]
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]
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.

database.add_object('GQ Lup B',
                    spectrum={'SINFONI': ('gqlupb_sinfoni_j.dat', None, 2500.)},
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.

line = species.EmissionLine(object_name='GQ Lup B',
                            wavel_range=(1.26, 1.32))

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

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

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.

Fitting continuum... [DONE]
Plotting continuum fit... [DONE]

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.

line.integrate_flux(wavel_int=(1.280, 1.285),
Plotting integrated line... [DONE]
Mean wavelength (nm): 1282.09 +/- 0.02
FWHM (km s-1): 202.63 +/- 8.61
Radial velocity (km s-1): 20.8 +/- 5.6
Line flux (W m-2): 1.43e-18 +/- 2.78e-20
Line luminosity (Lsun): 1.07e-06 +/- 2.06e-08
(1.4349177338040668e-18, 2.775071343210792e-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).

                  bounds={'gauss_amplitude': (0., 2e-15)},
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: 53479
[ultranest] Writing samples and results to disk ...
[ultranest] Writing samples and results to disk ... done
[ultranest]   logZ = -1.325e+07 +- 0.1514
[ultranest] Effective samples strategy satisfied (ESS = 1791.3, need >400)
[ultranest] Posterior uncertainty strategy is satisfied (KL: 0.46+-0.07 nat, need <0.50 nat)
[ultranest] Evidency uncertainty strategy is satisfied (dlogz=0.15, need <0.5)
[ultranest]   logZ error budget: single: 0.23 bs:0.15 tail:0.01 total:0.15 required:<0.50
[ultranest] done iterating.
Log-evidence = -13250634.12 +/- 0.33
Best-fit parameters (mean +/- std):
   - gauss_amplitude = 1.23e-15 +/- 1.08e-17
   - gauss_mean = 1.28e+00 +/- 5.89e-06
   - gauss_sigma = 4.29e-04 +/- 5.91e-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... [DONE]
Log-evidence = -13250634.12 +/- 0.33
Integrated autocorrelation time:
   - gauss_amplitude: fixed
   - gauss_mean: 1.03
   - gauss_sigma: 1.00
   - gauss_fwhm: 1.00
   - line_flux: fixed
   - line_luminosity: 1.00
   - line_eq_width: 1.00
   - line_vrad: 1.03
Plotting best-fit line profile... [DONE]

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.

                       title_fmt=['.2e', '.3f', '.3f', '.1f', '.2e', '.2e', '.2f', '.0f'],
                       offset=(-0.4, -0.35),
Median sample:
   - gauss_amplitude = 1.23e-15
   - gauss_mean = 1.28e+00
   - gauss_sigma = 4.29e-04
   - gauss_fwhm = 2.36e+02
   - line_flux = 1.32e-18
   - line_luminosity = 9.81e-07
   - line_eq_width = -3.72e+00
   - line_vrad = 1.50e+01
   - distance = 1.54e+02
Plotting the posterior... [DONE]

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