text:
# Mean-field solution of the parity-conserving kinetic phase transition in one dimension

## I Introduction

Phase transitions occurring away from thermodynamical equilibrium constitute one of the most challenging topics in statistical physics. They appear in a host of physical systems as well as in many models in biology, chemistry, sociology, etc. Given the lack of a general theory of non-equilibrium systems, we are still living a taxonomic era in this field: it would be highly desirable to reach a complete classification of the known phase transitions into universality classes, as a preliminary step to their full categorization, and identification of relevant features.

After more than twenty years of study, it has become clear that the degree of universality is much narrower away from equilibrium than it is in equilibrium transitions. Certainly, general properties such as symmetries, conservation laws, dimensionalities, etc., play a key role, as they do in equilibrium. But some other ingredients, such as microscopic dynamical details, hard core interactions, and the type of updating may, in some cases, influence the emerging critical behavior of non-equilibrium systems, making the task of theoreticians both stimulating and difficult.

One of the most robust and best studied non-equilibrium classes of transitions is directed percolation (DP). It includes an amazing variety of models and systems exhibiting a transition from an active to an absorbing phase Kin ; DK ; CP ; Reviews ; Marro ; IAS , and is well characterized at a field-theoretical level by Reggeon field theory (RFT) conjecture ; RFT .

DP is so robust, that identifying and classifying the nature of perturbations able to drive absorbing phase transitions away from this class has become a challenging task. Probably the best known instance of this is the $`DP2`$ class, where the presence of two perfectly symmetric absorbing states ($`Z_2`$ symmetry) is the main feature responsible for non-DP scaling equiv . This class includes a cellular automaton introduced by Grassberger et al., peter , interacting monomer-dimer models Park , non-equilibrium Ising models Menyhard ; Odor , monomer-monomer surface reaction models Brown , $`Z_2`$-symmetric generalizations of DP Hinrichsen  and of the contact process Inui , branching annihilating random walks with conserved parity TT ; Santos ; Bramson , and it has also been related to generalized versions of the Voter model Voter1 ; Voter2 . Extensive numerical simulations led to conjecture rational exponent values for the DP2 class Jensen , though this was later disproved by more exhaustive simulations Reviews .

Although the existence and robustness of the DP2 class are well established from a numerical view-point, a solid theoretical understanding is still missing. Bold attempts to write down and renormalize a field theory suitable for the DP2 class have been performed, but the results are not as satisfactory as they are for RFT tauber ; Korea . In particular, the renormalization is based on a clever, but somehow uncontrolled expansion around two-different critical dimensions. As we shall illustrate below, straightforward mean-field approaches (cluster approximations) fail to reproduce a phase transition if too small clusters are considered, and one has to resort to large-clusters, which make the calculation complicated and not very accurate, i.e. the convergence of the series is very slow (although results can be improved if combined with coherent-anomaly methods Odor ).

It is the purpose of this paper to shed some light on these issues, by examining models of Branching Annihilating Walks with two offspring (2-BAW) in one dimension: the 2-BAW with finite reaction rate exhibits a phase transition in the DP2 class 2-BAW . Our approach is based on a mean-field modification of the method of parity intervals, originally introduced for the study of annihilation reactions, $`A+A0`$ book ; parity\_ints ; masser ; lind . A similar approach has been previously employed in conjunction with the method of empty intervals book ; ddb ; Henkel  for the analysis of other intractable models that do not conserve parity approx . See parity\_ints ; masser  for a more detailed introduction to this method and its applications to different models.

The paper is organized as follows. The 2-BAW model is described in Section II, where we also perform simple (up to two sites) cluster approximations (which fail to capture the transition). In Section III we present the exact solution of the 2-BAW model in the limit of infinite reaction rate, using the method of parity-intervals. Although there is no phase transition in this limit, the analysis serves as a basis for the mean-field like approximation presented in Section IV, for the relevant case of finite reaction rates. There we derive the approximate steady-state solution for generic parameter values, and compare it with the outcome of Monte Carlo computer simulations. We conclude with a critical discussion of our results and further developments, in Section V.

## II The model and cluster approximations

The model is defined as follows. Each site of a one-dimensional lattice is either empty or singly occupied. The lattice is updated asynchronously: an occupied site is randomly chosen, and it is tried for diffusion, at rate $`\mathrm{\Gamma }`$ (probability $`\mathrm{\Gamma }/(\mathrm{\Gamma }+\mathrm{\Omega })`$), or branching, at rate $`\mathrm{\Omega }`$ (probability $`\mathrm{\Omega }/(\mathrm{\Gamma }+\mathrm{\Omega })`$); time is increased by $`1/N`$, where $`N`$ is the number of occupied (active) sites. In a diffusion attempt the particle is moved to one of its two nearest neighbors (target sites), with equal probabilities. If the target site is occupied, annihilation results with probability $`r`$: the move is effected and both particles, the diffuser and the target, are removed from the lattice. The move is rejected with probability $`1r`$, and the lattice state remains unchanged. In a branching attempt the particle gives birth to two new particles onto the nearest neighbor sites (target sites). If either, or both of the target sites is occupied, annihilation takes place with probability $`r`$: branching is effected, and the occupied target site(s) become empty. Once again, the move is rejected with probability $`1r`$, and the lattice state remains unchanged.

In the diffusion-controlled limit of infinite reaction rate ($`r=1`$), the 2-BAW model is known to possess only a steady absorbing phase Sudbury . However, for finite rates ($`r<1`$), a transition belonging to the DP2 class, from an absorbing state, at $`r<r_\mathrm{c}`$, to an active phase, at $`r>r_\mathrm{c}`$, is found in numerical simulations 2-BAW . We now demonstrate that simple mean-field approximations fail to capture this transition.

The simplest conceivable mean-field theory is that of the one-site approximation, obtained by neglecting all correlations between the states of different sites. For example, if the probability of one site being occupied is equal to the concentration $`c`$, then the probability of finding one occupied site followed immediately by an empty site is $`c(1c)`$. There are only three events that lead to a change in particle concentration: (i) annihilation of two particles by diffusion, $``$; (ii) creation of two particles via birth, $``$; and (iii) annihilation of two particles due to birth onto previously occupied sites, $``$. Note that the birth process $``$ (and its mirror-symmetric image) does not alter the concentration. However, processes (i) and (iii) occur with restricted probability $`r`$. Thus

$$\frac{d}{dt}c=2rc^2+2c(1c)^22rc^3,$$
(1)
where we took $`\mathrm{\Gamma }=\mathrm{\Omega }=1`$. The first, second and third terms on the r.h.s. correspond to (i), (ii) and (iii), respectively.

For the steady state, we set $`dc/dt=0`$ and choose the stable root (Fig. 1), $`c_\mathrm{s}=(2+r\sqrt{8r+r^2})/(22r)`$. We conclude that according to the one-site approximation the system always evolves to an active phase, with $`1/3<c_\mathrm{s}<1`$ (Fig.2).

The two-site approximation provides the next level of complexity. Let $`p`$ be the conditional probability for a site to be occupied, given that the adjacent site (to its left) is occupied. Let $`q`$ be the conditional probability for a site to be occupied, given that the site to its left is empty. Thus, the probability of the events $``$, $``$, $``$, $``$ is $`cp`$, $`c(1p)`$, $`(1c)q`$, $`(1c)(1q)`$, respectively. Note that since $``$ and $``$ are equally likely, we have,

$$c(1p)=(1c)q.$$
(2)
(The same relation may be derived from the fact that $`\mathrm{Pr}()=\mathrm{Pr}()+\mathrm{Pr}()`$.) Eq. (1) is now rewritten as

$$\frac{d}{dt}c=2rcp+2(1c)q(1p)2rcp^2,$$
(3)
where the terms on the r.h.s. correspond to the same events as before. Setting $`dc/dt=0`$ and using the relation (2), we obtain, for the steady state,

$$0=2rcp+2c(1p)^22rcp^2.$$

We conclude that either $`c_\mathrm{s}=0`$ or $`p_\mathrm{s}=(2+r\sqrt{8r+r^2})/(22r)`$. In the latter case, an additional evolution equation, for the event $``$, provides the missing value of $`c_\mathrm{s}`$: $`c_\mathrm{s}>0`$ for all $`r`$. The result agrees closely with that of the one-site approximation, with deviations not larger than 20 %. Thus, also the two-site approximation fails to predict the transition. Considering larger clusters, a phase transition can be generated Odor ; Attila , but a large set of coupled equations (which enlarges with cluster-size) has to be solved, the accuracy is not good, and the results are expected to be valid only up to the cluster-size. In the next section we present an alternative method intended to overcome these difficulties.

## III Parity intervals and exact solution for $`r=1`$

We now turn to a different approach, that of the method of parity intervals book ; parity\_ints ; masser ; lind . We first present the exactly soluble case of $`r=1`$, for which there is no transition. The exact approach followed here serves as a basis for a mean-field like approximation, for the more interesting case of $`r<1`$ — an approximation that does capture the transition and reproduces kinetic details surprisingly well.

Let $`G_n(t)`$ be the probability that (in an homogeneous system) an arbitrary segment of $`n`$ consecutive sites contains an even number of particles at time $`t`$. A site can be either empty or occupied by a single particle, so the probability that a site is occupied, i.e., the particle density, is

$$\rho (t)=1G_1(t).$$
(4)
Since the dynamic rules of the two-offspring BAW conserve parity, the only way that $`G_n`$ might change is when: (i) particles at the edge of the segment hop outside or branch, or (ii) particles just outside of the segment hop inside or branch. Let   and   represent segments of even and odd number of particles, respectively. Then, the changes of $`G_n`$ can be described schematically by the events:

$`{\displaystyle \frac{d}{dt}}G_n`$ $`=`$  n

+ n1

 n

 n1

superscript 𝑛

superscript 𝑛1

superscript 𝑛

superscript 𝑛1

\displaystyle{\mathop{\rule{30.00005pt}{4.30554pt}}^{n}{\shortstack{$\!\!\!\leftarrow$
\\
$\bullet$}}}+{\mathop{\framebox{\rule{30.00005pt}{0.0pt}\rule{0.0pt}{4.30554pt}}}^{n-1}{\shortstack{$\rightarrow\!\!$ \\
$\bullet$}}}-{\mathop{\framebox{\rule{30.00005pt}{0.0pt}\rule{0.0pt}{4.30554pt}}}^{n}{\shortstack{$\!\!\!\leftarrow$ \\
$\bullet$}}}-{\mathop{\rule{30.00005pt}{4.30554pt}}^{n-1}{\shortstack{$\rightarrow\!\!$ \\
$\bullet$}}} (5)
$`+`$  n

|

+ n1

|

 n

|

 n1

|

.superscript 𝑛

|

superscript 𝑛1

|

superscript 𝑛

|

superscript 𝑛1

|

\displaystyle{\mathop{\rule{30.00005pt}{4.30554pt}}^{n}{\shortstack{$\swarrow|\searrow$ \\
$\bullet$}}}+{\mathop{\framebox{\rule{30.00005pt}{0.0pt}\rule{0.0pt}{4.30554pt}}}^{n-1}{\shortstack{$\swarrow|\searrow$ \\
$\bullet$}}}-{\mathop{\framebox{\rule{30.00005pt}{0.0pt}\rule{0.0pt}{4.30554pt}}}^{n}{\shortstack{$\swarrow|\searrow$ \\
$\bullet$}}}-{\mathop{\rule{30.00005pt}{4.30554pt}}^{n-1}{\shortstack{$\swarrow|\searrow$ \\
$\bullet$}}}\;.
Arrows in the first four terms indicate the hopping of a particle to the left or right. In the last four terms, the arrows indicate branching. For example, the first term in (5) represents the event that a particle outside of an $`n`$-segment of odd parity jumps in, thus creating an $`n`$-segment of even parity.

For the case of immediate reactions, $`r=1`$, the processes indicated in (5) occur regardless of the state of target sites. It is easy to show (see book ; parity\_ints ; masser ; lind ) that:

$`F_n\mathrm{Pr}(\stackrel{n}{\text{ }\text{}})={\displaystyle \frac{1}{2}}[(1G_1)+(G_nG_{n+1})],`$ (6a)
$`H_n\mathrm{Pr}(\stackrel{n}{\text{ }})={\displaystyle \frac{1}{2}}[(1G_1)(G_nG_{n+1})].`$ (6b)
Using these relations, Eq. (5) becomes

