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:

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)[source] [github] [bitbucket]

RJ equivalent temperature (MS15 eqn 24)

pyspeckit.spectrum.models.lte_molecule.Jnu_cgs(nu, T)[source] [github] [bitbucket]

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

pyspeckit.spectrum.models.lte_molecule.generate_fitter(model_func, name)[source] [github] [bitbucket]

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)[source] [github] [bitbucket]

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

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

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

>>> 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)[source] [github] [bitbucket]

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)[source] [github] [bitbucket]

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)[source] [github] [bitbucket]

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

Mangum & Shirley 2015 eqn 31 Ntot/Nu = Q_rot / gu exp(E_u/k Tex)

Example:
>>> 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)[source] [github] [bitbucket]

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

\[Ntot = (3 k) / (8 pi^3 nu S mu^2 R_i) (Q/g) exp(E_u/k Tex) integ(T_R/f dv)\]

Eqn 31:

\[ \begin{align}\begin{aligned}Ntot/Nu = Q_rot / gu exp(E_u/k Tex)\\-> Ntot = Nu Q_rot / gu exp(E_u/k Tex) -> Nu = N_tot g / Qrot exp(-E_u / k Tex)\end{aligned}\end{align} \]

To get Nu of an observed line, then:

\[Nu Q_rot / gu exp(E_u/k Tex) = (3 k) / (8 pi^3 nu S mu^2 R_i) (Q/g) exp(E_u/k Tex) integ(T_R/f dv)\]

This term cancels:

\[Q_rot / gu exp(E_u/k Tex)\]

Leaving:

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

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

\[Nu = (3 k c) / (8 pi^3 nu^2 S mu^2 R_i) integ(T_R/f dnu)\]

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 Aul) / (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 Aul) / (64 pi^4 nu^3))) integ(T_R/f dnu) = (3 k c 64 pi^4 nu^3) / (8 pi^3 nu^2 3 h c^3 Aul) integ(T_R/f dnu) = (8 pi nu k / (Aul c^2 h)) integ(T_R/f dnu)\]

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 / (Aul c^3 h)) integ(T_R/f dv)\]