LTE Molecule Model

There is a tool for modeling the full LTE spectrum of a molecule. It uses the JPL or CDMS databases through their astroquery interfaces.

A very simple example looks like this:

import astropy.units as u
from pyspeckit.spectrum.models.lte_molecule import get_molecular_parameters, generate_fitter, generate_model

freqs, aij, deg, EU, partfunc = get_molecular_parameters('CH3OH',
                                                         fmin=90*u.GHz,
                                                         fmax=100*u.GHz)
def modfunc(xarr, vcen, width, tex, column):
    return generate_model(xarr, vcen, width, tex, column, freqs=freqs, aij=aij,
                          deg=deg, EU=EU, partfunc=partfunc)

fitter = generate_fitter(modfunc, name="CH3OH")

Details can be found in the API documentation:

LTE Molecule Modeling Tool

Uses astroquery to obtain molecular parameters.

Equations are based on Mangum & Shirley 2015 (2015PASP..127..266M)

Module API

pyspeckit.spectrum.models.lte_molecule.Jnu(nu, T)

RJ equivalent temperature (MS15 eqn 24)

pyspeckit.spectrum.models.lte_molecule.Jnu_cgs(nu, T)

RJ equivalent temperature (MS15 eqn 24) (use cgs constants for speed)

pyspeckit.spectrum.models.lte_molecule.generate_fitter(model_func, name)

Generator for hnco fitter class

pyspeckit.spectrum.models.lte_molecule.generate_model(xarr, vcen, width, tex, column, freqs, aij, deg, EU, partfunc, background=None, tbg=2.73, get_tau=False, get_tau_sticks=False)

Model Generator

Parameters

xarr : Quantity array [Hz]

The X-axis (frequency)

vcen : Quantity [km/s]

The central velocity

width : Quantity [km/s]

The 1-sigma Gaussian line width [not FWHM]

tex : Quantity [K]

Excitation temperature

column : Quantity [cm^-2]

The total column density of the molecule. If specified as a value < 30, it will be assumed to be the log10 of the column density.

freqs : Quantity array [Hz]

The central frequency of the lines to model. Subsequent parameters, aij, deg, EU, must have the same shape as freqs.

aij : array [log s^-1]

The Einstein A coefficients, assumed to be in units of log10 of the rate in 1/s

deg : array

The transition degeneracy

EU : Quantity array [erg]

The upper state energy in ergs

partfunc : function

The partition function. Can also be specified as an array of values in Kelvin.

background : None

An additional background field to include in the radiative transfer model. Must be in Kelvin. [Presently not correctly implemented]

tbg : Quantity [K]

The background brightness temperature, defaulting to 2.73 for the CMB temperature. This background will be treated as uniform with frequency and subtracted.

get_tau : bool, default False

If specified, the optical depth in each frequency bin will be returned

get_tau_sticks : bool, default False

If specified, the optical depth in each transition will be returned

Examples

Example case to produce a model:

from pyspeckit.spectrum.models.lte_molecule import get_molecular_parameters, generate_model, generate_fitter
freqs, aij, deg, eu, partfunc = get_molecular_parameters('CH3OH')
def modfunc(xarr, vcen, width, tex, column):
    return generate_model(xarr, vcen, width, tex, column, freqs=freqs, aij=aij,
                          deg=deg, EU=EU, partfunc=partfunc)

fitter = generate_fitter(modfunc, name="CH3OH")

pyspeckit.spectrum.models.lte_molecule.get_molecular_parameters(molecule_name, tex=50, fmin=<Quantity 1. GHz>, fmax=<Quantity 1. THz>, catalog='JPL', **kwargs)

Get the molecular parameters for a molecule from the JPL or CDMS catalog

(this version should, in principle, be entirely self-consistent)

Parameters

molecule_name : string

The string name of the molecule (normal name, like CH3OH or CH3CH2OH, but it has to match the JPL catalog spec)

tex : float

Optional excitation temperature (basically checks if the partition function calculator works)

catalog : ‘JPL’ or ‘CDMS’

Which catalog to pull from

fmin : quantity with frequency units

fmax : quantity with frequency units

The minimum and maximum frequency to search over

Examples

>>> from pyspeckit.spectrum.models.lte_molecule import get_molecular_parameters
>>> freqs, aij, deg, EU, partfunc = get_molecular_parameters('CH2CHCN',
...                                                          fmin=220*u.GHz,
...                                                          fmax=222*u.GHz,
                                                            )
>>> freqs, aij, deg, EU, partfunc = get_molecular_parameters('CH3OH',
...                                                          fmin=90*u.GHz,
...                                                          fmax=100*u.GHz)

pyspeckit.spectrum.models.lte_molecule.line_tau(tex, total_column, partition_function, degeneracy, frequency, energy_upper, einstein_A=None)

Given the excitation temperature of the state, total column density of the molecule, the partition function, the degeneracy of the state, the frequency of the state, and the upper-state energy level, return the optical depth of that transition.

This is a helper function for the LTE molecule calculations. It implements the equations

\[\int \tau_\nu d\nu = \frac{c^2}{8 \pi \nu^2} A_{ij} N_u \left[ \exp\left( \frac{h \nu}{k_B T_{ex}} \right) - 1 \right]\]

or

\[\tau_\nu = \frac{c^2}{8 \pi \nu^2} A_{ij} N_u \phi_\nu \left[ \exp\left( \frac{h \nu}{k_B T_{ex}} \right) - 1 \right]\]

where

\[N_{u} = N_{tot} \frac{g_u}{Q} \exp\left(\frac{-E_u}{k_B T_{ex}} \right)\]

based on Equations 11 and 29 of Mangum & Shirley 2015 (2015PASP..127..266M)

The return value is therefore

\[\tau_\nu / \phi_\nu\]