$$\frac{d}{dt}G_n(t)=2(\mathrm{\Gamma }+\mathrm{\Omega })(G_{n1}2G_n+G_{n+1}).$$
(7)
$`\mathrm{\Gamma }`$ and $`\mathrm{\Omega }`$ are the rates of hopping and branching, respectively, and the factor of 2 accounts for the events in (5) taking place also at the left edge of the interval. The case of $`n=1`$ requires a special equation, since $`G_0`$ is undefined. Taking into account all the ways $`G_1`$ might change, we find

$$\frac{}{t}G_1(t)=2\mathrm{\Gamma }(12G_1+G_2)+2\mathrm{\Omega }(G_2G_1).$$
(8)
Thus, Eq. (7) may be understood to be valid also for $`n=1`$, provided that one uses the boundary condition

$$G_0=\frac{\mathrm{\Gamma }+\mathrm{\Omega }G_1}{\mathrm{\Gamma }+\mathrm{\Omega }}.$$
(9)
Additionally, since the $`G_n`$ are probabilities, we have

$$0G_n(t)1.$$
(10)
Eq. (7), with the boundary conditions (9), (10) may be analyzed exactly for a variety of initial conditions parity\_ints . Here we merely observe that the only steady state solution supported by these equations is $`G_n=1`$, $`n=1,2,\mathrm{}`$, corresponding to the absorbing state $`\rho _\mathrm{s}=0`$. Indeed, the 2-BAW model with $`r=1`$ lacks an active phase and exhibits no transition.

## IV Parity intervals approximation for the general case

Consider now the case of finite reaction probability, $`r<1`$. A transition about some critical probability value $`r_\mathrm{c}`$ is known to take place, from the absorbing state into an active (non-empty) phase 2-BAW , and our goal is to capture this transition, if only in an approximate fashion. We shall assume that $`\mathrm{\Gamma }=\mathrm{\Omega }=1`$, without loss of generality remark .

Consider the first process on the r.h.s. of (5). For $`r<1`$, there is a difference in the reaction rate depending on whether the target site is empty or occupied:

n

= n1

+r n1

.superscript 𝑛

superscript 𝑛1

𝑟superscript 𝑛1

\mathop{\rule{30.00005pt}{4.30554pt}}^{n}{\shortstack{$\!\!\!\leftarrow$\\
 $\bullet$}}=\mathop{\rule{30.00005pt}{4.30554pt}}^{n-1}\circ{\shortstack{$\!\!\!\leftarrow$\\
 $\bullet$}}+r\,\mathop{\framebox{\rule{30.00005pt}{0.0pt}\rule{0.0pt}{4.30554pt}}}^{n-1}\bullet{\shortstack{$\!\!\!\leftarrow$\\
 $\bullet$}}\;.
Unfortunately, events such as $`\text{ }`$ and $`\text{ }\text{}`$ cannot be expressed in closed form in terms of the $`G_n`$. We therefore rewrite the process in a way that the terms associated with $`r<1`$, and which cannot be expressed in closed form, appear as a perturbation, proportional to $`s1r`$:

n

superscript 𝑛

\displaystyle\mathop{\rule{30.00005pt}{4.30554pt}}^{n}{\shortstack{$\!\!\!\leftarrow$\\
 $\bullet$}} = n1

+r n1

absentsuperscript 𝑛1

𝑟superscript 𝑛1

\displaystyle=\mathop{\rule{30.00005pt}{4.30554pt}}^{n-1}\circ{\shortstack{$\!\!\!\leftarrow$\\
 $\bullet$}}+r\,\mathop{\framebox{\rule{30.00005pt}{0.0pt}\rule{0.0pt}{4.30554pt}}}^{n-1}\bullet{\shortstack{$\!\!\!\leftarrow$\\
 $\bullet$}}
$`=\stackrel{n1}{\text{ }}+\stackrel{n1}{\text{ }\text{}}(1r)\stackrel{n1}{\text{ }\text{}}`$
$`=\stackrel{n}{\text{ }}s\stackrel{n1}{\text{ }\text{}}.`$
Likewise, the remainder of the diffusion events in (5) may be rewritten as
n1

= n1s n1,\displaystyle\mathop{\framebox{\rule{30.00005pt}{0.0pt}\rule{0.0pt}{4.30554pt}}}^{n-1}\shortstack{$\rightarrow\!\!$\\
$\bullet$}=\mathop{\framebox{\rule{30.00005pt}{0.0pt}\rule{0.0pt}{4.30554pt}}}^{n-1}\bullet-s\,\mathop{\framebox{\rule{30.00005pt}{0.0pt}\rule{0.0pt}{4.30554pt}}}^{n-1}\bullet\bullet\;, (11b)
n

= ns n1,\displaystyle\mathop{\framebox{\rule{30.00005pt}{0.0pt}\rule{0.0pt}{4.30554pt}}}^{n}\shortstack{$\!\!\!\leftarrow$\\
$\bullet$}=\mathop{\framebox{\rule{30.00005pt}{0.0pt}\rule{0.0pt}{4.30554pt}}}^{n}\bullet-s\,\mathop{\rule{30.00005pt}{4.30554pt}}^{n-1}\bullet\bullet\;, (11c)
n1

= n1s n1.\displaystyle\mathop{\rule{30.00005pt}{4.30554pt}}^{n-1}\shortstack{$\rightarrow\!\!$\\
$\bullet$}=\mathop{\rule{30.00005pt}{4.30554pt}}^{n-1}\bullet-s\,\mathop{\rule{30.00005pt}{4.30554pt}}^{n-1}\bullet\!\bullet\;. (11d)
The branching events too require close inspection of the target sites. For example, we rewrite the first branching event of (5) in a perturbative fashion:

n

|

superscript 𝑛

|

\displaystyle\mathop{\rule{30.00005pt}{4.30554pt}}^{n}{\shortstack{$\swarrow|\searrow$ \\
 $\bullet$}} = n1

|

+r n1

|

+r n1

|

+r n1

|

\displaystyle=\mathop{\rule{30.00005pt}{4.30554pt}}^{n-1}\circ{\shortstack{$\swarrow|\searrow$ \\
$\bullet$}}\circ+r\,\mathop{\framebox{\rule{30.00005pt}{0.0pt}\rule{0.0pt}{4.30554pt}}}^{n-1}\bullet{\shortstack{$\swarrow|\searrow$ \\
$\bullet$}}\circ+r\,\mathop{\rule{30.00005pt}{4.30554pt}}^{n-1}\circ{\shortstack{$\swarrow|\searrow$ \\
$\bullet$}}\bullet+r\,\mathop{\framebox{\rule{30.00005pt}{0.0pt}\rule{0.0pt}{4.30554pt}}}^{n-1}\bullet{\shortstack{$\swarrow|\searrow$ \\
$\bullet$}}\bullet
$`=(1r)\stackrel{n1}{\text{ }}+r\stackrel{n1}{\text{ }}+r\stackrel{n1}{\text{ }\text{}}+r\stackrel{n1}{\text{ }}+r\stackrel{n1}{\text{ }\text{}}`$
$`=s\stackrel{n1}{\text{ }}+(1s)\stackrel{n}{\text{ }}.`$
The remainder of the branching terms are similarly expressed as
n1

|

=s n2+(1s) n1,\displaystyle\mathop{\framebox{\rule{30.00005pt}{0.0pt}\rule{0.0pt}{4.30554pt}}}^{n-1}{\shortstack{$\swarrow|\searrow$ \\
 $\bullet$}}=s\,\mathop{\framebox{\rule{30.00005pt}{0.0pt}\rule{0.0pt}{4.30554pt}}}^{n-2}\circ\!\bullet\!\circ+(1-s)\,\mathop{\framebox{\rule{30.00005pt}{0.0pt}\rule{0.0pt}{4.30554pt}}}^{n-1}\bullet\;, (12b)
n

|

=s n1+(1s) n,\displaystyle\mathop{\framebox{\rule{30.00005pt}{0.0pt}\rule{0.0pt}{4.30554pt}}}^{n}{\shortstack{$\swarrow|\searrow$ \\
 $\bullet$}}=s\,\mathop{\framebox{\rule{30.00005pt}{0.0pt}\rule{0.0pt}{4.30554pt}}}^{n-1}\circ\!\bullet\!\circ+(1-s)\,\mathop{\framebox{\rule{30.00005pt}{0.0pt}\rule{0.0pt}{4.30554pt}}}^{n}\bullet\;, (12c)
n1

|

=s n2+(1s) n1.\displaystyle\mathop{\rule{30.00005pt}{4.30554pt}}^{n-1}{\shortstack{$\swarrow|\searrow$ \\
 $\bullet$}}=s\,\mathop{\rule{30.00005pt}{4.30554pt}}^{n-2}\circ\!\bullet\!\circ+(1-s)\,\mathop{\rule{30.00005pt}{4.30554pt}}^{n-1}\bullet\;. (12d)
So far everything is exact. In order to proceed, we approximate the terms proportional to $`s`$ in the simplest possible way, by neglecting correlations. Thus, for example, we write

$$\mathrm{Pr}(\stackrel{n1}{\text{ }\text{}})\mathrm{Pr}(\stackrel{n1}{\text{ }\text{}})\mathrm{Pr}()=F_{n1}(1G_1),$$
(13)
for the problematic term in (11a), and similar expressions for the ones in (11b)–(11d). For the problematic terms proportional to $`s`$ in (12), we introduce the notation

$$\mathrm{Pr}(\stackrel{n}{\text{ }\text{}})=f_n,\mathrm{Pr}(\stackrel{n}{\text{ }})=h_n,$$

and

$$\mathrm{Pr}()=x,\mathrm{Pr}()=\mathrm{Pr}()=y,\mathrm{Pr}()=z,$$

and approximate in the same spirit as above,

$$\mathrm{Pr}(\stackrel{n1}{\text{ }})\mathrm{Pr}(\stackrel{n1}{\text{ }})\mathrm{Pr}()=h_{n1}y,$$
(14)
and similarly for the other terms. Since $`y+z=G_1`$, $`x+z=G_2`$, and $`x+2y+z=1`$, it follows that $`x=\frac{1}{2}(12G_1+G_2)`$, $`y=\frac{1}{2}(1G_2)`$, and $`z=\frac{1}{2}(1+2G_1+G_2)`$ are expressible in terms of the $`G_n`$ in closed form. Since $`f_n+F_n=G_n`$ and $`h_n+H_n=1G_1`$, it follows that $`f_n=\frac{1}{2}(1+G_1+G_n+G_{n+1})`$ and $`h_n=\frac{1}{2}(1+G_1G_nG_{n+1})`$ are too given in terms of the $`G_n`$ in closed form.

We are now ready to write down a closed evolution equation for $`G_n`$. Starting from Eq. (5), using the representations of Eqs. (11) and (12) with their respective approximations, collecting terms and rearranging, we obtain

$$\frac{d}{dt}G_n=\frac{1}{2}s(1G_2)G_{n2}+(23s+2sG_1)G_{n1}+\frac{1}{2}(8+7s4sG_1+sG_2)G_n+(2s)G_{n+1},$$
(15)
valid for $`n3`$. As expected, this reduces to (7), with $`\mathrm{\Gamma }=\mathrm{\Omega }=1`$, in the limit $`s0`$. The equation for $`n=2`$ is exactly the same as (15), provided that one adopts the boundary condition

$$G_0=1.$$
(16)
$`G_1`$ requires a separate evolution equation that we find by considering all the events that contribute to $`dG_1/dt`$, and write down, in the spirit of Eq. (15), as

$`{\displaystyle \frac{d}{dt}}G_1`$ =2(

+

+

|

|

)\displaystyle=2(\bullet\shortstack{$\!\!\!\leftarrow$\\
 $\bullet$}+\shortstack{$\rightarrow\!\!$\\
 $\bullet$}-\circ\shortstack{$\!\!\!\leftarrow$\\
 $\bullet$}+\bullet\shortstack{$\swarrow|\searrow$ \\
 $\bullet$}-\circ\shortstack{$\swarrow|\searrow$ \\
 $\bullet$})
