Wednesday 25 December 2013

Atomic partition functions in IDL

Many books on radiative transfer and spectroscopy introduce the partition function ($U$, PF) as a statistical weight or probability that an atom is in the ionization stage $r$ at the temperature $T$. It is defined as the sum over all excitation states ($i$) of their statistical weights ($g_i$) multiplied by an
exponential factor (cf. Rutten, 2003, p.30, Eq.289):
$$U_r \equiv \sum_i g_i \mathrm{e}^{-\frac{\chi_{r, i}}{kT}},$$ where $\chi_{r, i}$ is the excitation potential of the state $i$ in the ionization stage $r$, $k$ is the Boltzmann constant) and $T$ is the temperature.

The trouble with PF is that there is an infinite number of the energy levels. No matter how small their contribution becomes (the exponential factor is decreasing with the excitation potential), PF computed after this definition goes to infinity. Farther, that means that the probability of an atom to occupy any finite state would be zero. Since this is clearly unphysical, thus there is a need to define a reliable procedure to cut off the partition function. There are many approached described in literature. However, they are usually rather non-trivial and, therefore, computing partition function even for a limited range of temperature and for a specific atomic element can be a very tedious job. The results of these computations are often published as the coefficients of the best polynomial fits for a range of temperature values. The IDL procedure that you can find below is written to load these tables from several different sources (see below) and to interpolate them for given a temperature, an atomic specie and an ionization stage.


In thermodynamical equilibrium (TE) distribution of atoms over their bound states is set by the Boltzmann distribution.
$$ \frac{n_i}{n_j} = \frac{g_i}{g_j}\mathrm{e}^{-\chi_{ij}/kT}, $$ where $n_i$ and $n_j$ are number densities of the particles in the bound states $i$ and $j$, $g_i$ and $g_j$ are the corresponding statistical weights, $\chi_{ij}$ is the energy difference between the two states, $k$ is Boltzmann constant and $T$ is the equilibrium temperature.

The distribution of the atoms over the ionization stages in TE is set by the Saha law:
$$ \frac{N_{r+1}}{N_{r}} = \frac{1}{n_e} \frac{2 U_{r+1}}{U_{r}} \left(\frac{2 \pi m_e k T}{h^2}\right)^{3/2} \mathrm{e}^{-\chi_{r}/kT}, $$ where $N_{r+1}$ and $N_{r}$ are the number density of two successive ionization stages, $n_e$ is the electron density, $ U_{r+1}$ and $ U_{r}$ are the partition functions of the corresponding stages, while other quantities have their common meanings.

In spectroscopy the partition functions appear only in the Saha ionization equation. However, the concept of the partition function is much broader and more general. In thermodynamics, the partition function is a link between the description of the microstates and the macroscopic quantities. Here is a little reminder of that.


Partition function in statistical mechanics

In statistical mechanics, a statistically large system within a fixed volume and with a fixed number of particles in thermal equilibrium represents a (micro)canonical ensemble. The system can occupy many microstates ($s= 1, 2, ...$), each of them with a corresponding energy $E_s$. In quantum statistical mechanics, the energies $E_s$ are discrete. For such a system, a quantity called canonical partition function can be defined:
$$Z = \sum_s \mathrm{e}^{- \beta E_s}, $$
where $\beta = 1/(kT)$ is the inverse temperature. If the microstates are degerate, i.e. if there is more microstates that correspond to the same energy $E_i$, then their contributions become $g_i \exp{-\beta E_i}$ and the canonical partition function is
$$Z = \sum_i g_i \mathrm{e}^{- \beta E_i}, $$
where $g_i$ is the number of microstates $s$ with energy $E_i$.

The canonical partition function is closely related to the probability $P_s$ that the system is in the microstate $s$:
$$ P_s = \frac{1}{Z} \mathrm{e}^{-\beta E_s}, $$ so that $\sum P_s$ is normalized to 1.