The line profile function is, for a Gaussian, given by eqn A1:

\[\phi_\nu = \frac{1}{\sqrt{2 \pi} \sigma} \exp \left[ -\frac{(\nu-\nu_0)^2}{2 \sigma^2} \right]\]

pyspeckit.spectrum.models.lte_molecule.line_tau_cgs(tex, total_column, partition_function, degeneracy, frequency, energy_upper, einstein_A)

Given the excitation temperature of the state, total column density of the molecule, the partition function, the degeneracy of the state, the frequency of the state, and the upper-state energy level, return the optical depth of that transition.

Unlike line_tau(), this function requires inputs in CGS units.

This is a helper function for the LTE molecule calculations. It implements the equations

\[\int \tau_\nu d\nu = \frac{c^2}{8 \pi \nu^2} A_{ij} N_u \left[ \exp\left( \frac{h \nu}{k_B T_{ex}} \right) - 1 \right]\]

or

\[\tau_\nu = \frac{c^2}{8 \pi \nu^2} A_{ij} N_u \phi_\nu \left[ \exp\left( \frac{h \nu}{k_B T_{ex}} \right) - 1 \right]\]

where

\[N_{u} = N_{tot} \frac{g_u}{Q} \exp\left(\frac{-E_u}{k_B T_{ex}} \right)\]

based on Equations 11 and 29 of Mangum & Shirley 2015 (2015PASP..127..266M)

The return value is therefore

\[\tau_\nu / \phi_\nu\]

The line profile function is, for a Gaussian, given by eqn A1:

\[\phi_\nu = \frac{1}{\sqrt{2 \pi} \sigma} \exp \left[ -\frac{(\nu-\nu_0)^2}{2 \sigma^2} \right]\]

pyspeckit.spectrum.models.lte_molecule.ntot_of_nupper(nupper, eupper, tex, Q_rot, degeneracy=1)

Given an N_upper, E_upper, tex, Q_rot, and degeneracy for a single state, give N_tot

Mangum & Shirley 2015 eqn 31

\[N_{tot}/N_u = Q_{rot} / g_u \exp(E_u/k T_{ex})\]

Example:

>>> import astropy.units as u
>>> tex = 50*u.K
>>> kkms = 100*u.K*u.km/u.s
>>> from pyspeckit.spectrum.models import lte_molecule
>>> freqs, aij, deg, EU, partfunc = lte_molecule.get_molecular_parameters(molecule_name='HNCO v=0', molecule_name_jpl='HNCO', fmin=87*u.GHz, fmax=88*u.GHz)
>>> nupper = lte_molecule.nupper_of_kkms(kkms, freqs, 10**aij)
>>> ntot = lte_molecule.ntot_of_nupper(nupper, EU*u.erg, tex, Q_rot=partfunc(tex), degeneracy=deg)

pyspeckit.spectrum.models.lte_molecule.nupper_of_kkms(kkms, freq, Aul, replace_bad=None)

Mangum & Shirley 2015 eqn 82 gives, for the optically thin, Rayleigh-Jeans, negligible background approximation:

\[N_{tot} = (3 k) / (8 \pi^3 \nu S \mu^2 R_i) (Q/g) \exp(E_u/k T_{ex}) \int(T_R/f dv)\]

Eqn 31:

\[ \begin{align}\begin{aligned}N_{tot}/N_u = Q_{rot} / g_u \exp(E_u/k T_{ex})\\-> N_{tot} = N_u Q_{rot} / g_u \exp(E_u/k T_{ex}) -> Nu = N_{tot} g / Q_{rot} \exp(-E_u / k T_{ex})\end{aligned}\end{align} \]

To get Nu of an observed line, then:

\[Nu Q_{rot} / g_u \exp(E_u/k T_{ex}) = (3 k) / (8 \pi^3 \nu S \mu^2 R_i) (Q/g) \exp(E_u/k T_{ex}) \int(T_R/f dv)\]

This term cancels:

\[Q_{rot} / g_u \exp(E_u/k T_{ex})\]

Leaving:

\[N_u = (3 k) / (8 \pi^3 \nu S \mu^2 R_i) \int(T_R/f dv)\]

$int(T_R/f dv)$ is the optically thin integrated intensity in K km/s dnu/nu = dv/c [doppler eqn], so to get $int(T_R dnu)$, sub in $dv = c/nu dnu$

\[N_u = (3 k c) / (8 \pi^3 \nu^2 S \mu^2 R_i) \int(T_R/f d\nu)\]

We then need to deal with the S mu^2 R_i term. We assume R_i = 1, since we are not measuring any hyperfine transitions (R_i is the hyperfine degeneracy; eqn 75) Equation 11:

\[A_{ul} = 64 \pi^4 \nu^3 / (3 h c^3) |\mu_{ul}|^2\]

Equation 62:

\[ \begin{align}\begin{aligned}|\mu_{ul}|^2 = S \mu^2\\-> S \mu^2 = (3 h c^3 A_{ul}) / (64 \pi^4 \nu^3)\end{aligned}\end{align} \]

Plugging that in gives

\[Nu = (3 k c) / (8 \pi^3 \nu^2 ((3 h c^3 A_{ul}) / (64 \pi^4 \nu^3))) \int(T_R/f d\nu) = (3 k c 64 \pi^4 \nu^3) / (8 \pi^3 \nu^2 3 h c^3 A_{ul}) \int(T_R/f d\nu) = (8 \pi \nu k / (A_{ul} c^2 h)) \int(T_R/f d\nu)\]

which is the equation implemented below. We could also have left this in dv units by substituting du = nu/c dv:

\[= (8 \pi \nu^2 k / (A_{ul} c^3 h)) \int (T_R/f dv)\]

[edit this page on github]

[edit this page on bitbucket]