$`2\{(1s)x+[(1G_1)sx]y+(1s)x[(1s)y+syG_1]\}`$
$`=26G_1+4G_2+s(2+5G_14G_2+G_1G_2).`$
To compute the steady state, we set $`dG_n/dt=0`$. Eq. (15) yields a recursion relation for the $`G_n`$, with constant coefficients (that depend partly on $`G_1`$ and $`G_2`$). Applying the ansatz $`G_n=\lambda ^n`$ we find a cubic equation for $`\lambda `$, with roots

$$\lambda _0=1,\lambda _\pm =\frac{45s+4sG_1sG_2\pm \sqrt{(45s+4sG_1sG_2)^2+8s(2s)(1G_2)}}{84s}.$$
(18)
Thus $`G_n`$ has the general solution

$$G_n=A\lambda _+^n+B\lambda _{}^n+C,$$

where $`A`$, $`B`$, $`C`$ are constants to be determined from boundary conditions. Since $`0<s,G_1,G_2<1`$ (in the active phase) it follows that $`|\lambda _\pm |<1`$, and $`C=lim_n\mathrm{}G_nG_{\mathrm{}}`$. Suppose that the initial distribution of particles is random, at density $`\rho _0`$, then $`G_n(t=0)=\frac{1}{2}+\frac{1}{2}(12\rho _0)^n`$, and $`G_{\mathrm{}}(0)=1/2`$ parity\_ints . Because our model conserves parity, it follows that $`C=1/2`$. Furthermore, the boundary condition (16) implies $`B=\frac{1}{2}A`$. The remaining coefficient, $`A`$, could be found from the relations

$`G_1=A\lambda _++({\displaystyle \frac{1}{2}}A)\lambda _{}+1/2,`$ (19a)
$`G_2=A\lambda _+^2+({\displaystyle \frac{1}{2}}A)\lambda _{}^2+1/2,`$ (19b)
$`G_2={\displaystyle \frac{2s2+(65s)G_1}{44s+sG_1}}.`$ (19c)
The first two relations are required for self-consistence, since $`\lambda _\pm `$ depend on $`G_1`$ and $`G_2`$ (as well as on $`s`$). Eq. (19c) is derived from (IV), with $`dG_1/dt=0`$. Because Eqs. (19) are unwielding, in practice we set a numerical value for $`G_1`$ and search for $`(A,s,G_2)`$ that satisfy (19), using Mathematica. We thus find a kinetic phase transition which we next compare to simulations.

The critical point $`r_\mathrm{c}=1s_\mathrm{c}`$, and the behavior of the order parameter in its vicinity, can be obtained analytically, by perturbing (19) about the values $`A=0`$, $`G_1=1`$. We thus find $`r_\mathrm{c}=\frac{1}{6}\sqrt{33}\frac{1}{2}=0.457427\mathrm{}`$ and $`\rho _\mathrm{s}(r_\mathrm{c}r)^\beta `$, with $`\beta =1`$. The value of $`\beta =1`$ is in accordance with the expected from a mean-field like solution.

We can also obtain analytically the behavior of the order parameter in the neighborhood of $`r=0`$. The answer, $`c_\mathrm{s}1\sqrt{2r}`$, is the same as the one obtained from the one-site and two-site cluster approximations, and seems to fit the numerical data from simulations exactly. Clearly, the role of fluctuations diminishes as $`r`$ decreases away from the critical point, and it is conceivable that their impact is negligible in the limit of zero reaction.

### Comparison to Simulations

In order to test the limit of validity of the above calculations, we have performed a Monte Carlo simulation of the 2-BAW model. We use a system size $`L=10^4`$, and implement the previously discussed rules. The initial configuration consists of a randomly half-occupied lattice; we have verified that the obtained long-time asymptotic results are insensitive to variations of the initial condition. Plotting the average density as a function of time in a double logarithmic plot, we identify the critical point with the value of $`r`$ for which a separatrix, between curves that tend asymptotically to a constant (active phase) and curves that bend downwards (absorbing phase) converging asymptotically to a $`t^{1/2}`$ decay exp  is obtained. We thus find $`r_c=0.495(10)`$. From the slope of the plot the density decay critical exponent $`\theta =0.28(1)`$ (extending for about four decades), is determined. It is in excellent agreement with previous measurements in the DP2 class Reviews . The stationary density values are represented in Fig. 2.

A point worth emphasizing is that finite-size corrections differ from those usually encountered in systems with absorbing states: Rather than finding a larger stationary order parameter for smaller systems, as is common, here one observes a faster decay to the absorbing phase. Owing to this the exponent ratio $`\beta /\nu `$ cannot be determined from finite-size scaling analysis.

We have also measured the order-parameter critical exponent, by plotting the average stationary density as a function of the distance to the critical point (Fig. 3). Our finding of $`\beta =0.92(5)`$ is again in good agreement with the commonly accepted value Reviews . Other exponents can certainly be measured, but with $`\theta `$ and $`\beta `$ we can already guarantee that the model is indeed in the DP2 class.

A more stringent test of the mean-field parity interval approximation is obtained by comparing its predictions for $`G_n`$ to the numerical results of Monte Carlo simulations, for different values of the control parameter $`r`$ (Fig. 4). In all cases, we expect $`G_n1/2`$ for large $`n`$, since $`G_{\mathrm{}}=1/2`$ and parity is conserved. The relevant question is how well the approximation captures the $`G_n`$ for $`n`$ small. Because the typical distance between particles grows as one approaches criticality, we plot the results against $`\nu =\rho n`$ rather than $`n`$, in order to compare them better.

For small values of $`r`$ (as shown in Fig. 4), deep into the active phase, spatial correlations play a minimal role and the theoretical prediction matches simulations remarkably well, for all $`\nu `$. As $`r`$ is raised toward the critical point, correlations play a larger role and the fit worsens. For example for $`r=0.38`$ ($`\rho =0.2`$) the theoretical curve fits experiments well only up to $`\nu 2`$ (and, of course, also at $`\nu 1`$). Closer still to the critical point, at $`r=0.45`$ ($`\rho =0.1`$), the theoretical prediction is good only up to $`\nu 1`$. The large deviations observed for $`\nu 1`$ in this last case indicate that the state of the system near the DP2 transition, despite being dilute, is highly ordered and correlated.

The approximate curves, give converge to the exact value for $`G_n1/2`$ for asymptotically large values of $`\nu `$.

## V Discussion and future perspectives

We have presented an approximation, based on the method of parity intervals, for the analysis of the DP2 transition observed in the 2-BAW model in one dimension with finite reaction rates. The key ingredient of our approach is that it respects the parity conservation implicit in the dynamical rules of the 2-BAW model: parity conservation is here responsible for the $`Z_2`$ symmetry that underlies transitions in the DP2 class, and provides a surprisingly good description of the transition at a mean-field level.

An interesting finding is the fact that spatial correlations play no role in the limit of zero reaction rate, $`r0`$. In this “saturation” limit, the dependence of the order parameter upon the critical field, $`\rho 1\sqrt{2r}`$, is captured even by the simplest one-site cluster approximation.

Extensive numerical simulations of the DP2 transition suggest a value of $`\beta =0.92(2)`$ for the order parameter critical exponent Reviews . On the other hand, an exact value $`\beta =1`$ was conjectured by Jensen and later supported by numerical simulations combined with Padé approximants Inui . Regrettably, we have little to add to this controversy. Our mean-field results away from criticality ($`r_\mathrm{c}r0.1`$) are consistent with $`\beta =0.92`$, while $`\beta 1`$ as we approach the transition (Fig. 3). However, the large spatial correlations observed near criticality are not faithfully modeled by our mean-field like approach (Fig. 4), and the values we obtain for $`r_\mathrm{c}r0.1`$ are unreliable. It would be desirable to improve our approximation so as to better model regions closer to $`r_\mathrm{c}`$. Our attempts in this direction have been futile, so far. Arguably, the simplest improvement would be to replace the approximation of Eq. (13) by the more accurate:

$$\mathrm{Pr}(\stackrel{n1}{\text{ }\text{}})\frac{\mathrm{Pr}(\stackrel{n1}{\text{ }\text{}})\mathrm{Pr}()}{\mathrm{Pr}()},$$

and likewise elsewhere. However, on employing this scheme, instead of an improved result we find that the transition disappears. A similar phenomenon is known to occur also in mean-field cluster approximations, where sometimes increasing the cluster size yields lower quality predictions.

Other interesting open prospects include using the parity interval approximation for the analysis of dynamical aspects of the DP2 transition. We have here systematically assumed that all time derivatives are zero, thereby accessing the steady state alone. Analyzing the very same equations with the full time dependence built in should yield a prediction for the order parameter time decay. Additionally, the equations could be modified to describe a spatially inhomogeneous system. In that case, one could study the dynamics of spreading.

###### Acknowledgements.

We acknowledge financial support from the Spanish MCyT (FEDER) under project BFM2001-2841. DbA thanks the NSF for partial support, under contract no. PHY-0140094, as well as the warm hospitality of the Instituto de Física Teórica y Computacional Carlos I, University of Granada, during a crucial phase of the project. We also gratefully thank G. Odor and A. Szolnoki, for very useful comments on the cluster mean field calculations that have helped to improve the paper.


format: text

text:
# Dielectronic recombination data for dynamic finite-density plasmas

## 1 Introduction

Dielectronic recombination (DR) is the dominant electron–ion recombination process in many astrophysical and laboratory plasmas. It plays an important role in determining both the level populations and the ionization balance of both high- and low-temperature non-LTE plasmas over a wide range of electron densities. Dielectronic recombination can be viewed as a two-step process. Firstly, a free electron excites an electron in an ion $`X^{+z}`$, say, and in the process transfers sufficient of its energy that it is captured into an autoionizing state of the ion $`X^{+z1}`$. The process is reversible since the total energy of the system remains conserved. However, if an electron (either the captured electron or an electron in the parent core) then makes a spontaneous radiative transition that leaves the ion in a non-autoionizing (bound) state then the recombination can be viewed as complete, at least in the low-density (coronal) limit. Dielectronic recombination can take place via many intermediate autoionizing states – indeed, entire Rydberg series. It was the significance of the effective statistical weight of so many available states that led Burgess (1964) to recognize the importance of dielectronic recombination in the first instance. A consequence of this is that a huge population structure model is required in principle for further progress. In astrophysics and fusion, the problem has usually been made manageable by simply summing over all final states so as to produce a total dielectronic recombination rate coefficient. In combination with radiative recombination, this is the effective recombination coefficient from the point of view of an ionization state only – the so-called ‘coronal picture’. In the latter, excited-state populations are depleted exclusively by spontaneous radiative transitions and are small compared to those of ground states – with which they are in quasi-static equilibrium. Collisional processes are negligible, except with ground state targets. In turn, the ionization state is determined by balancing the effective recombination rate of (the ground state of) $`X^{+z}`$ against the collisional ionization rate of (the ground state of) $`X^{+z1}`$.

For dynamic finite-density plasmas, there are therefore two critical limitations to the coronal approximation which must be addressed. Both require a re-appraisal of the strategy for computing dielectronic recombination, which is the objective of this paper. Firstly, from the atomic point of view, a dynamic plasma is one for which the timescale of change in plasma parameters (especially electron temperature, $`T_\mathrm{e}`$, and electron density, $`N_\mathrm{e}`$) is comparable with the lifetime of metastable populations of its constituent ions. This situation is a concern, for example, for impurities near plasma contacted surfaces in fusion (Summers et al. 2002) and for active solar events (Lanza et al. 2001). The population of an ionization stage may no longer be assumed to be concentrated in the ground state of each ion. A significant population is found in the metastables, and they are not in quasi-static equilibrium with the ground state. This means that such metastables may be the starting point for recombination events, and so the time-evolution of their populations must be tracked in the same manner as for the ground states. We call this the generalized collisional–radiative (GCR) picture (Summers & Hooper 1983).

The second issue is that of finite-density effects. The coronal (zero-density limit) picture assumes that, following the two-step dielectronic recombination process, the resultant (non-autoionizing) excited-state (electron) radiatively cascades back down to the ground state without collisional disruption. At finite electron densities, this radiative cascade can be interrupted by further electron collisions which redistribute the population – in particular, ionization (possibly stepwise) out of excited states, which reduces the effective dielectronic recombination rate. Collisional–radiative modelling removes the limitations of the coronal model, but at the cost of much more elaborate excited population calculations. These in turn require much more detailed dielectronic recombination data, and in an easily accessible form. The onset of these density effects on the ‘post-DR’ population structure depends markedly on ion charge, but can be significant even at electron densities as low as $`N_\mathrm{e}10^8\mathrm{cm}^3`$ typical of the solar corona. At much higher densities, $`N_\mathrm{e}10^{14}\mathrm{cm}^3`$, redistributive collisions can interrupt the two-step dielectronic process itself.