The total energy of the system $E$ (commonly labeled as $U$ or $<E>$) is the sum of energy of every possible microstate multiplies by the probability of the system to be in it:
$$E = \sum_s\,P_s\, E_s = \sum_s \frac{1}{Z}\, E_s\,\mathrm{e}^{-\beta E_s}= \frac{1}{Z} \sum_s E_s\,\mathrm{e}^{-\beta E_s}.$$ The sum on the right hand side is exactly:
$$\frac{\partial Z}{\partial \beta} = \frac{\partial}{\partial \beta}\sum_s  \mathrm{e}^{-\beta E_s}  =\sum_s  \frac{\partial}{\partial \beta} \mathrm{e}^{-\beta E_s}  = - \sum_s  E_s\, \mathrm{e}^{-\beta E_s}.$$
Therefore, the total energy of the system is
\begin{equation}E = -\frac{1}{Z}\frac{\partial Z}{\partial \beta} = -\frac{\partial \ln Z}{\partial \beta}.\end{equation} The entropy of the system is defined by the probabilities $P_s$:
\begin{eqnarray}S &=&-k \sum_s P_s \ln P_s = -k \sum_s \frac{1}{Z} \mathrm{e}^{-\beta E_s} \ln\left( \frac{1}{Z} \mathrm{e}^{-\beta E_s}\right) =\nonumber\\
&=& k \sum_s \frac{1}{Z} \mathrm{e}^{-\beta E_s} \ln Z + k \sum_s \frac{1}{Z} \mathrm{e}^{-\beta E_s} \beta E_s = \nonumber\\
&=& k \frac{1}{Z} \ln Z \sum_s  \mathrm{e}^{-\beta E_s}  + k \frac{1}{Z} \beta \sum_s  E_s \mathrm{e}^{-\beta E_s}  = k(\ln Z + \beta E). \nonumber\\
&=& k\left(\ln Z - \beta \frac{\partial \ln Z}{\partial \beta} \right).\nonumber \end{eqnarray}
In addition, we substitute $\beta$ with $1/kT$:
$$S = k\left(\ln Z - \frac{1}{kT} \frac{\partial \ln Z}{\partial (kT)^{-1}}\right)= k\left( \ln Z - \frac{1}{kT} kT^2 \frac{\partial \ln Z}{\partial T}\right),$$ \begin{equation}S =  \frac{\partial (kT\ln Z)}{\partial T}.\end{equation}

The free energy that is $A = E - T\,S$ in terms of $Z$ becomes:
$$A =  -\frac{\partial \ln Z}{\partial \beta} - T\,\frac{\partial (kT\ln Z)}{\partial T} =
kT^2\frac{\partial \ln Z}{\partial T} - kT\,\left(\frac{T \partial \ln Z}{\partial T} +\ln Z \right),$$
\begin{equation}A = - kT\,\ln Z. \end{equation}
It is clear now that the partition function in spectroscopy is just one implementation of the canonical PF. Furthermore, the internal energy that is commonly used quantity in the radiative MHD problems, is clearly the total energy $E$ derived above. The same holds for the entropy and the Helmholtz free energy. These definitions can be found in virtually any textbook on thermodynamics and statistical physics, on Wikipedia as well. Still, my favorite introduction to this topic remains Chandrasekhar's Introduction to the Study of Stellar Structure (1958). Papers of Mihalas et al on the so-called MHD equation of state are on the higher level, but nicely written (Mihalas et al, 1988ApJ...331..794H, 1988ApJ...331..815M, 1988ApJ...332..261D, 1990ApJ...350..300M).


IDL routines

On the other hand, in practice, it has a finite values as the levels above certain threshold are affected by the collisions and do not contribute to the sum any more.

Numerical computations that take into account detailed atomic physics and use various approximations are performed by various authors and groups. The coefficients of the best polynomial fits are usually tabulated for a range of temperature values.

There is a number of papers that provide the partition function values. However, they vary in terms of the periodic-system coverage, of the ionization stages that are included and in the form of the interpolation polynomial. I collected data from five different sources and sorted them out in IDL so that it is easy to load and interpolate partition function for a given element, ionization stage and temperature. If you're interested in the details of the solution behind any of these partition functions, please check the original papers.

The sources are:

Gray, 20052005oasp.book.....Gpf_gray.dat
Irwin, 19811981ApJS...45..621Ipf_irwin.sav
Sauval and Tatum, 19841984ApjS...56..193Spf_sauval_and_tatum.dat
Wittmann, 19741974SoPh...35...11Wpf_wittmann.pro
Cowleyafter 2001A&A...374..265Wpf_cowley.dat

If you use any of the data in a publication, please acknowledge it by citing the corresponding paper.

There is an IDL routine for each of the sources. There is also a file with the coefficients for each of them (except for Wittmann's PF - these coefficients are hardcoded in the corresponding routine). The IDL routines simply read the tables and do the interpolation. When the data is missing data, PF is set to -1.

Instead of calling particular IDL functions, it is recommended to call the pf.pro function that takes the source of PF data as an input and calls the corresponding function.

The tarball is here, pf.tar. The individual functions are:

pf.pro
partition_gray.pro
partition_irwin.pro
partition_sauval_and_tatum.pro
partition_wittmann.pro
partition_cowley.pro


Installation

Before running the code for the first time, copy all the functions into a directory within your IDL path. Open pf_path.pro and set the path to that exact directory. That's all.


Example

Let's compare PFs from different sources for the neutral manganese:
IDL> z = 25
IDL> t = 2000. + DINDGEN(6)*2000.
IDL> pfw = PF(t, z, data = 'w')
IDL> pfi = PF(t, z, data = 'i')
IDL> pfg = PF(t, z, data = 'g')
IDL> pfs = PF(t, z, data = 's')
IDL> pfc = PF(t, z, data = 'c')


No comments:

Post a Comment