Therefore, we require dielectronic recombination from the metastable states as well as the ground. Secondly, we require final-state resolved data, i.e. we need to know the specific level that each two-step recombination ends-up in – the generalized collisional–radiative population rate equations govern the subsequent time-evolution of these states. (Incidentally, the first requirement imposes an additional requirement on the second, namely, that we now require dielectronic recombination into metastable parent final-states that lie above the ionization limit so that the collisional–radiative modelling recovers the Saha–Boltzman populations and ionization fractions in the high-density limit.) As we have already indicated, there are very many possible final states accessible to the dielectronic recombination process which would seem to make for unacceptably large tabulations. However, it is not only impracticable, but unnecessary, to treat each final state in the same manner for the purposes of collisional–radiative modelling. The techniques of matrix condensation and projection reduce the effective number of high-lying states by the progressive bundling of representative states and project the full influence of the high-lying states down onto a fully-resolved low-lying set. (It is the low-level set which is the focus of detailed spectral analysis.) With this in mind, we can tailor our dielectronic recombination tabulations to reflect this situation.

Of course, many people have calculated both partial and total dielectronic recombination rate coefficients and it is impractical to list them all here. A useful starting point is the compilation of total recombination rate coefficients from the literature by Mazzotta et al. (1998) who then used this data to compute the coronal ionization balance for all elements up to Ni. When the results of ab initio calculations are not available then much use, and abuse, is made of the General Formula of Burgess (1965). Partial dielectronic recombination rate coefficients also abound in the literature in connection with the study of particular physical problems: for example, at low temperatures, where only a few autoionizing states contribute significantly (Nussbaumer & Storey 1983); satellite lines (e.g. Bely-Dubau et al. 1979); laser produced plasmas, where the electron density and/or charge state is high enough that the non-LTE populations are concentrated in a limited number of low-lying states (see e.g. Abdallah & Clark 1994).

We have found it helpful for application to prepare and handle dielectronic recombination data in a hierarchy of increasing sophistication which we call baseline, level 1 and level 2. Baseline data are those produced using the Burgess general formula (GF), or with the techniques and state-selective programs associated with it. The generality of dielectronic data in use today are still from the GF. The background codes to the GF have a capability significantly beyond that of the GF, but are less well known. We return to these in more detail in section 2.3. The level 1 approach was introduced in the early ‘ninties’ to support the GCR modelling of light elements such as Be, B and C. These species are used in ‘light-element strategies’ for plasma facing wall components in fusion technology. For such elements, it is sufficient to use LS-coupled atomic structure and term population modelling. The relevant metastable populations in this case are terms such as C$`{}_{}{}^{2+}(2\mathrm{s}2\mathrm{p}^3\mathrm{P})`$. The level 2 approach, which is the main purpose of the present paper series, is concerned with the need to handle medium and heavy species in both fusion and astrophysics and to handle more extreme environments. In this approach, we aim to work with levels rather than terms and to use an intermediate coupling scheme, based-on the use of the Breit–Pauli Hamiltonian (Badnell 1997).

There are several reasons why the level 2 approach is now necessary: astrophysical spectral diagnostics tend to be based on levels rather than terms and the competition between autoionization and radiation makes it difficult to partition term-resolved dielectronic recombination data over levels; nuclear spin–orbit mixing is important even for low-$`Z`$ ions, e.g. carbon, because only a weak mixing of LS-forbidden autoionization rates with LS-allowed can give rise to ‘forbidden’ autoionization rates that are comparable with the dominant radiative rates (see Nussbaumer & Storey 1984, Badnell 1988); as $`Z`$ increases further then LS-forbidden radiative rates start to become significant; finally, dielectronic recombination via fine-structure transitions is completely absent in LS-coupling, giving rise to a large underestimate of the low-temperature dielectronic recombination rate coefficient in some iso-electronic sequences (see Savin et al. 1997).

The goal of this work is to calculate multi-configuration intermediate coupling dielectronic recombination rate coefficients from the (ground plus) metastable levels of an ion to all possible final states, resolved by level, and/or bundling, appropriate for generalized collisional–radiative modelling. We will cover elements applicable to astrophysics and magnetic fusion viz. He, Li, Be, B, C, N, O, F, Ne, Na, Mg, Al, Si, P, S, Cl, Ar, Ca, Ti, Cr, Fe, Ni, Zn, Kr, Mo and Xe. The first phase of the work will be the H- through Ne-like sequences. Level 1 LS-coupling data for many elements of these sequences was calculated by Badnell (1991–92, unpublished) and incorporated into the atomic database part of the Atomic Data and Analysis Structure (ADAS) and is routinely drawn into the generalized collisional–radiative part of ADAS (see Summers 2001). So, we already have a clear pathway through to the complete utilization of the detailed level 2 data that we will produce. The second phase of the work will cover the Na-like through Ar-like sequences, piloted initially by further (level 1) LS-coupling calculations to extend the 1991–92 data. Note that ADAS uses year numbers for the introduction of new approximations bringing substantive contributions to the database. The LS-coupled work of Badnell above has the year number ’93’. The third phase will focus on the remaining sequences of particular elements of interest, e.g. Fe. Even with the compactification of the partial dielectronic recombination data, along the lines already indicated, full publication in a paper journal is impractical, and not especially useful. So, the entire data will be made available via the World Wide Web (see section 5). The organization of the dielectronic data product follows the data format specifications of the ADAS Project. Dielectronic data are assigned to the data format adf09 (Summers 2001).

The plan for the remainder of the paper is as follows: in section 2.1 we review the generalized collisional–radiative approach encapsulated in ADAS and which influences our approach to handling dielectronic recombination data. In section 2.2 we describe and justify the theoretical approach that we take to calculate level 2 data. In section 2.3 we review in some detail the essence of the Burgess approach. It will be shown that this remains of importance for a full exploitation of the new work, for example, as applied to the $`l`$-redistribution of autoionizing states. Also, it is necessary to assess the progress in the precision of new dielectronic data in comparison with the baseline data. In section 3 we discuss the experimental situation for verifying dielectronic recombination data and the role of external fields. Some comparisons of results of our theoretical approach with high-resolution experimental results from storage rings are given. In section 4 we address derived data and we present some illustrative comparisons of GCR effective coefficients obtained using baseline, level 1 and level 2 data from dielectronic calculations. Also, we illustrate metastable-resolved ionization fractions. In section 5 we give more detail of the organization of the database and the computational implementation of its production. We finish with a short summary.

## 2 Theory

### 2.1 Generalized collisional–radiative modelling

Consider ions $`X^{+z}`$ of element $`X`$ of charge state $`z`$. We separate the levels of $`X^{+z}`$ into metastable levels $`X_\rho ^{+z}`$, indexed by Greek indices, and excited levels, indexed by Roman indices. The metastable levels include the ground level. We assume that the excited levels $`X_i^{+z}`$ are populated by excitation from all levels, $`\rho `$ and $`j`$, of $`X^{+z}`$, by ionization from the metastable levels of $`X_\mu ^{+z1}`$ and recombination from the metastable levels of $`X_\nu ^{+z+1}`$. The dominant population densities of these ions in the plasma are denoted by $`N_\rho ^z`$, $`N_\mu ^{z1}`$ and $`N_\nu ^{z+1}`$. The excited-state population densities, $`N_i^z`$, are assumed to be in quasi-static equilibrium with respect to the metastable populations. Thus,

$`0`$ $`=`$ $`{\displaystyle \underset{\sigma }{}}C_{i\sigma }^zN_\sigma ^z+{\displaystyle \underset{j}{}}C_{ij}^zN_j^z+{\displaystyle \underset{\mu }{}}S_{i\mu }^{z1}N_\mu ^{z1}+{\displaystyle \underset{\nu }{}}R_{i\nu }^{z+1}N_\nu ^{z+1},`$ (1)
where $`C_{ij}^z`$ are elements of the collisional–radiative matrix defined by

$`N_\mathrm{e}C_{ij}^z=A_{ji}^\mathrm{r}+N_\mathrm{e}q_{ji}^\mathrm{e},j>i,`$ (2)
where $`N_\mathrm{e}`$ is the electron density, $`q_{ji}^\mathrm{e}`$ is the electron-impact de-excitation rate coefficient and $`A_{ji}^\mathrm{r}`$ is the spontaneous radiative rate, both for the $`ji`$ transition in the ion $`X^{+z}`$. The equivalent expressions for upward transitions ($`j<i`$) and for $`i`$ and/or $`j`$ replaced by a metastable index follow trivially. The diagonal element $`C_{ii}^z`$ denotes the loss rate coefficient from the excited state $`i`$ and is given by

$`C_{ii}^z={\displaystyle \underset{ji}{}}C_{ji}^zS_i^z,`$ (3)
where

$`S_i^z={\displaystyle \underset{j}{}}S_{ji}^z`$ (4)
is the total ionization rate coefficient out of $`i`$. Finally, $`S_{i\mu }^{z1}`$ is the partial ionization rate coefficient out of metastable $`\mu `$ and $`R_{i\nu }^{z+1}`$ is the partial recombination rate coefficient out of metastable level $`\nu `$, both into level $`i`$ of the ion $`X^{+z}`$. $`S_{i\mu }^{z1}`$ includes contributions from both direct ionization and excitation-autoionization. $`R_{i\nu }^{z+1}`$ includes contributions from three-body, radiative and dielectronic recombination.

Solving for $`N_j^z`$, we have

$`N_j^z`$ $`=`$ $`{\displaystyle \underset{\sigma ,i}{}}\left(C^z\right)_{ji}^1C_{i\sigma }^zN_\sigma ^z{\displaystyle \underset{\mu ,i}{}}\left(C^z\right)_{ji}^1S_{i\mu }^{z1}N_\mu ^{z1}`$ (5)
$`{\displaystyle \underset{\nu ,i}{}}\left(C^z\right)_{ji}^1R_{i\nu }^{z+1}N_\nu ^{z+1}`$
$``$ $`{\displaystyle \underset{\sigma }{}}N_\mathrm{e}^\mathrm{X}_{j\sigma }^zN_\sigma ^z+{\displaystyle \underset{\mu }{}}N_\mathrm{e}^\mathrm{I}_{j\mu }^zN_\mu ^{z1}+{\displaystyle \underset{\nu }{}}N_\mathrm{e}^\mathrm{R}_{j\nu }^zN_\nu ^{z+1}`$
$``$ $`{\displaystyle \underset{\sigma }{}}{}_{}{}^{\mathrm{X}}N_{j\sigma }^{z}+{\displaystyle \underset{\mu }{}}{}_{}{}^{\mathrm{I}}N_{j\mu }^{z}+{\displaystyle \underset{\nu }{}}{}_{}{}^{\mathrm{R}}N_{j\nu }^{z},`$
where $`{}_{}{}^{\mathrm{X}}N_{j\sigma }^{z}`$, $`{}_{}{}^{\mathrm{I}}N_{j\mu }^{z}`$ and $`{}_{}{}^{\mathrm{R}}N_{j\nu }^{z}`$ are the effective populations of $`j`$ due to excitation, ionization and recombination from their respective metastables, and $`{}_{}{}^{\mathrm{X},\mathrm{I},\mathrm{R}}_{j\beta }^{z}`$ the corresponding coefficients. It is here that the connection with spectral analysis is made. The total emissivity in the line $`jk`$ is given by

$`\epsilon _{jk}^z`$ $`=`$ $`N_j^zA_{jk}^\mathrm{r},`$ (6)
with the populations $`N_j^z`$ given by equation (5). The corresponding photon emissivity coefficients are defined by

$`{}_{}{}^{\mathrm{X},\mathrm{I},\mathrm{R}}𝒫𝒞_{\beta ,𝒿𝓀}^𝓏`$ $``$ $`{}_{}{}^{\mathrm{X},\mathrm{I},\mathrm{R}}_{j\beta }^{z}A_{jk}^\mathrm{r},`$ (7)
respectively. Thus, the contributions to the total $`jk`$ emissivity from excitation, ionization and recombination are given by

$`{}_{}{}^{\mathrm{X},\mathrm{I},\mathrm{R}}\epsilon _{jk}^{z}`$ $`=`$ $`N_\mathrm{e}{\displaystyle \underset{\beta }{}}{}_{}{}^{\mathrm{X},\mathrm{I},\mathrm{R}}𝒫𝒞_{\beta ,𝒿𝓀}^𝓏𝒩_\beta ^{𝓏,𝓏\mathcal{1},𝓏+\mathcal{1}},`$ (8)
respectively.

The dynamic metastable populations $`N_\rho ^z`$ of $`X^z`$ satisfy

$`{\displaystyle \frac{1}{N_\mathrm{e}}}{\displaystyle \frac{\mathrm{d}N_\rho ^z}{\mathrm{d}t}}`$ $`=`$ $`{\displaystyle \underset{\sigma }{}}\left\{C_{\rho \sigma }^zN_\sigma ^z+{\displaystyle \underset{j}{}}C_{\rho j}^z{}_{}{}^{\mathrm{X}}N_{j\sigma }^{z}\right\}`$ (9)
$``$ $`{\displaystyle \underset{\sigma }{}}\left\{C_{\sigma \rho }^zN_\rho ^z+{\displaystyle \underset{j}{}}C_{\sigma j}^z{}_{}{}^{\mathrm{X}}N_{j\rho }^{z}\right\}`$
$`+`$ $`{\displaystyle \underset{\mu }{}}\left\{S_{\rho \mu }^{z1}N_\mu ^{z1}+{\displaystyle \underset{m}{}}S_{\rho m}^{z1}{}_{}{}^{\mathrm{X}}N_{m\rho }^{z1}\right\}`$
$``$ $`{\displaystyle \underset{\nu }{}}\left\{S_{\nu \rho }^zN_\rho ^z+{\displaystyle \underset{j}{}}S_{\nu j}^z{}_{}{}^{\mathrm{X}}N_{j\rho }^{z}\right\}`$
$`+`$ $`{\displaystyle \underset{\nu }{}}\left\{R_{\rho \nu }^{z+1}N_\nu ^{z+1}+{\displaystyle \underset{j}{}}C_{\rho j}^z{}_{}{}^{\mathrm{R}}N_{j\nu }^{z}\right\}`$
$``$ $`{\displaystyle \underset{\mu }{}}\left\{R_{\mu \rho }^zN_\rho ^z+{\displaystyle \underset{m}{}}C_{\mu m}^{z1}{}_{}{}^{\mathrm{R}}N_{m\rho }^{z1}\right\}`$
$`+`$ $`{\displaystyle \underset{\sigma }{}}\left\{Q_{\rho \sigma }^zN_\sigma ^z+{\displaystyle \underset{m}{}}S_{\rho m}^{z1}{}_{}{}^{\mathrm{R}}N_{m\sigma }^{z1}\right\}`$
$``$ $`{\displaystyle \underset{\sigma }{}}\left\{Q_{\sigma \rho }^zN_\rho ^z+{\displaystyle \underset{m}{}}S_{\sigma m}^{z1}{}_{}{}^{\mathrm{R}}N_{m\rho }^{z1}\right\}`$
$``$ $`{\displaystyle \underset{\sigma }{}}X_{\mathrm{CD}:\sigma \rho }^{zz}N_\sigma ^z{\displaystyle \underset{\sigma }{}}X_{\mathrm{CD}:\rho \sigma }^{zz}N_\rho ^z`$
$`+`$ $`{\displaystyle \underset{\mu }{}}S_{\mathrm{CD}:\mu \rho }^{z1z}N_\mu ^{z1}{\displaystyle \underset{\nu }{}}S_{\mathrm{CD}:\rho \nu }^{zz+1}N_\rho ^z`$
$`+`$ $`{\displaystyle \underset{\nu }{}}\alpha _{\mathrm{CD}:\nu \rho }^{z+1z}N_\nu ^{z+1}{\displaystyle \underset{\mu }{}}\alpha _{\mathrm{CD}:\rho \mu }^{zz1}N_\rho ^z`$
$`+`$ $`{\displaystyle \underset{\mu }{}}Q_{\mathrm{CD}:\sigma \rho }^{zz}N_\sigma ^z{\displaystyle \underset{\sigma }{}}Q_{\mathrm{CD}:\rho \sigma }^{zz}N_\rho ^z,`$
which defines the generalized collisional–radiative excitation ($`X_{\mathrm{CD}}`$), ionization ($`S_{\mathrm{CD}}`$), recombination ($`R_{\mathrm{CD}}`$) and parent metastable cross-coupling ($`Q_{\mathrm{CD}}`$) rate coefficients. (We note that $`Q_{\rho \sigma }^z0`$ initially.) This set of equations, together with those for $`N_j^z`$, are sufficient to solve the low-level problem – those levels with principal quantum number $`nn_\mathrm{c}`$, say. In the absence of dielectronic recombination, or at sufficiently high electron densities (e.g. $`10^{18}\mathrm{cm}^3`$, say, as found in laser-produced plasmas) then $`n_\mathrm{c}`$ can be small enough for this to be a complete solution since all higher levels are in collisional LTE, with a Boltzman population distribution. However, the presence of dielectronic recombination in low- to medium-density plasmas means that $`n_\mathrm{c}`$ can be prohibitively large ($`500`$, say).

We now consider a projection-condensation approach that allows for the effect of the high-level populations ($`\overline{n}n>n_\mathrm{c}`$) on the low-level populations ($`nn_\mathrm{c}`$). We work in the bundled-$`n`$ picture. Here the populations are grouped according to their parent level and principal quantum number. We assume that the high-level populations are in quasi-static equilibrium with the low-level populations and adjacent stage metastables. Thus, for each parent $`\tau `$, the high-level populations (denoted by $`\tau \overline{n}`$) satisfy

$`0`$ $`=`$ $`{\displaystyle \underset{n}{}}C_{\tau \overline{n},\tau n}^zN_{\tau n}^z+{\displaystyle \underset{\overline{n}^{}}{}}C_{\tau \overline{n},\tau \overline{n}^{}}^zN_{\tau \overline{n}^{}}^z`$ (10)
$`+`$ $`{\displaystyle \underset{\mu }{}}S_{\tau \overline{n},\mu }^{z1}N_\mu ^{z1}+{\displaystyle \underset{\nu }{}}R_{\tau \overline{n},\nu }^{z+1}N_\nu ^{z+1},`$
which is of the same form as equation (1) for the low-level excited-states. Thus, the low-level populations satisfy equations of the same form as (9). This is a full solution for all levels, in the bundled-$`n`$ picture. It includes direct couplings between the low-level populations (e.g. $`nn^{}`$) and indirect couplings via the high-level populations (e.g. $`n\overline{n}n^{}`$). However, we already have a description of the low-level problem in the fully-resolved picture, given by equations (1) and (9). We can supplement these equations with the indirect couplings of the bundled-$`n`$ picture, expanded over the low-level set using level weighting factors, $`\omega _{ij}`$. This projection corresponds to solving equations (1) and (9) with $`C`$, $`S`$, $`R`$ and $`Q`$ (including their appearance in equation (5)) replaced by $`𝒞`$, $`𝒮`$, $``$ and $`𝒬`$ where $`𝒞=C+^{\mathrm{ind}}C`$, $`𝒮=S+^{\mathrm{ind}}S`$, $`=R+^{\mathrm{ind}}R`$ and $`𝒬=Q+^{\mathrm{ind}}Q`$, where

$`{}_{}{}^{\mathrm{ind}}C_{\rho \sigma }^{z}`$ $`=`$ $`\omega _{\rho ,\tau n}\omega _{\tau n^{},\sigma }{\displaystyle \underset{\overline{n},\overline{n}^{}}{}}C_{\tau n,\tau \overline{n}^{}}^z\left(C^z\right)_{\tau \overline{n}^{},\tau \overline{n}}^1C_{\tau \overline{n},\tau n^{}}^z,`$ (11)
$`{}_{}{}^{\mathrm{ind}}S_{\nu ^{}i}^{z}`$ $`=`$ $`\omega _{\nu n,i}{\displaystyle \underset{\overline{n},\overline{n}^{}}{}}S_{\nu ^{},\nu \overline{n}^{}}^z\left(C^z\right)_{\nu \overline{n}^{},\nu \overline{n}}^1C_{\nu \overline{n},\nu n}^z,`$ (12)
$`{}_{}{}^{\mathrm{ind}}R_{i\nu ^{}}^{z+1}`$ $`=`$ $`\omega _{i,\nu n}{\displaystyle \underset{\overline{n},\overline{n}^{}}{}}C_{\nu n,\nu \overline{n}^{}}^z\left(C^z\right)_{\nu \overline{n}^{},\nu \overline{n}}^1R_{\nu \overline{n},\nu ^{}}^{z+1}`$ (13)
and

$`{}_{}{}^{\mathrm{ind}}Q_{\rho \sigma }^{z}`$ $`=`$ $`\omega _{\rho ,\tau n}\omega _{\tau n^{},\sigma }{\displaystyle \underset{\overline{n},\overline{n}^{}}{}}S_{\rho ,\tau \overline{n}^{}}^{z1}\left(C^{z1}\right)_{\tau \overline{n}^{},\tau \overline{n}}^1R_{\tau \overline{n},\sigma }^z.`$ (14)
Although we now have a complete solution in terms of the fully-resolved low-level and bundled-$`n`$ high-level picture, one further step is of practical significance. In order to span a wide range of electron densities it is necessary to treat very large principal quantum numbers in order to reach the collision limit at low densities. It is not necessary to treat each $`n`$ individually. Rather, a set of representative $`n`$-values can be used instead. If $`N_{\overline{n}}`$ denotes the bundled-$`n`$ populations for $`\overline{n}=n_\mathrm{c}+1,n_\mathrm{c}+2,\mathrm{}`$, and $`𝒩_{\overline{\overline{n}}}`$ denotes a subset of them, then the two are related via $`N_{\overline{n}}=\omega _{\overline{n}\overline{\overline{n}}}𝒩_{\overline{\overline{n}}}`$, where $`\omega _{\overline{n}\overline{\overline{n}}}`$ are the interpolation coefficients. Substituting for $`N_{\overline{n}}`$ into equations (1) and (9) yields a condensed set of equations for $`𝒩_{\overline{\overline{n}}}`$ and the $`𝒞`$, $`𝒮`$, $``$ and $`𝒬`$ obtained from the condensed set of equations are identical in form to those obtained from the full set of equations.

We note here that we have made an assumption viz. that, following dielectronic capture, the autoionizing state is not perturbed by a further collision before it either autoionizes or radiates. This is not due to a limitation of ADAS but rather a choice that we have made (and defined within the adf09 specification) so as to make the general collisional–radiative problem tractable over a wide range of electron densities. Working explicitly with autoionization and radiative rates and bound and non-bound states rather than partial dielectronic recombination rate coefficients and (mostly) bound states vastly increases the data requirements, in general. Although our collisional–radiative model goes over to the correct LTE limit at high electron densities, there is a density range ($`10^{16}`$ cm<sup>-3</sup>, found in laser-produced plasmas) where the levels of spectroscopic interest have non-LTE populations that are influenced by non-LTE populations of autoionizing levels that are themselves collisionally redistributed. We describe an approximate solution in section 2.3 below.

We are now in a position to spell out our requirements of recombination data:
(i) We require recombination data from all metastable levels, not just the ground.
(ii) We require recombination data into particular final states.
(iii) We require the final-state to be level-resolved for $`nn_\mathrm{c}`$ and parent-level-resolved bundled-$`n`$ for $`n>n_\mathrm{c}`$.
(iv) Parent metastable cross-coupling means that we require recombination into metastable autoionizing final states.
(v) We need only produce data for a representative set of $`\overline{n}`$.
The ADAS adf09 data specification incorporates all of these requirements. Finally, we note that use of total zero-density ground-state recombination rate coefficients is quite unsafe for the collisional–radiative modelling of dynamic finite-density plasmas.

### 2.2 Dielectronic recombination rate coefficient modelling

We have already noted that the partial recombination rate coefficient ($`R_{i\nu }^{z+1}`$) includes contributions from three-body, radiative and dielectronic recombination. In ADAS, three-body recombination rate coefficients are obtained from electron-impact ionization rate coefficients, via detailed balance. This also ensures that the correct Saha–Boltzman limit is reached at high electron densities. Since three-body recombination is separate from dielectronic and radiative recombination, it is not necessary to consider if further. However, quantum mechanically, dielectronic and radiative recombination are indistinguishable processes which interfere with each other. In practice (see Pindzola et al. 1992), this interference is a very small effect and can safely be neglected for our purposes. This is the independent processes approximation whereby dielectronic and radiative recombination can be considered separately and is the approach taken by the database aspect of ADAS. Separate data files exist for dielectronic (adf09) and radiative recombination (adf08)and they can be updated independently. Our focus is dielectronic recombination. Details of the ADAS data status for radiative recombination can be found in Summers (2001).

In the isolated resonance approximation, the partial dielectronic recombination rate coefficient $`\alpha _{i\nu }^{z+1}`$ from an initial metastable state $`\nu `$ into a resolved final state $`i`$ of an ion $`X^{+z}`$ is given by

$`\alpha _{i\nu }^{z+1}`$ $`=`$ $`\left({\displaystyle \frac{4\pi a_0^2I_\mathrm{H}}{k_\mathrm{B}T_\mathrm{e}}}\right)^{3/2}{\displaystyle \underset{j}{}}{\displaystyle \frac{\omega _j}{2\omega _\nu }}\mathrm{e}^{E_c/k_\mathrm{B}T_\mathrm{e}}`$ (15)
$`\times `$ $`{\displaystyle \frac{\underset{l}{}A_{j\nu ,E_cl}^\mathrm{a}A_{ji}^\mathrm{r}}{_hA_{jh}^\mathrm{r}+_{m,l}A_{jm,E_cl}^\mathrm{a}}},`$
where $`\omega _j`$ is the statistical weight of the $`(N+1)`$-electron doubly-excited resonance state $`j`$, $`\omega _\nu `$ is the statistical weight of the $`N`$-electron target state and the autoionization ($`A^\mathrm{a}`$) and radiative ($`A^\mathrm{r}`$) rates are in inverse seconds. Here, $`E_c`$ is the energy of the continuum electron (with angular momentum $`l`$), which is fixed by the position of the resonances, and $`I_\mathrm{H}`$ is the ionization potential energy of the hydrogen atom (both in the same units of energy), $`k_\mathrm{B}`$ is the Boltzman constant, $`T_\mathrm{e}`$ the electron temperature and $`(4\pi a_0^2)^{3/2}=6.6011\times 10^{24}`$ cm<sup>3</sup>. The effect of interacting resonances on dielectronic recombination has been investigated by Pindzola et al. (1992) and can safely be neglected, at least in the absence of external electric and magnetic fields (see section 3 below). While autoionization rates can be determined (within the isolated resonance approximation) via the fitting of resonances calculated in a close-coupling approximation, or via the extrapolation of threshold close-coupling collision strengths using the correspondence principle, it is usual now to introduce a further approximation – that of using distorted waves, i.e. the autoionization rates are calculated via perturbation theory using the Golden Rule (Dirac 1930). This is the only approximation that we have made so far that may need to be reconsidered in certain cases. In low-charge ions, a perturbative distorted wave calculation may give inaccurate autoionization rates compared to those calculated in a close-coupling approximation. However, this only has a direct effect on the partial dielectronic recombination rate if the autoionization rates do not ‘cancel-out’ between the numerator and denominator of (15) – typically, autoionization rates are orders of magnitude larger than radiative rates.

One could obtain (some) partial dielectronic recombination data from an $`R`$-matrix photoionization calculation, on making use of detailed balance, either in the absence of radiation damping (Nahar & Pradhan 1994) or with its inclusion (Robicheaux et al. 1995, Zhang et al. 1999). (One must take care not to double count the radiative recombination contribution in the modelling now.) However, this cannot provide us with a complete set of partial recombination rate coefficients since it is only possible to compute photoionization from (i.e. photorecombination to) a relatively low-lying set of states – up to $`n10`$, say. Total recombination rates are obtained by supplementing the photorecombination data with high-$`n`$ ‘close-coupling’ dielectronic recombination rate coefficients calculated using Bell & Seaton (1985) or Hickman’s (1984) approach. This is based on the radiative-loss term from a unitary S-matrix and does not, and cannot, resolve recombination into a particular final state, which is essential for collisional–radiative modelling. Indeed, even when summed-over all final states, errors can still result. This has been demonstrated explicitly by Gorczyca et al. (2002) in the case of Fe<sup>17+</sup>. They found that only the IPIRDW approach could reproduce the measured dielectronic recombination cross section of Savin et al. (1997, 1999) for high Rydberg states. Thus, our initial goal is to generate complete data sets within the independent processes and isolated resonance using distorted waves (IPIRDW) approximation. Subsequently, selective upgrades from $`R`$-matrix data may be made via, for example, the RmaX network which can be viewed as a progression of the Iron Project (Hummer et al. 1993) and which is focusing on X-ray transitions – see, for example, Ballance et al. (2001). We note that while the Opacity Project (Seaton 1987) calculated a large amount of photoionization data, which in principle could be used for recombination (via detailed balance), unfortunately, only total photoionization cross sections were archived, i.e. summed-over the final electron continuum, and so it is impossible to apply detailed balance and so it cannot be used as a source of recombination data.

We use the code autostructure (Badnell 1986, Badnell & Pindzola 1989, Badnell 1997) to calculate multi-configuration intermediate coupling energy levels and rates within the IPIRDW approximation. The code can make use both of non-relativistic and semi-relativistic wavefunctions (Pindzola and Badnell 1990). The low-$`n`$ problem is no different from the one of computing atomic structure. The high-$`n`$ problem requires some discussion. The mean radius of a Rydberg orbital scales as $`n^2`$ and so it rapidly becomes impossible to calculate an explicit bound orbital (for $`n>20`$, say) and some approximation must be made. We note that (Seaton 1983)

$`\underset{n\mathrm{}}{lim}\left({\displaystyle \frac{\pi n^3}{2z^2}}\right)^{3/2}P_{nl}(r)=F_{kl}(r)|_{k=0},`$ (16)
where the bound orbitals are normalized to unity and the continuum orbitals to $`\pi \delta (kk^{})`$, here $`k^2=E_c`$(Ry). The approach taken in autostructure is to make use of (16) at finite $`n`$, i.e. to approximate the bound orbital by a suitably normalized zero-energy continuum orbital, for $`n>15+l^2/4`$ (evaluated in integer arithmetic). A further refinement is to evaluate the (true) continuum orbital for the incident electron at $`E_c+z^2/n^2`$ Rydbergs (instead of $`E_c`$) so as to maintain the same transition energy. (In the Bethe approximation, the free–free dipole acceleration integral is proportional to $`\mathrm{\Delta }\epsilon ^2`$ times a slowly-varying-with-energy dipole-length integral, where $`\mathrm{\Delta }\epsilon `$ is the transition energy.) Actually, the true continuum orbital is calculated at about 15 energies per $`l`$ and the one- and two-body bound–free integrals are interpolated at the required energy, which is given by energy conservation. (The use of a zero-energy continuum orbital means that long-range free–free integrals arise and these are treated using the techniques of distorted wave scattering theory, see Badnell 1983.) For each $`nl`$, autostructure reforms both the $`N`$\- and $`(N+1)`$-electron Hamiltonians and diagonalizes them (separately) to reform the rates. Only for H-like ions have we found it necessary to treat all $`l`$ at the same time, and treat only each $`n`$ separately. Typically, each $`n`$ is calculated explicitly until no new continua can open-up and then only a representative set of $`n`$, up to $`n=999`$, is used. This approach avoids the extrapolation of low-$`n`$ autoionization rates or, even worse, partial dielectronic recombination rate coefficients to high-$`n`$. Use of (16) is, in effect, an interpolation since it is exact in the limit $`n\mathrm{}`$ and any error is bounded at the lowest-$`n`$ by knowledge of the ‘exact’ result obtained from using $`P_{nl}`$ directly, rather than $`F_{kl}`$.

autostructure is implemented within ADAS as ADAS701. It produces the raw autoionization and radiative rates. To produce partial dielectronic recombination rate coefficients, according to the prescription of Sect. 2.1, requires further non-trivial organization of the raw data. In particular, radiative transitions between highly-excited Rydberg states are computed hydrogenically and added-in during a ‘post-processing’ exercise with the code adasdr, which is implemented with ADAS as ADAS702. Also, observed energies for the core and parent levels are used at this stage to ensure accurate positioning of the resonances and, hence, accurate low-temperature rate coefficients. adasdr outputs directly the adf09 file for use by ADAS. Separate adf09 files are produced for different ‘core-excitations’ ($`nn^{}`$), e.g. $`12`$, $`22`$ and $`23`$ for Li-like ions. This enables selective upgrades of the adf09 database.

### 2.3 Exploitation of the Burgess–Bethe approach

We are concerned with how the precise calculations described above relate to other calculations and, in particular, to those commonly used in astrophysics. Our baseline calculation is based on the methodologies of Burgess, which represent what can be achieved without recourse to the detail of the above sections. The Burgess GF itself was in fact a functional fit to extended numerical calculations. The associated code, with extensions, we call the ‘Burgess–Bethe general program’ (BBGP). It will be shown in this sub-section how the BBGP can be used to obtain a working model for the $`l`$-redistribution of doubly-excited states and, hence, provide a correction to accurate, but unredistributed, dielectronic data so as to model the dynamic part of the plasma microfield. Also, our baseline, calculated using the BBGP, will allow an assessment of the typical error that is present in the general dielectronic modelling in astrophysics to date.

In the LS-coupled term picture, introduce a set $`P`$ of parent terms $`\gamma _p^{2S_p+1}L_p`$ of energy $`E_p`$ relative to the ground parent term, indexed by $`p`$. Suppose that the excited parents are those with $`p_{min}pP`$. The metastable parents, which are the initial metastables for recombination and the final parents on which the recombining excited $`nl`$-electron is built, are the subset $`1p<p_{min}=M`$. Let $`p^{}`$ denote an initial parent with the incident electron denoted by $`k^{}l^{}`$. We wish to re-establish the expressions used by Burgess in his development and it is helpful to work in $`z`$-scaled dimensionless coordinates. Then, introducing $`z_{\mathrm{eff}}`$, the effective charge, the collision strength for a dipole excitation of $`\gamma _p^{2S_p+1}L_p`$ to $`\gamma _p^{}{}_{}{}^{2S_p+1}L_{p}^{}`$, evaluated in the Bethe approximation, is given by

$`\mathrm{\Omega }((S_pL_p^{})k^{}l^{},(S_pL_p)klSL)=48((z_{\mathrm{eff}}+1)^2/z_{\mathrm{eff}}^4)`$
$`\times \left({\displaystyle \frac{I_\mathrm{H}}{\mathrm{\Delta }ϵ_{pp^{}}}}\right)\left({\displaystyle \frac{(2S+1)(2L+1)}{2(2S_p+1)}}\right)l_>\left\{\begin{array}{ccc}L_p^{}& L_p& 1\\ l& l^{}& L\end{array}\right\}^2`$ (19)
$`\times (2S_p+1)(2L_p^{}+1)f_{S_pL_p^{}S_pL_p}\left|F_{k^{}l^{}}\left|\rho ^2\right|F_{kl}\right|^2`$ (20)
in terms of the parent oscillator strength $`f`$, where $`l_>=\mathrm{max}(l,l^{})`$. For dielectronic recombination, it is convenient to express this in terms of the Einstein A-coefficient for the parent transition and to analytically continue the collision strength to negative energies for the $`kl`$ electron, that is as $`\kappa =k/z_{\mathrm{eff}}i/n`$. Formally, we identify a band of free-electron energies $`\mathrm{d}E/I_\mathrm{H}=z_{\mathrm{eff}}^2\mathrm{d}ϵ/I_\mathrm{H}`$, with the separation between $`n`$-shells $`2z_{\mathrm{eff}}^2I_\mathrm{H}/n^3`$, and let

$`\left|F_{k^{}l^{}}\left|\rho ^2\right|F_{kl}\right|^2\mathrm{d}(ϵ/I_\mathrm{H})(\pi /16)(ϵ/I_\mathrm{H})^4\left|F_{k^{}l^{}}\left|\rho \right|P_{nl}\right|^2.`$ (21)
Then, the (dielectronic) resonance-capture cross section is given by

$`Q^\mathrm{c}((S_pL_p^{})k^{}l^{},(S_pL_p)klSL)\mathrm{d}(E/I_\mathrm{H})\left({\displaystyle \frac{(z_{\mathrm{eff}}+1)^2}{z_{\mathrm{eff}}^4}}\right)`$
$`\times \left({\displaystyle \frac{6\pi ^2a_0^3}{\alpha ^4c}}\right)\left({\displaystyle \frac{\mathrm{\Delta }ϵ_{pp^{}}}{I_\mathrm{H}}}\right)\left({\displaystyle \frac{I_\mathrm{H}}{ϵ}}\right)\left({\displaystyle \frac{(2S+1)(2L+1)}{2(2S_p+1)}}\right)l_>`$
$`\times \left\{\begin{array}{ccc}L_p^{}& L_p& 1\\ l& l^{}& L\end{array}\right\}^2\left({\displaystyle \frac{A^\mathrm{r}(S_pL_pS_pL_p^{})}{(z_{\mathrm{eff}}+1)^4}}\right)\left|F_{k^{}l^{}}\left|\rho \right|P_{nl}\right|^2.`$ (24)
The inverse (Auger) rate coefficient is obtained by invoking detailed balance as

$`A^\mathrm{a}((S_pL_p)nlSL(S_pL_p^{})k^{}l^{})=\left({\displaystyle \frac{(z_{\mathrm{eff}}+1)^2}{z_{\mathrm{eff}}^2}}\right)\left({\displaystyle \frac{3}{2\alpha ^3}}\right)`$
$`\times \left({\displaystyle \frac{\mathrm{\Delta }ϵ_{pp^{}}}{I_\mathrm{H}}}\right)(2L_p^{}+1)l_>\left\{\begin{array}{ccc}L_p^{}& L_p& 1\\ l& l^{}& L\end{array}\right\}^2\left({\displaystyle \frac{A^\mathrm{r}(S_pL_pS_pL_p^{})}{(z_{\mathrm{eff}}+1)^4}}\right)`$ (27)
$`\times \left|F_{k^{}l^{}}\left|\rho \right|P_{nl}\right|^2.`$ (28)
For the generalized collisional–radiative modelling of light element ions, it is convenient to use $`L`$-averaged doubly-excited levels, but still resolved by spin $`S`$, whereas the BBGP treatment uses $`SL`$-averaged levels. The corresponding Auger rates are then

$`A^\mathrm{a}((S_pL_p)nlS(S_pL_p^{})k^{}l^{})=\left({\displaystyle \frac{(z_{\mathrm{eff}}+1)^2}{z_{\mathrm{eff}}^2}}\right)\left({\displaystyle \frac{1}{2\alpha ^3}}\right)\left({\displaystyle \frac{\mathrm{\Delta }ϵ_{pp^{}}}{I_\mathrm{H}}}\right)`$
$`\times \left({\displaystyle \frac{l_>}{(2l+1)}}\right)\left({\displaystyle \frac{A^\mathrm{r}(S_pL_pS_pL_p^{})}{(z_{\mathrm{eff}}+1)^4}}\right)\left|F_{k^{}l^{}}\left|\rho \right|P_{nl}\right|^2,`$ (29)
in both cases. The matching resonance-capture coefficients are obtained by detailed balance, or by summing and averaging over the resolved expression for $`Q^\mathrm{c}`$ given above.

Turning to the radiative decay of the doubly-excited resonant states in the $`LS`$-resolved picture: the spontaneous emission coefficient, with a passive spectator in the $`nl`$ shell, is given by

$`A^\mathrm{r}((S_pL_p)nlSL(S_pL_p^{})nlSL^{})=(2L_p+1)(2L^{}+1)`$
$`\times \left\{\begin{array}{ccc}L_p^{}& L_p& 1\\ L& L^{}& l\end{array}\right\}^2A^\mathrm{r}(S_pL_pS_pL_p^{}),`$ (32)
where $`A^\mathrm{r}(S_pL_pS_pL_p^{})`$ is the parent-core spontaneous transition probability. The $`L`$\- and $`LS`$-averaged probabilities are both simply equal to $`A^\mathrm{r}(S_pL_pS_pL_p^{})`$. The BBGP method exploits the fact that, in the dipole case,

$`{\displaystyle \frac{A^\mathrm{a}}{A^\mathrm{r}}}`$ $`=`$ $`{\displaystyle \frac{l_>\left|F_{k^{}l^{}}\left|\rho \right|P_{nl}\right|^2}{2\alpha ^3(2l+1)(z_{\mathrm{eff}}+1)^2z_{\mathrm{eff}}^2}},`$ (33)
and efficient recurrence relations are available for the generation of hydrogenic bound–free radial integrals for all parameter values. It is clear that the Bethe approximation for the partial collision strengths can be substantially in error for $`0l2`$. Burgess introduced correction factors for the lowest partial collision strengths, based-on a comparison with more sophisticated collision calculation results that were available at the time. More precisely, introduce

$`\mathrm{𝑐𝑜𝑟}_l`$ $`=`$ $`{\displaystyle \underset{l^{}}{}}\mathrm{\Omega }((S_pL_p^{})k^{}l^{},(S_pL_p)k|_{=0}l)/`$ (34)
$`{\displaystyle \underset{l^{}}{}}\mathrm{\Omega }^{\mathrm{Bethe}}((S_pL_p^{})k^{}l^{},(S_pL_p)k|_{=0}l).`$
The general formula for zero-density total dielectronic recombination rate coefficients used a fixed-set of $`\mathrm{𝑐𝑜𝑟}_l`$ for all parent transitions. These were based-on parent $`1\mathrm{s}2\mathrm{p}`$ transitions and it is the case that the inclusion of corrections is most significant for parent $`\mathrm{\Delta }n1`$ transitions. To exploit the BBGP method beyond its use for the general formula, we must establish the population equations of the doubly-excited levels. For $`LS`$-averaged levels, the number densities expressed in terms of their deviations, $`b_{p,nl}`$, from Saha–Boltzmann, and referred to the initial parent $`p^{}`$, are given by

$`N_{p,nl}=N_\mathrm{e}N_p^{}^+8\left[{\displaystyle \frac{\pi a_0^2I_\mathrm{H}}{k_\mathrm{B}T_\mathrm{e}}}\right]^{3/2}{\displaystyle \frac{\omega _{p,nl}}{\omega _p^{}}}e^{E/kT_\mathrm{e}}b_{p,nl}.`$ (35)
Then, in the BBGP baseline zero-density limit, with only resonant capture from the $`p^{}`$ parent balanced by Auger breakup and radiative stabilization back to the same parent, we have

$`b_{p,nl}`$ $`=`$ $`\left({\displaystyle \frac{_l^{}A^\mathrm{a}(p,nlp^{}k^{}l^{})}{_l^{}A^\mathrm{a}(p,nlp^{}k^{}l^{})+A^\mathrm{r}(p,nlp^{},nl)}}\right).`$ (36)
In the extended BBGP program, we can also include resonant-capture from initial metastables other than the ground, dipole-allowed collisional redistribution between adjacent doubly-excited $`l`$-substates by secondary ion- and electron-impact, and losses by ‘alternate’ Auger break-up and parent radiative transition pathways. The population equations for the $`l`$-substates of a doubly-excited $`n`$-shell become

$`\left(N_\mathrm{e}q_{nl1nl}^\mathrm{e}+N^{z_{\mathrm{eff}}}q_{nl1nl}^{z_{\mathrm{eff}}}\right)N_{p,nl1}`$
$`+({\displaystyle \underset{l^{}=l\pm 1}{}}N_\mathrm{e}q_{nlnl^{}}^\mathrm{e}+{\displaystyle \underset{l^{}=l\pm 1}{}}N^{z_{\mathrm{eff}}}q_{nlnl^{}}^{z_{\mathrm{eff}}}`$
$`+{\displaystyle \underset{p_1=1}{\overset{p1}{}}}{\displaystyle \underset{l^{}=l1}{\overset{l+1}{}}}A_{p,nlp_1,\kappa l^{}}^\mathrm{a}+{\displaystyle \underset{p_1=1}{\overset{p1}{}}}A_{p,nlp_1,nl}^\mathrm{r})N_{p,nl}`$
$`\left(N_\mathrm{e}q_{nl+1nl}^\mathrm{e}+N_{z_{\mathrm{eff}}}q_{nl+1nl}^{z_{\mathrm{eff}}}\right)N_{p,nl+1}`$
$`=N_\mathrm{e}{\displaystyle \underset{p_2=1}{\overset{M}{}}}{\displaystyle \underset{l^{}=l1}{\overset{l+1}{}}}q_{p_2,\kappa l^{}p,nl}^\mathrm{c}N_{p_2}+{\displaystyle \underset{p_1=p+1}{\overset{P}{}}}A_{p_1,nlp,nl}^\mathrm{r}N_{p_1,nl}.`$ (37)
These equations may be solved progressively downwards through the levels built on excited parents, terminating with levels built on the ‘final’ ground and metastable parents. The calculations yield state-selective dielectronic recombination coefficients to levels built on each metastable parent, together with Auger rates, from levels built on metastables above the ground parent, to lower metastable parents. The results from the above solution, at zero density, we call baseline data.

We note some details of the implementation:

#### 2.3.1 Collisional rates

Ion and electron $`l`$-redistributive cross sections are evaluated following the method of Pengelly & Seaton (1964). This is to be viewed as a very simplified treatment of the dynamic part of the plasma microfield associated with ions which move closer than their neighbours to the target. The quasi-static part of the ion microfield, which is approximately 1/3 of the total, is ignored at this level of analysis (see section 3 for further discussion of field effects in an experimental context).

#### 2.3.2 Energy levels

The energy differences $`E_{p,nlS}E_{p,nl\pm 1S}`$ are critical to the problem of doubly-excited-state redistribution. They are small, tending to zero at large $`l`$. Thus, the redistributive cross-sections are very large, remaining finite in the degenerate level limit only because of the finite radiative lifetime of the target (which is short for resonant states) or through screening of the projectile by nearest neighbours. For the $`l=0,1,2`$ waves of the $`n`$-shell spectator electron, quantum defect expansions of the form $`\mu _{p,nlS}=a_0+a_1/n^2`$ are usually available. For the higher $`l(>2)`$, the large number of energies required can be estimated more economically from the dipole polarizabilities, $`\alpha _p^{\mathrm{pol}}`$, of the parents, following Edlen(1964), as

$`E_{p,nl}`$ $`=`$ $`(z_1/n)^2[1+5.23504^4(z_1/n)^2(n/(l+1/2)+`$ (38)
$`3/4)]+\alpha ^{\mathrm{pol}}_p(z_1/n)^4(3n^2l(l+1))/`$
$`(2n(l1/2)l(l+1/2)(l+1)(l+3/2)).`$
Small-scale autostructure runs have been used to prepare these data for the baseline. We use two-term quantum defect expansions fitted at $`n\mathrm{}`$ and $`n=10`$.

#### 2.3.3 Bethe correction factors and radiative transition probabilities

For our comparisons, it is important that the inputs for the BBGP baseline calculations are consistent with the autostructure calculations for the level 1 and level 2 results. Small-scale autostructure runs have been used to prepare these inputs for the baseline. Low partial-wave collision strengths are now widely available (e.g. from autostructure), such that the $`\mathrm{𝑐𝑜𝑟}_l`$ can be prepared fairly easily. Within the ADAS Project, efficient subroutines for bound–free integrals and complete BBGP calculations are available together with sets of $`\mathrm{𝑐𝑜𝑟}_l`$ for the principal types of transitions. They are available to those for whom direct utilization of level 1 and level 2 data in full GCR modelling is not an option.

The results presented in this overview paper are illustrative only. Fig. 1 contrasts zero-density total dielectronic recombination coefficients ($`\alpha _\mathrm{d}^{\mathrm{tot}}`$) calculated in the GF and the BBGP baseline approximations with those of the level 1 and level 2 computations reported here. In particular, we note that level 1 and 2 data are required to describe the recombination at low temperatures and that the level 2 data provides a noticeable refinement over the level 1 results. Figs. 2a,b illustrate the partial recombination into $`n`$-shells and the population structure of the $`l`$-subshells of a representative doubly-excited $`n`$-shell. Fig. 2a shows the very good convergence of BBGP to level 1 data with increasing completeness of alternate Auger pathways. Fig. 2b shows the effects of collisional redistribution at finite-density. A ratio of the sum over $`l`$-substates at a given density to that at zero density yields a BBGP finite-density adjustment factor of the total $`n`$-shell capture at zero density. The consistency between the BBGP, level 1, and level 2 approaches allows us to use this adjustment factor on the level 1 and level 2 data. In advanced generalized collisional–radiative modelling, the BBGP finite-density redistributive code acts as an interface between the extraction of state-selective zero-density dielectronic data from the ADAS adf09 database and its entry into the GCR population codes, corrected for finite-density doubly-excited state redistribution. Note also that routine semi-automatic comparisons as, illustrated here, provide the theoretical uncertainty estimate with which we can tag each dielectronic datum.

## 3 Experimental validation of dielectronic recombination data and the role of fields

The last decade has seen an enormous amount of experimental activity in the area of dielectronic recombination. In particular, heavy-ion storage rings coupled with electron-coolers have provided a wealth of data for partial dielectronic recombination. (The partial here is by the intermediate resonance state rather than the final state.) Most data is of the ‘bundled-$`n`$’ form, but some $`l`$-resolution is possible for very low-lying states. The iso-electronic sequences studied range almost exclusively from H-like through to Na-like and nuclear charges have ranged between $`Z`$=2 and 92. In all cases, in the absence of external fields, there is rarely any significant disagreement with theory, i.e. outside of the experimental uncertainty. A few of the more recent, typical, comparisons between experiment and the results of IPIRDW calculations include: Böhm et al. (2002), Savin et al. (2002a,b), Brandau et al. (2002). In Fig. 3, we show representative comparisons of dielectronic recombination data for O$`{}_{}{}^{5+}+\mathrm{e}^{}`$O<sup>4+</sup> and Fe$`{}_{}{}^{18+}+\mathrm{e}^{}`$Fe<sup>17+</sup>, calculated in the IPIRDW approximation with autostructure and which illustrates the level of accuracy that can be expected of the theoretical data.

One major area of uncertainty is the role of external fields on dielectronic recombination, and it is this more than anything that renders pointless efforts to compute (zero-density) field-free data to an accuracy of better than $``$20% , say. It has long been known that the high Rydberg states that frequently dominate the dielectronic recombination process can be Stark-mixed by weak electric fields (Burgess & Summers 1969), in particular the plasma microfield (Jacobs et al. 1976), and so increase the partial rate coefficients by factors of 2, or 3, or more, over a wide range of $`n`$. Recently, the picture has been further complicated by the discovery that magnetic fields, when crossed with an electric field, strongly affect the electric field enhancement – by reducing it in most cases (see Robicheaux et al. 1997, Bartsch et al. 1999, Schippers et al. 2000 and Böhm et al. 2001). While this suppression of the electric field enhancement is advantageous towards the use of field-free dielectronic recombination data, it is disadvantageous in terms of trying to compute field-dependent data for plasma modelling. Previously, it appeared that a reasonable approach would be to use the values of the plasma microfield (which in turn depends on the plasma density) for the electric field strength for use in the generation of field dependent (i.e. density dependent) data as input to collisional–radiative modelling. This in itself ignored any further (e.g. external) electric fields that might be present in the plasma environment, beyond the plasma microfield. The recognition of the importance of magnetic fields as well makes a comprehensive solution to dielectronic recombination in a plasma a distant goal and partial data accurate to $``$20% as meaningful as necessary. Furthermore, field enhancement is sensitive to interacting resonances as well (see Robicheaux et al. 1998) unlike the field-free case. We do note again that high Rydberg states in a finite density plasma are brought into LTE by (electron) collisions. A preliminary study by Badnell et al. (1993) showed that the effect of the plasma microfield on the density-dependent effective recombination rate coefficient was suppressed by collisions driving high Rydberg states into LTE – larger values for the microfield, which lead to larger enhancements of the zero-density rate coefficient, corresponds to denser plasmas for which collisions drive more states into LTE.

## 4 Comparisons and assessments of the derived data

As described in section 2.1, the emissivity and generalized collisional–radiative (GCR) coefficients depend inter alia on the fundamental dielectronic cross section data. There are two issues of concern in assessing these derived theoretical data, viz. the relative contributions to the effective coefficients coming from many different direct and indirect pathways and, secondly, estimation of the uncertainty in the theoretical data, which may be treated as a ‘working error’ in the interpretation of spectral observations from plasmas.

For the effective photon emissivity coefficients ($`𝒫𝒞`$s), it is firstly to be noted that the relative importance of the contribution from excitation ($`{}_{}{}^{\mathrm{X}}𝒫𝒞`$) and recombination ($`{}_{}{}^{\mathrm{R}}𝒫𝒞`$) is directly proportional to the ionization balance fractional abundances of the (metastable) ‘driver’ populations. The recombination part is most significant in transiently recombining plasmas and it is on this part only that we focus here. The partitioning of the collisional–radiative matrix described in section 2.1 allows us to contrast the direct capture, capture coming via the complete set of resolved low-levels and capture via the high bundled-$`n`$ quantum shells, which are treated by projection. The relative contributions depend differentially on density since the projection part is suppressed selectively at higher densities. Also, electron temperature and the recombining ion charge influence the relative importance of the dielectronic and radiative recombination contributions and the role of the more highly-excited levels. In Fig. 4, we show the main effects with some illustrative results from ADAS for the C iii $`2\mathrm{s}2\mathrm{p}^3\mathrm{P}2\mathrm{s}3\mathrm{d}^3\mathrm{D}`$ multiplet at 459.6 Å and the Ne vii $`2\mathrm{s}2\mathrm{p}^3\mathrm{P}2\mathrm{s}3\mathrm{d}^3\mathrm{D}`$ multiplet at 106.1 Å.

Fig. 5 illustrates the main features of the generalized collisional–radiative recombination coefficients for O$`{}_{}{}^{4+}+\mathrm{e}^{}`$O<sup>3+</sup>. In the GCR term-resolved picture for light elements, there are four coefficients associated with the pairings of the $`2\mathrm{s}^2{}_{}{}^{1}\mathrm{S}`$ & $`2\mathrm{s}2\mathrm{p}^1\mathrm{S}`$ and $`2\mathrm{s}^22\mathrm{p}^2\mathrm{P}`$ & $`2\mathrm{s}2\mathrm{p}^2{}_{}{}^{4}\mathrm{P}`$ metastable terms in the recombining and recombined systems, respectively. As the radiative and three-body processes are included, the low temperature and high density behaviours reflect these contributions. The finite-density suppression of the coefficient for the ground parent case and the effect of alterative Auger channels are both pronounced and also depend on the ion charge. These effects require GCR modelling. The simpler stage-to-stage picture introduces a significant and, generally, unquantifiable error. It is to be noted that the intermediate-coupled dielectronic recombination data of this project also sustains production of fine-structure-resolved metastable GCR coefficients appropriate to medium- and heavy-weight elements.

The present paper’s main concern is with ensuring the quality and completeness of dielectronic data for plasma modelling and not with all the consequential modelling of populations and ion distributions in plasmas. There are, however, two points to draw attention to. Firstly, it is well known that dielectronic recombination shifts equilibrium ionization balance fraction curves to higher temperatures. This is most pronounced for the ions with one or two electrons outside of closed shells and produces a characteristic ‘piling-up’ of these stages. It is also these ion fractions which show most markedly the effect of finite density reduction of the collisional–radiative coefficients. This cannot be ignored for moderately ionized systems in plasmas with $`N_\mathrm{e}`$10<sup>10</sup> cm<sup>-3</sup>. These effects are shown in Fig. 6 for oxygen. Secondly, most plasma transport models work only with whole ionization stage populations. The present work, however, sustains the metastable-resolved GCR picture. GCR coefficients and ionization balance fractional abundances must be bundled back to the ionization stage for such models, at the expense of precision. Fig. 6b illustrates the resolved picture for the beryllium-like ionization stage of oxygen. The simplest bundling strategy imposes equilibrium fractions on the metastable populations relative to the ground, as is used for the stage-to-stage fractional abundances in Fig. 6a.

Modern good practice requires an estimate of uncertainty in derived theoretical data so that meaningful deductions can be drawn from the comparison with observations. It is unfortunately the case that most theoretical dielectronic data has no error associated with it. Because of the relative complexity of dielectronic recombination and the many contributions, agreements between different theories and with observations sometimes appear fortuitously and do not reflect the underlying reliability (see Savin et al. 2002a). For the present derived GCR coefficients and $`{}_{}{}^{\mathrm{R}}𝒫𝒞`$, we outline our approach to procuring a relevant ‘working error’.

In the ADAS project (Summers et al. 2002), a distinction is made between ‘locked’ parameters, as distinct from ‘search’ parameters, in the optimized fitting of models to observations. Search parameters return a fit uncertainty or confidence level, the locked parameters must carry an error with them. An effective rate coefficient is such a locked parameter. Its uncertainty, called the cumulative statistical error, is computed from the errors of the fundamental reaction rates as follows: Monte Carlo samples are made of all the individual reaction coefficients, within their (assumed) independent Gaussian uncertainty distributions, and the derived coefficient calculated. The process is repeated many times until statistics are built up. The accumulated results are fitted with a Gaussian variance.

The key issue then is the starting point of uncertainties in the fundamental component dielectronic coefficients. The BBGP codes described in section 2.3 have been arranged to generate adf09 baseline files. Such a file is differenced with the matching level 1 file to provide an error estimate for the baseline values and may be stored in a .err file exactly paralleling the naming of the actual .dat file. In like manner, the level 1 file may be differenced with the level 2 file (averaged-over fine-structure) to provide the .err file for level 1. We treat this also as a conservative error for level 2. It is emphasized that this is not a confident absolute error, but a (hopefully) helpful appraisal of the theoretical data. It is most appropriate for the $`n`$\- and $`nl`$\- shell bundled data. The experimental comparisons of the type discussed in section 3 indicate that a minimum uncertainty $``$20% is appropriate for the term and level selective dielectronic data. Use of such .err files is not yet a common practice and its handling within a projection matrix framework is complex. $`{}_{}{}^{\mathrm{R}}𝒫𝒞^z`$ error surfaces are shown in Fig. 7 using the ADAS procedure, but propagating error only from the state-selective part. The full handling of error will be treated in a separate paper.

## 5 Structure and access to the database

The complete set of dielectronic recombination data (both ‘partial’ and ‘total’) will be publically available as adf09 files from the Controlled Fusion Atomic Data Center at the Oak Ridge National Laboratory, USA (http://www-cfadc.phys.ornl.gov/). These data files are simple ascii text in a formatted organization. The layout differs slightly between the level 1 LS-coupled and level 2 intermediate-coupled forms. Fig. 8 summarizes the intermediate-coupled form. The document appxa\_09.pdf of Appendix A of the ADAS user manual (available at: http://adas.phys.strath.ac.uk) provides the detailed description. This includes a summary of the sublibraries and their current status, content of the data lines and the meanings of all parameters, together with some samples of the format. An example of a pathway to a member is /../adf09/jc00#li/jc00#li\_ne7ic23.dat which distinguishes the producer initials ‘$`jc`$’, year number ‘00’, recombining ion iso-electronic sequence ‘$`li`$’, element ‘$`ne`$’, coupling ‘$`ic`$’ and parent $`n=2n=3`$ transition group ‘23’.

We have found it convenient to archive also the driver data sets, which initiate the level 1 and level 2 calculations, as other ADAS data formats (adf27 and adf28). These have pathways which parallel adf09. Various ADAS codes execute the primary calculations and the subsequent collisional–radiative modelling. The flow of calculation is summarized in Fig. 9.

## 6 Summary

We have described the goals and methodology behind a programme to calculate a comprehensive database of dielectronic recombination data for the collisional–radiative modelling of dynamic finite-density plasmas and illustrated its use in such environments. The first phase of the program covering H- through Ne-like sequences is under way and illustrative results, comparisons, and total (zero-density) rates will be the subject of a series of papers to be submitted to A&A in the near future, e.g. O-like ions, Zatsarinny et al. (2003).