text:
# Interference effects in the photorecombination of argonlike Sc3+ ions: Storage-ring experiment and theory

## I Introduction

Dielectronic recombination (DR)

$$e^{}+A^{q+}[A^{(q1)+}]^{}A^{(q1)+}+h\nu $$
(1)
is a two-step electron-ion collision process, where in a first step a multiply-excited intermediate state is created by a resonant dielectronic capture (DC, inverse Auger) and in a second step that intermediate state decays by photon emission. Another possible recombination mechanism is radiative recombination (RR)

$$e^{}+A^{q+}A^{(q1)+}+h\nu $$
(2)
which by emission of a photon directly proceeds to a non-autoionizing state of the recombined ion. When both the initial and the final state including the photons emitted are the same, RR and DR are indistinguishable and they may interfere, with the signature of interference being an asymmetric DR resonance line profile in the energy differential recombination cross section Alber1984 ; Badnell1992 ; Zimmermann1997 ; Behar2000b . Asymmetric line shapes have experimentally so far been observed only in partial recombination cross sections of highly charged uranium ions Knapp1995a . Such an observation is interesting from a fundamental point of view since it questions the widely used independent processes approximation (IPA) in the calculation of RR+DR recombination cross sections.

One motivation for investigating the recombination of argonlike Sc3+(3s23p6) ions with free electrons is the theoretical predication by Gorczyca *et al.* Gorczyca1997  of a strongly asymmetric line shape for the 3p53d$`^{22}`$F DR resonance caused by interference between the recombination pathways

$`\mathrm{e}^{}+\mathrm{Sc}^{3+}(3\mathrm{s}^23\mathrm{p}^6)`$ $``$ $`\mathrm{Sc}^{2+}(3\mathrm{s}^23\mathrm{p}^53\mathrm{d}^2)`$
$``$ $``$ (3)
$`\mathrm{Sc}^{2+}(3\mathrm{s}^23\mathrm{p}^63\mathrm{d})+\mathrm{h}\nu .`$
Theoretically they predicted a Sc2+(3s23p53d$`^{22}`$F) resonance position of 14.6 eV and a width of 1.8 eV. The large width of the resonance is explained by the fact that it decays predominantly via a very fast Super-Coster-Kronig transition (reverse of the horizontal arrow in Eq. 3) with a calculated rate of $`2.84\times 10^{15}`$ s-1.

A first attempt to experimentally determine the line shape of the Sc2+(3p53d$`^{22}`$F) resonance by measuring the total Sc3+ recombination rate coefficient at a heavy-ion storage ring failed due to the rather large statistical uncertainty of that measurement Schippers1999a  (see also Sec. II). However, distinct differences between theoretical and experimental peak positions were found and attributed to the fact that a correct theoretical description of the highly correlated 3p53d2 configurations is extremely difficult. Using Cowan’s Cowan1981  atomic structure code we find that the Sc2+(3p53d2) configuration gives rise to 45 fine structure components distributed over an energy interval of about 20 eV width. In order to account for correlation effects Hansen and Quinet Hansen1996  considered the interaction of 10 initial and 16 final configurations in the calculation of 3p63d $``$ 3p53d2 transition energies in Ca+. Still their result for the Ca+(3p63d$`^2`$D $``$ 3p53d$`^{22}`$F) transition energy deviates 0.76 eV from the experimental result (29.34 eV) determined by photoionization of metastable Ca+(3p63d) ions Mueller2001 . This again illustrates the difficulty to accurately describe inner shell transitions in atomic or ionic systems with open 3d-subshells. In this situation DR measurements (and complementary photoionization measurements) can yield valuable spectroscopic information on such systems.

Here we present new experimental as well as theoretical results for the total Sc3+ recombination rate coefficient. The experiment is a twofold extension of our previous work Schippers1999a . The new data have much smaller statistical uncertainties and they are extended to the much wider energy range 0–45 eV (previous range 12.2–18.2 eV) including all DR series limits due to 3p $``$ 3d and 3p $``$ 4s core excitations.

We have also carried out both nonperturbative and perturbative calculations for the $`e^{}`$ \+ Sc3+ photorecombination. The nonperturbative method is based on the rigorous continuum-bound transition theory and the close-coupling R-matrix approach Zhao2000 . The perturbative evaluation is a treatment including radiative recombination (RR), dielectronic recombination (DR) and their interference. Our calculations reveal that three adjacent resonances mask the low energy side of the 3p53d22F resonance. This makes the experimental observation of the predicted Gorczyca1997  interference between the 3p53d2(3F) 2F DR resonance and the continuous RR cross section difficult. Moreover, we find that several experimental DR resonances not reproduced by the LS-coupling R-matrix calculations of Gorczyca *et al.* Gorczyca1997  are due to relativistic effects.

In Sec. II and Sec. III we outline our experimental and theoretical methods. Experimental and theoretical results are presented and discussed in Secs. IV.1 and IV.2, respectively. In Sec. IV.3 possible origins for the asymmetry of the 3p53d2(3F) 2F DR resonance are explored by theoretically considering interference between DR and RR as well as interference between overlapping DR resonances. Our main results are summarized in Sec. V.

## II Experiment

The experiment was performed at the heavy ion storage ring TSR of the Max Planck Institut für Kernphysik (MPIK) in Heidelberg. Details of the experimental setup and the data reduction procedures can be found in Refs. Kilgus1992 ; Lampert1996 ; Schippers2000b ; Schippers2001c . The 45Sc3+ ions were accelerated by the MPIK tandem booster facility to their final energy and subsequently injected into the TSR. The rather low charge to mass ratio of $`q/A=1/15`$ leads to the rather unfavorable condition (see below) that the ion energy is limited by the maximum bending power of the TSR dipole magnets. With the maximum magnetic rigidity of $`B\rho =1.4`$ Tm the highest accessible ion energy can be calculated from

$$E_\mathrm{i}/A931.5\text{MeV/u}\left[\sqrt{1+\left(\frac{q}{A}\frac{B\rho }{3.107\text{Tm}}\right)^2}1\right]$$
(4)
to $`E_\mathrm{i}/A420`$ keV/u. In the storage ring the circulating 18.9 MeV Sc3+ ion beam was merged with the magnetically guided electron beam of the electron cooler. During electron-cooling, the electrons have to move with the same average velocity as the ions. This condition defines the cooling energy $`E_\mathrm{c}=(E_\mathrm{i}/A)(m_\mathrm{e}/\mathrm{u})230`$ eV. In order to increase the phase space density of the already stored ions and thereby free phase space for the next injection pulse from the accelerator, electron-cooling was allowed to be effective for 2 seconds after each injection of ions. These ‘ecool-stacking’ cycles Grieser1991  were repeated three times until the accumulated ion current stabilized at values of up to 5 $`\mu `$A. Before data taking was started a prolonged cooling interval of 5 s after the last injection allowed the ion beam to shrink to a final diameter of about 2 mm as verified with a beam profile monitor based on residual gas ionization Hochadel1994a .

During the measurement the electron energy $`E_\mathrm{e}`$ in the laboratory frame was stepped through a preset range of values different from $`E_\mathrm{c}`$ thus introducing non-zero relative energies

$$E_{\mathrm{rel}}=(\sqrt{E_\mathrm{e}}\sqrt{E_\mathrm{c}})^2$$
(5)
between ions and electrons. In the data analysis a relativistically correct expression Schippers2000b  was used instead of Eq. 5. Recombined Sc2+ ions were counted with $`100\genfrac{}{}{0pt}{}{+0}{3}`$% efficiency as a function of cooler voltage on a single particle detector Rinn1982  located behind the first dipole magnet downstream of the electron cooler. The dipole magnet bends the circulating Sc3+ ion beam onto a closed orbit and separates the recombined Sc2+ ions from that orbit. In between two measurement steps the cooler voltage was first set back to the cooling value in order to maintain the ion beam quality and then set to a reference value at $`E_{\mathrm{rel}}=45`$ eV which is chosen to lead to a relative velocity where the electron-ion recombination coefficient is only due to a negligible RR rate. Under this condition the recombination rate measured at the reference point monitors the background signal. Choosing short time intervals of 10 ms duration for dwelling on the measurement, cooling and reference voltages ensured that the experimental environment did not change significantly in between signal and background measurement. An additional interval of 1.5 ms after each change of the cooler voltage allowed the power supplies to reach the preset values before data taking was started.

At the rather low ion energy of 420 keV/u the cross section for electron capture from residual gas molecules is estimated — taking the measured residual gas composition into account and using a semi-empirical formula Schlachter1983  for the charge capture cross section $`\sigma ^{(\mathrm{capt})}`$ — to be as large as $`\sigma ^{(\mathrm{capt})}=4\times 10^{18}`$ cm2. Electron capture from residual gas molecules is therefore expected to contribute significantly to the measured recombination signal as a background even at the TSR residual gas pressure of only $`5\times 10^{11}`$ mbar. This already proved to be all the more the case in a previous TSR experiment with Sc3+ ions Schippers1999a  where due to a technical defect in one of the TSR dipole power supplies a reduced maximum rigidity of only $`B\rho 1.2`$ Tm had been available. Consequently, the ion energy had been limited to 300 keV/u. Because of the rapid increase of $`\sigma ^{(\mathrm{capt})}`$ with decreasing ion energy and because the experiment had been carried out in early summer, i. e., close to a maximum in the slight seasonal variation of the TSR residual gas pressure, the measured signal to background ratio had been of the order of 1/100. Under these conditions the Sc3+ recombination rate coefficient could only be measured over a limited energy range with large statistical uncertainties. In the present measurements, that were carried out in December 2000, the signal to background ratio was a factor of about 10 higher as compared to the previous experiment.

The electron-ion recombination coefficient $`\alpha (E_{\mathrm{rel}})`$ is obtained as the background-subtracted recombination count rate normalized to the electron current and to the number of stored ions Kilgus1992 . The estimated systematical uncertainties are 15% for the absolute value of the measured rate coefficient and less than 2% for the relative energy in the energy range under study.

## III Theory

### III.1 Nonperturbative approach of continuum-bound transitions

Starting from the rigorous continuum-bound transition theory, Davies and Seaton Davies1969  discussed the process of emission of radiation into the optical continuum due to radiative capture of an electron by an ion and derived a general formalism including radiation damping. The details of their formalism and our numerical scheme for application can be found in Ref. Davies1969  and Ref. Zhao2000 , respectively. Nevertheless, a brief review of the formalism and our method follows.

Using the close-coupling R-matrix approach (see Zhao2000  for details) we start from a treatment of the electron-ion collision process that neglects any interaction with the radiation field. The resulting wave functions are used as a basis for setting up equations including the interaction with the radiation field in the approximation that only the electric-dipole terms are retained and that the radiation field is restricted to one-photon and no-photon states. The exact solutions for the probability amplitudes of the time-dependent matrix equations involved can be obtained by the application of a Laplace transform. They are expressed in terms of a scattering matrix $`\chi `$ with the partitioning

$$\chi =\left(\begin{array}{cc}\chi _{\mathrm{ee}}& \chi _{\mathrm{ep}}\\ \chi _{\mathrm{pe}}& \chi _{\mathrm{pp}}\end{array}\right)$$
(6)
with $`\chi _{\mathrm{ee}}`$ representing the submatrix for electron-electron scattering allowing for radiative decays, $`\chi _{\mathrm{ep}}`$ that for photoionization, $`\chi _{\mathrm{pe}}`$ that for electron capture with the emission of a photon, and $`\chi _{\mathrm{pp}}`$ that for photon-photon scattering. $`\chi _{\mathrm{ee}}`$ and $`\chi _{\mathrm{pe}}`$ are written as

$`\chi _{\mathrm{ee}}`$ $`=`$ $`S[12\pi ^2D(1+Z)^1D^{}],`$ (7)
$`\chi _{\mathrm{pe}}`$ $`=`$ $`2\pi i(1+Z)^1D^{},`$ (8)
respectively, where $`S`$ is the usual electron-electron scattering matrix in the absence of interaction with radiation fields, and $`D`$ is the reduced dipole matrix ($`D^{}`$ denotes its hermitian conjugate) with its matrix element in the form

$$D_{\gamma J,\gamma ^{}J^{}}=\left(\frac{2\omega ^3\alpha ^3}{3\pi }\right)^{1/2}\frac{\gamma JR\gamma ^{}J^{}}{(2J+1)^{1/2}}$$
(9)
where $`\alpha `$ is the fine-structure constant, $`\omega `$ is the photon energy in units of hartrees, and $`R=_i𝐫_i`$ is the dipole operator with the summation extending over all atomic electrons. The quantum numbers $`\gamma J`$ and $`\gamma ^{}J^{}`$, respectively, specify the continuum and bound states of the atomic system, and $`J(J^{})`$ are the total angular momenta. The wave function of the continuum electron is normalized per unit hartree. In Eqs. 7 and 8 the matrix

$$Z(\mathrm{\Omega })=i\pi 𝑑E\frac{D^{}(E)D(E)}{(E\mathrm{\Omega })}$$
(10)
is related to radiation damping. The variable $`\mathrm{\Omega }`$ denotes the total energy of the ion+photon system. In the usual first-order theory damping is neglected and $`\chi _{\mathrm{pe}}=2\pi iD`$. For the calculation of the PR cross section with damping we employ the numerical method developed by Zhao et al. Zhao2000  for the evaluation of the principal value of the integral appearing in Eq. 10.

In the present calculation the parameters for the 1s, 2s, 2p, 3s, 3p orbitals of scandium were taken from the compilation of Clementi and Roetti Clementi1974 , and the 3d orbital was optimized on the 3p53d 1,3Po, 1,3Do, 1,3Fo states (weighted) by using the CIV3 code of Hibbert Hibbert1975 . We also evaluate PR cross sections for further target states including 4s, 4p, 4d orbitals. All these orbitals were optimized in a way similar to the 3d orbital. The symmetries $`{}_{}{}^{2S+1}L_{}^{\pi }`$ with $`L5`$ were included in our calculations. Here, S and L are the total spin and orbital angular momentum quantum numbers, respectively, and $`\pi `$ denotes the parity. In our calculation we included configurations of odd parity only, since these are expected to lead to the dominant resonances in the energy region of our primary interest around the 3p53d2(3F) 2F resonance.

### III.2 Perturbation theory

On the basis of the principle of detailed balance, the photorecombination (PR) cross sections (in atomic units) from an initial continuum state $`j`$ to a final bound state $`f`$ may be written in the form Sobelman

$$\sigma _{jf}^{\mathrm{PR}}=\frac{g_f}{g_j}\frac{\alpha ^2\omega ^2}{2ϵ}\sigma _{fj}^{\mathrm{PI}}$$
(11)
where $`g_j`$ and $`g_f`$ are the statistical weights of the initial ion core and the final recombined ion, $`\omega `$ and $`ϵ`$ are the photon energy and the free electron energy in hartrees, and $`\sigma _{fj}^{\mathrm{PI}}`$ denotes the photoionization (PI) cross section from state $`f`$ to state $`j`$. In first order perturbation treatment it is

$$\sigma _{fj}^{\mathrm{PI}}=\frac{4\pi ^2\alpha \omega }{3}|M_{fj}|^2$$
(12)
where $`M_{fj}`$ is the PI matrix element of the corresponding transition with the continuum state being normalized per unit energy (hartree). According to continuum-bound configuration-interaction theory Fano1961 , the perturbative PI matrix element $`M_{fj}`$ in low orders may be written in the form (see also Gorczyca1997 ; Mitnik1999 ),

$`M_{fj}`$ $`=`$ $`j|R|f\left(1i{\displaystyle \underset{d}{}}{\displaystyle \frac{\mathrm{\Gamma }_{dj}^\mathrm{a}/2}{\mathrm{\Delta }_d+i\mathrm{\Gamma }_d/2}}\right)`$ (13)
$`+{\displaystyle \underset{d}{}}{\displaystyle \frac{j|V|dd|R|f}{\mathrm{\Delta }_d+i\mathrm{\Gamma }_d/2}}`$
where $`\mathrm{\Gamma }_d`$ is the summation of radiative and autoionization widths of the resonance state $`d`$, $`\mathrm{\Gamma }_{dj}^\mathrm{a}`$ is the autoionization width from state $`d`$ to state $`j`$, $`\mathrm{\Delta }_d=ϵϵ_d`$, in which $`ϵ_d`$ is the energy level of resonance $`d`$, $`V`$ is the electron-electron interaction. It should be noted that in the usual independent-processes and isolated-resonance approximations all cross terms in $`|M_{fj}|^2`$ are omitted that lead to interference between resonant and non-resonant recombination channels as well as to the interaction between resonances. In the present investigation these cross terms are taken into account where necessary.

## IV Results and Discussion

### IV.1 Experimental results

Fig. 1 shows the measured Sc3+(3s23p$`^{61}`$S) recombination rate coefficient. The most prominent feature is the peaklike structure around $`E_{\mathrm{rel}}=40`$ eV. It can be attributed to mainly unresolved high $`n`$ Rydberg DR resonances attached to 3p53d$`^1`$P and 3p54s$`^{1,3}`$P core excitations. Vertical lines indicate DR resonance positions $`E(n)`$ as calculated with the Rydberg formula $`E(n)=E(\mathrm{})q^2/n^2`$ with $`=13.606`$ eV, $`q=3`$, and the series limits $`E(\mathrm{})`$ taken from Tab. 1. It should be noted that we do not observe Rydberg resonances to arbitrary high $`n`$, since recombined Sc2+(3p6$`nl`$) ions in weakly bound high $`n`$ states are field-ionized by motional electric fields in the charge-analyzing dipole magnet and therefore do not reach the recombination detector. Applying a detailed model of the field-ionizing properties of our apparatus that takes into account the radiative decay of high Rydberg states on their way from the electron cooler to the charge-analyzing magnet as well as state-selective field-ionization probabilities Schippers2001c  we estimate the highest detected Rydberg state to be approximately $`n=20`$.

In principle, Rydberg series attached to further 3p53d core excitations (Tab. 1) could have been expected to be visible. The respective core-excitation energies are marked by vertical arrows in Fig. 1. However, no strong DR resonances belonging to such series are observed. Apparently, only dipole allowed (under $`LS`$-coupling conditions) 1S $``$1P core excitations give rise to strong DR channels. This has already been observed for iso-electronic Ti4+ ions Schippers1998  where DR proceeds almost exclusively via resonances attached to 3p53d$`^1`$P excitations. In contrast to the present finding for Sc3+, in the DR spectrum of Ti4+ no DR resonances due to 3p$``$4s excitations were found. One strong hint to the origin of this difference is the fact that in Sc3+ the 3p $``$ 4s excitation energies are about 1 eV lower than the 3p53d$`^1`$P excitation energy (Tab. 1), while in Ti4+ they are higher by up to 6 eV. Although being iso-electronic Sc3+ and Ti4+ do have a markedly different structure.

The recombination rate coefficient rises sharply towards zero relative energy (lower inset of Fig. 1). This rise is due to RR. The RR rate coefficient is rather small as can be expected from the fact that RR into the K, L, M1 and M2,3 shells of argonlike Sc3+ is not possible. This is in contrast to e. g. bare Li3+ ions where these channels yield the major contribution (43% at $`E_{\mathrm{rel}}=10^6`$ eV) to the total hydrogenic RR cross section Burgess1964b  summed up to $`n=20`$. The hydrogenic total Sc3+ rate coefficient, i. e. the sum of the $`nl`$-differential Li3+ RR cross sections ranging from the 3d-subshell to the experimental cutoff quantum number $`n=20`$ convoluted with the experimental electron-energy distribution (see e. g. Ref. Schippers2001c ), is plotted as the full line in the lower inset of Fig. 1. It agrees with the experimental curve within the experimental uncertainties. In spite of the long measuring time of 60 hours, in particular in the energy range of 8–17 eV, the statistical uncertainty did not become low enough to bring out all structures. At the present level of statistical uncertainty it cannot be decided whether the faint structures in the energy range 1–7.5 eV are due to weak DR resonances or due to unaccounted background modulations.

The DR resonances visible in the upper inset of Fig. 1 are due to 3p53d2 and 3p53d 4$`l`$ doubly excited configurations. These configurations are considerably shifted down in energy as compared to the Rydberg energies shown by vertical lines in Fig. 1 and the energies of the various terms must be taken from detailed calculations. As already mentioned in the introduction the terms belonging to the 3p53d2 configuration straddle an energy range of about 20 eV. Also the 3p53d 4$`l`$ configurations can be expected to be subject to strong fine structure interactions and hence to be spread out over several electron volts. The energy range between about 8 and 26 eV (upper inset in Fig. 1) was searched preferentially in order to detect the broad 3p53d22F Super-Coster-Kronig resonance predicted by theory. The measurement reproduces the strong double peak structure near 15 eV found already in the previous experiment Schippers1999a  (scan range 12–18 eV) but now also reveals considerable additional structure both above and below these peaks. While much of this structure seems to arise from a superposition of several narrow peaks, a single broad resonance (with some narrower structure superimposed on its low energy slope) seems to dominate the recombination rate near 12.3 eV. This signal bears strong similarity with the predicted 3p53d22F Super-Coster-Kronig resonance Gorczyca1997  and we take this feature as a key to assigning a considerable part of the observed spectrum. This interpretation is also supported by the argument that for the isoelectronic system Ti4+ Schippers1998  the 3p53d22F resonance was observed about 3 eV below the prediction obtained by the same theoretical approach Gorczyca1998a , so here a position of $`12`$ eV can be expected.

A closeup of the experimental spectrum is shown in Fig. 2, where the resonances in the energy range 10.9–16.5 eV have been fitted by either Lorentzian line profiles

$$L(E)=\frac{A}{\mathrm{\Gamma }_d}\frac{2}{\pi }\frac{1}{1+\epsilon ^2}$$
(14)
or an asymmetric Fano profile

$$F(E)=\frac{A}{Q^2\mathrm{\Gamma }_d}\frac{2}{\pi }\left[\frac{(Q+\epsilon )^2}{1+\epsilon ^2}1\right]$$
(15)
with $`\epsilon =2(EE_d)/\mathrm{\Gamma }_d`$. The term -1 in the square brackets ensures that $`F(E)0`$ for $`E\pm \mathrm{}`$ and the normalization factors have been chosen such that $`F(E)L(E)`$ in the limit $`Q\mathrm{}`$. Conversely, small absolute values of the asymmetry parameter $`Q`$ lead to significantly asymmetric line shapes. Assuming that the dominant contribution to the resonance width $`\mathrm{\Gamma }_d`$ stems from autoionization and neglecting contributions from neighboring resonances, $`Q`$ can be expressed in terms of the matrix elements occurring in Eq. 13, i. e.,

$$Q=\frac{2j|V|dd|R|f}{\mathrm{\Gamma }_dj|R|f}$$
(16)
as can be shown by simple algebraic manipulations on $`|M_{fj}|^2`$. For each resonance the parameters varied in the fit were the peak area $`A`$, the resonance energy $`E_d`$, the resonance width $`\mathrm{\Gamma }_d`$, and in case of the 3p53d$`^{22}`$F resonance the asymmetry parameter $`Q`$. In the fit both $`L(E)`$ and $`F(E)`$ have been convoluted (see appendix) by a Gaussian representing the experimental energy-distribution function with its full width at half maximum (FWHM) $`\mathrm{\Delta }E4(E_dk_\mathrm{B}T_{}\mathrm{ln}2)^{1/2}`$ in the energy region of interest Kilgus1992 . The longitudinal electron-beam temperature $`k_\mathrm{B}T_{}`$ ($`k_\mathrm{B}`$ denotes the Boltzmann constant) has also been allowed to vary during the fit. From the fit we obtain $`k_\mathrm{B}T_{}=0.16\pm 0.03`$ meV. Further fit results are summarized in Tab. 2.

The most important result of the fit is the determination of the 3p53d$`^{22}`$F resonance energy and width to $`12.31\pm 0.03`$ eV ($`\pm 2\%`$ systematic uncertainty) and $`0.89\pm 0.07`$ eV, respectively. For the asymmetry parameter of this resonance we obtain $`Q=6.3\pm 1.8`$. A major source of uncertainty beyond the errors given is the background from collisions with residual gas molecules. Under the conditions of the present experiment the background causes the major fraction of the measured recombination count rate (see Sec. II). A subtraction of this high background level from the measured signal is difficult, especially, if changes of the residual gas pressure have occurred on the time scale of the switching between signal and background measurement. In this case a small fraction of the background may not have been properly subtracted possibly leading, e. g., to the nonzero recombination signal in the energy range 15.6–17 eV. In principle, this signal might also be caused by a group of weak DR-resonances, however, this assumption is not supported by our theoretical calculations (see below). When allowing a constant background level as an additional free parameter in the fit a considerably different asymmetry parameter $`Q=10.4\pm 4.1`$ is obtained with the fitted background level amounting to $`6.0\pm 0.9\times 10^{13}`$ cm3s-1. On the other hand the resonance position and width change only within their errors as listed in Tab. 2 to $`12.35\pm 0.03`$ eV and $`0.85\pm 0.07`$, respectively.

The significance of our fit result for the asymmetry parameter may also be assessed by fitting the 3p53d$`^{22}`$F resonance with a *symmetric* Lorentzian instead of an *asymmetric* Fano line shape. From such a Lorentzian fit we obtain a $`\chi ^2`$ value that is higher by only 5% as compared to the Fano fit, i. e., the experimental evidence for an asymmetric line shape is not very strong.

Fig. 3 shows the comparison of the present experimental result with the theoretical result of Gorczyca *et al.* Gorczyca1997  in the region of the 3p53d$`^{22}`$F resonance. As already inferred from our previous work Schippers1999a  the theoretical resonance positions do not match the experimental ones. For the comparison presented in Fig. 3 the theoretical energy scale has been shifted by 2.2 eV in order to match the experimental position of the broad 3p53d$`^{22}`$F resonance. Moreover, the absolute theoretical rate coefficient has been normalized to the experimental height of this resonance by a multiplication factor of 0.63. After these adjustments a basic similarity between predicted and measured spectrum is found. In particular, a broad feature is observed strongly resembling the calculated 3p53d22F resonance. However, in the experimental spectrum considerably more resonances are visible than in the theoretical one. The theoretically predicted single peaks at 11.6 and 15.1 eV appear to be split into two components. Also, the predicted relative resonance strengths deviate from the experimental findings. The theoretical 3p53d$`^{22}`$P peak at 15.1 eV is much smaller than the experimental doublet.

### IV.2 Theoretical results

The result of our LS-coupling–R-Matrix calculation of the total Sc3+ PR cross section including damping is plotted in Fig. 4a. It can be noticed that there is a very wide resonance at about 13.6 eV. As already predicted by Gorczyca *et al.* Gorczyca1997 , it is due to DR via the 3p53d2(3F) 2F state. However, our resonance energy is by about 1 eV closer to the experimental value of 12.29 eV (see Tab. 2) than that of Gorczyca *et al*. This difference originates from the use of different sets of basis functions in both calculations. Table 3 lists measured and theoretical resonance energies including further theoretical results Tiwary1983 ; Altun1999 . The relatively large differences between the various theoretical results illustrates again the difficulty to exactly describe the many-electron atomic system under consideration.

In order to be able to more accurately compare the theoretical result with the measured merged-beams rate coefficient we convoluted the R-matrix cross section with the experimental electron energy distribution function, that can be parameterized by the longitudinal and transverse (with respect to the electron beam direction) electron beam temperatures $`T_{}`$ and $`T_{}`$, respectively Kilgus1992 . For the convolution result displayed in Fig. 4b we have used $`k_\mathrm{B}T_{}=0.15`$ meV (cf., Fig. 2) and $`k_\mathrm{B}T_{}=10`$ meV. Comparing with the experimental data in Fig. 1 we find that the R-matrix calculation does not reproduce the experiment in every detail. For example, the calculated resonance peaks in the energy range 7–17 eV are too high, whereas those beyond 19 eV are too low. Moreover, some of the observed peaks at lower energies are not reproduced by the R-matrix calculation. The disagreement at high energies can be traced back to the neglect of even-parity configurations and to the neglect of relativistic effects by restricting the R-Matrix calculation to $`LS`$-coupling. This also explains the missing of peaks at lower energies. Using Cowan’s code Cowan1981  for the calculation of DR resonances in the framework of perturbation theory (Sec. III.2) we find that due to relativistic effects (mainly spin-orbit interaction) also $`LS`$-forbidden excitations to 3p53d2(3F) 2D as well as to 3p53d(1D)4s 2D and 3p53d(3D)4s 2D doubly excited states may contribute to the PR cross section. It should be emphasized that resonances arising from relativistic effects have in general small Auger widths if the atomic number $`Z`$ is not large. As has been explained in some detail by Zhao and Shirai Zhao2001  they, nevertheless, may still play an important role in recombination if their radiative lifetime is of the order of the radiative lifetimes of the dominant DR resonances.

As shown in Fig. 5, a much better agreement with the experimental data is achieved when the the 2D resonances mediated by spin-orbit effects are taken into account in addition to the R-Matrix cross section. For the generation of the incoherent sum of the R-Matrix + perturbation-theory rate coefficient shown in Fig. 5 we have used the 2D Auger and radiative rates as calculated by Cowan’s code. The resonance energy, however, has been determined by using the experimentally observed energy splittings between the 2D resonances and neighboring resonances obtained from R-Matrix theory. For example, the experiment yields 3p53d2(3P) 2P and 3p53d2(3F) 2D resonance energies of 14.8 and 15.3 eV, respectively (Tab. 2). The level difference equals 0.5 eV. Since our R-matrix calculation yields a 3p53d2(3P) 2P resonance energy of 16.4 eV, we chose 16.9 eV to be the 3p53d2(3F) 2D resonance energy. Moreover, for the comparison in Fig. 5 the theoretical merged-beams rate coefficient was multiplied by a factor 0.55 and shifted by 1.13 eV to lower energies. The difference between experimental and theoretical cross sections of nearly a factor 2 is found to be independent of the inclusion of 4s, 4p and 4d orbitals in the set of target states.

### IV.3 Line shape of the 3p53d2(3F) 2F resonance

Clearly, the theoretical line shape of the broad 3p53d2(3F) 2F resonance is asymmetric. As already mentioned in the introduction, Gorczyca *et al.* Gorczyca1997  have attributed this asymmetry to interference between RR and DR pathways. However, this interpretation may be doubted in view of the fact that in the relevant energy range the RR cross section is orders of magnitude smaller than the DR cross section. Naively, one would expect that a significant interference between DR and RR requires cross sections of similar magnitude. It must be kept in mind, however, that the asymmetry parameter $`Q`$ is also governed by the natural line width $`\mathrm{\Gamma }_d`$. Inspection of Eq. 16 shows that a large value of $`\mathrm{\Gamma }_d`$ can make up for a small RR transition matrix element $`j|R|f`$. Fig. 6a shows the RR, DR and interference contributions to the recombination cross section in the energy region of the 3p53d2(3F) 2F resonance as calculated by perturbation theory. With the width of the 3p53d2(3F) 2F resonance due to its Super-Coster-Kronig decay channel being extraordinarily large ($`\mathrm{\Gamma }_d=1.06`$ eV obtained from Cowan’s atomic-structure code Cowan1981 ) the asymmetry in the merged-beams rate coefficient is retained even when the RR cross section is artificially lowered by a factor of 10 (Fig. 6b). The theoretical asymmetry parameter changes from $`Q=3.97`$ to $`Q=12.6`$. Our experimental findings of $`\mathrm{\Gamma }_d=0.89`$ eV and of $`Q=6.3`$ (see Tab. 2) correspond to these theoretical results.

Apart from interference between DR and RR another source of distorted line shapes is the interaction between overlapping DR resonances of the same symmetry Griffin1994 . From our R-Matrix calculation we find two further resonances of 2F symmetry, namely the resonances 3p53d(3F)4s 2F and 3p53d(1F)4s 2F located at 11.00 and 12.71 eV, respectively. As shown in Fig. 7a, especially the latter is overlapping with the broad 3p53d2(3F) 2F resonance. In Figs. 7b and 7c we show how the line shape of the 3p53d2(3F) 2F resonance depends on its (artificially varied) distance to the 3p53d(1F)4s 2F resonance. Clearly the shapes of both resonances are affected by their mutual separation. The experimental observation of this effect, however, is hindered by the presence of the 2D resonances close to the 3p53d(1F)4s 2F resonance. Under this circumstance the experimental resolution does presently not suffice to reveal the theoretically predicted subtle resonance-resonance interference.

## V Conclusions

In the experimental energy range 0-45 eV photorecombination of Sc3+(3s23p6) ions is dominated by the 3p53d 1P $`nl`$ and 3p54s 1P $`nl`$ Rydberg series of DR resonances. Significant differences to the recombination of isoelectronic Ti4+ ions Schippers1998  are found that are due to subtle differences in the electronic structure of Sc3+ and Ti4+. Level energies, widths and resonance strengths were determined for a number of doubly excited Sc2+(3s23p53d2) and Sc2+(3s23p53d 4s) states. Furthermore, we find a small experimental evidence for an asymmetry of the 3p53d2(3F) 2F resonance line shape. A more precise determination of the line shape from the experimental data would require less statistical uncertainties especially in the high energy tail of the resonance. In the present measurement the high background level from electron capture in Sc3+ collisions with residual gas molecules prevented this requirement from being fulfilled. Moreover, the peaks on the low energy side of the 3p53d$`^{22}`$F resonance leave additional room for ambiguity.

Starting from rigorous continuum-bound transition theory the Sc3+ photorecombination cross section has been calculated using both the R-matrix method and a perturbative treatment. We found that the line shape of the 3p53d2(3F) 2F resonance is influenced by both interacting resonances and interference between RR and DR. In spite of the fact that the Sc3+ RR cross section is almost negligibly small, DR-RR interference becomes noticeable due to the large 3p53d2(3F) 2F resonance width of 0.89 eV. Moreover, we have shown that relativistic effects play an important role in the recombination of Sc3+. Three experimentally observed strong 3p53d22D and 3p53d4s 2D resonances can only be reproduced theoretically by invoking relativistic effects. Our R-matrix calculations as well as the calculations of Gorczyca *et al.* Gorczyca1997  that both have been restricted to LS-coupling fail in reproducing these strong experimental peaks.

Our results show the limits of the independent-processes and isolated-resonances approximations usually made in photorecombination calculations. A more stringent test of our theoretical results could in principle be provided by the measurement of high-resolution photoionization cross sections of Sc2+(3p63d$`^2`$D) ions. Results of such an experiment will be reported elsewhere Schippers2002 .

###### Acknowledgements.

Technical support by the Heidelberg accelerator and TSR groups is gratefully acknowledged. The experimental work has been supported by the Deutsche Forschungsgemeinschaft (DFG, grant no. Mu 1068/8). A.C and M.E.B were supported by a NATO collaborative grant (no. PST/CLG 976362) and L.B.Z was supported by the National Natural Science Foundation through grant no. 19974006.

## Appendix A Convolution of a Fano-profile with a Gaussian

The convolution of a Fano profile (Eq. 15) with a normalized Gaussian

$$G(E)=\frac{2}{\mathrm{\Delta }E}\sqrt{\frac{\mathrm{ln}2}{\pi }}\mathrm{exp}\left[\frac{4(\mathrm{ln}2)E^2}{(\mathrm{\Delta }E)^2}\right]$$
(17)
with $`\mathrm{\Delta }E`$ being its full width at half maximum (FWHM), yields

$`C(E)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}F(E^{})G(E^{}E)𝑑E^{}`$
$`=`$ $`A{\displaystyle \frac{2}{Q^2\mathrm{\Gamma }\sqrt{\pi }}}\left[{\displaystyle \frac{1}{\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{[Qy+(tx)]^2}{(tx)^2+y^2}}e^{t^2}𝑑t1\right]`$
$`=`$ $`A{\displaystyle \frac{2\sqrt{\mathrm{ln}2}}{\mathrm{\Delta }E\sqrt{\pi }}}\left[\left(1{\displaystyle \frac{1}{Q^2}}\right)\mathrm{}[w(z)]{\displaystyle \frac{2}{Q}}\mathrm{}[w(z)]\right]`$
where the definitions

$`t`$ $`=`$ $`{\displaystyle \frac{2\sqrt{\mathrm{ln}2}(E^{}E)}{\mathrm{\Delta }E}},`$ (19)
$`x`$ $`=`$ $`{\displaystyle \frac{2\sqrt{\mathrm{ln}2}(E_dE)}{\mathrm{\Delta }E}},`$ (20)
$`y`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }\sqrt{\mathrm{ln}2}}{\mathrm{\Delta }E}},`$ (21)
$`z`$ $`=`$ $`x+iy`$ (22)
have been used and where $`w(z)`$ denotes the complex error function Abramowitz1964 . Its real and imaginary parts can efficiently be calculated with a numerical algorithm by Humlicek Humlicek1979 . It should be noted that for $`Q\mathrm{}`$ Eq. A yields the Voigt profile, i. e. the convolution of a Lorentzian (Eq. 14) with a Gaussian.


format: text

text:
# Dynamic Hydrogen Ionization

## 1 Introduction

The solar chromosphere is the region about 1.5 Mm thick between the temperature minimum (at about 4000 K) and the transition region to coronal temperatures of several million Kelvin. Traditionally, its structure has been determined by semi-empirical fitting of temporally and spatially averaged continua and line intensities (e.g., Vernazza et al., 1981; Maltby et al., 1986; Fontenla et al., 1993). However, the chromosphere is actually a very dynamic region. In this paper we investigate the hydrogen ionization structure of a dynamic chromosphere.

It has long been known that hydrogen ionization and excitation in the solar atmosphere is not in equilibrium for the local temperature and density (Thomas, 1948). Usually, hydrogen ionization is calculated from the condition of “statistical equilibrium”, that is, the equality of the ionization and recombination rates for the local temperature, density and radiation. Statistical equilibrium assumes infinitely fast rates and an instantaneous adjustment to the local thermodynamic and radiation state. However, if the local state changes in time, or if there is a flow through a inhomogeneous region, and if the time scale to reach ionization or excitation equilibrium is longer than the dynamic times scale, then it is necessary to solve the population rate equations,

$$\frac{dn_i}{dt}=\mathrm{gains}\mathrm{losses},$$

for each species $`i`$ (Joselyn et al., 1979). This latter is the case in the solar atmosphere heated by any intermittent process (such as shocks or nanoflares) (Elzner, 1975; Klein et al., 1976; Kneer & Nakagawa, 1976; Poletto, 1979; Kneer, 1980). The ionization and excitation state then depends on its history as well as the instantaneous temperature, density and radiation. However, the importance of slow chemical reaction rates in the quiet solar atmosphere has generally been neglected since these early papers were written. The active solar atmosphere (e.g. flares, prominence formation and coronal mass ejections) has long been recognized to require a time dependent analysis (McClymont & Canfield, 1983; Doschek, 1984; Fisher et al., 1985; Heinzel, 1991; Abbett & Hawley, 1999; Sarro et al., 1999; Ciaravella et al., 2001; Ding et al., 2001; Lanza et al., 2001). The same issue of ionization and other chemical reaction rates being slow compared to the rate of dynamical changes, and thus requiring the solution of the species rate equations rather than statistical equilibrium, arises in many other areas of astrophysics where conditions are changing in time. Examples are: the interstellar medium (Lyu & Bruhweiler, 1996; Joulain et al., 1998), HII regions (Rodriguez-Gaspar & Tenorio-Tagle, 1998; Richling & Yorke, 2000), planetary nebula (Schmidt-Voigt & Koeppen, 1987; Frank & Mellema, 1994; Marten & Szczerba, 1997), novae (Hauschildt et al., 1992; Beck et al., 1995), supernova (Kozma & Fransson, 1998), and the reionization of the intergalactic medium and the Lyα forest (Ikeuchi & Ostriker, 1986; Shapiro & Kang, 1987; Shapiro et al., 1994; Ferrara & Giallongo, 1996; Giroux & Shapiro, 1996; Zhang et al., 1997).

The solar chromosphere is not well represented by a static mean structure. We show that because the chromosphere is dynamic the average hydrogen ionization is significantly different from the ionization state of the mean atmosphere. There are two reasons for this. First, the mean of any quantity that depends non-linearly on the atmospheric properties is not the same as that quantity determined from the mean atmosphere. Second, the finite rates of ionization and recombination reduce the response of the ionization to dynamic variations in the atmospheric state. Hydrogen does not have time to reach its equilibrium ionization. The result is a hydrogen ionization fraction that is higher than obtained for the mean of the dynamic atmosphere. This conclusion is independent of just what that mean state is. Thus analyzing observations on the basis of a static model chromosphere leads to very different conclusions from an analysis based on a dynamic model chromosphere. This significantly alters the interpretation of observations.

In this paper, we first describe (section 2) the numerical simulations used to study the dynamic hydrogen ionization. We then show that the dominant hydrogen ionization process is photoionization from the second level and we investigate the processes that populate the second level (section 3). There is a very rapid equilibration of the 2nd and all higher levels with the continuum, with a slow collisional leakage of electrons from the ground state to the 2nd level, or visa versa, depending on whether hydrogen is ionizing or recombining. Next we calculate the time scale for ionization/recombination and its relation to the eigenvalues of the rate matrix (section 4). The ionization and recombination rates are slow and not in equilibrium, so that the level populations change in time. We demonstrate that there is only one relaxation time scale to approach equilibrium, not separate ionization and recombination times scales. We show (section 5) that because of the slow ionization and recombination rates, the dynamic ionization fluctuates much less than the statistical equilibrium ionization. As a result, line and continuum intensity variations do not mimic the underlying dynamics. Further, the mean ionization fraction tends to represent the maximum statistical equilibrium ionization, and is much higher than is obtained from the statistical equilibrium of the mean atmosphere. We conclude (section 6) with a reiteration of why the chromosphere must be dynamic based on fundamental physical principles and observations and a statement of which results are robust and why.

## 2 Numerical Simulations

To properly model the dynamic solar chromosphere one has to perform radiation-hydrodynamic simulations taking into account the non-local, non-linear rate equations for all important species. Such a self-consistent radiation-hydrodynamic modeling of the solar chromosphere has been performed by Carlsson & Stein (1992, 1994, 1995, 1997a, 1997b) and we summarize the methodology here.

We solve the one-dimensional equations of mass, momentum, energy and charge conservation together with the non-LTE radiative transfer and population rate equations, implicitly on an adaptive mesh. Advection is treated using Van Leer’s (1977) second order upwind scheme to ensure stability and monotonicity in the presence of shocks. An adaptive mesh is used (Dorfi & Drury, 1987) in order to resolve the regions where the atomic level populations are changing rapidly (such as in shocks). The equations are solved simultaneously and implicitly to ensure self-consistency and stability in the presence of radiative energy transfer, stiff population rate equations, and to have the time steps controlled by the rate of change of the variables and not by the small Courant time for the smallest zones. A linearization method is used to solve the radiative transfer (Scharmer & Carlsson, 1985) but with a penta-diagonal approximate lambda operator (Rybicki & Hummer, 1991) instead of the global Scharmer operator. The effects of non-equilibrium ionization, excitation, and radiative energy exchange from several atomic species (H, He and Ca) on fluid motions and the effect of motion on the emitted radiation from these species are calculated. We model hydrogen and singly ionized calcium by 6 level atoms and helium with a 9 level atom. For helium we collapse terms to collective levels and include the $`1s^2`$, $`2s`$ and $`2p`$ terms in the singlet system and the $`2s`$ and $`2p`$ terms in the triplet system of neutral helium, and the $`1s`$, $`2s`$ and $`2p`$ terms of singly ionized helium. In addition we include doubly ionized helium. We include in detail all transitions between these levels. For singly ionized calcium they are the H and K resonance lines, the infrared triplet and the photoionization continua from the five lowest levels. We use 31–101 frequency points in each line and 4–23 frequency points in each continuum; a total of 1424 frequency points. Continua from elements other than H, He and Ca are treated as background continua in LTE, using the Uppsala atmospheres program (Gustafsson, 1973).

The upper boundary is a corona at 106K with a transmitting boundary condition. Incident radiation from the corona is included which causes ionization in the helium continua in the upper chromosphere. The lower boundary condition, at 500 km below $`\tau _{500}=1`$, is also transmitting. Waves are driven through the atmosphere by a piston located at the bottom of the computational domain. The piston velocity is chosen to reproduce a 3750 second sequence of Doppler-shift observations in an Fe I line at $`\lambda `$396.68 nm in the wing of the Ca H-line (Lites, Rutten & Kalkofen 1993) This line is formed about 260 km above $`\tau _{500}=1`$.

The initial atmosphere (Fig. 1) is in radiative equilibrium above the convection zone (for the processes we consider) without line blanketing and extends 500 km into the convection zone, with a time constant divergence of the convective energy flux (on a column mass scale) calculated with the Uppsala code without line blanketing. Note that the transition region in the initial atmosphere occurs at a lower height than in standard models. This is because the lower temperature means a smaller pressure scale height and a less extended atmosphere. The atmospheric dynamics moves it outwards on average. The transition region occurs at a smaller pressure than in the VALIIIC model. This is set by the amount of conductive flux at the upper boundary and is a free parameter in the model (follows from the location of the upper boundary at the fixed temperature of 106K). The mean structure of the dynamic atmosphere (Fig. 1) has a low (5000 K) temperature throughout most of the chromosphere, with a temperature rise in the upper chromosphere produced by absorption of coronal radiation in the helium continua.

The difference between the calculations reported here and the ones in the references cited above is that we now include helium, a corona and transition region, the incident radiation from the corona, extend the calculations deeper and have a transmitting lower boundary condition. This makes it possible to discuss the upper chromosphere, in particular the ionization of hydrogen, in some more detail than before.

The validity of these simulations can be checked by comparing their predictions with observations. There is both agreement with some observations and disagreement with other observations. CO molecular line observations indicate that there is no temperature rise in the low chromosphere, consistent with the simulation. The simulations reproduce many details in observations of the calcium H-line (Carlsson & Stein, 1997a). The Ca grains are due to waves that steepen into shocks around a height of 1 Mm. However, the cores of the simulated H and K lines are darker and the bright points are brighter than is observed. The Lyman continuum intensity variation is also larger in the simulations than what is observed. There is a possible disagreement with SUMER observations that the cores of NI, OI, CI and CII lines, formed in the mid and upper chromosphere, show emission everywhere. Although some of these line cores are formed below the height where coronal radiative heating raises the mean chromospheric temperature, there is some contribution from that region of enhanced temperature. Preliminary calculations (that are non-trivial because effects of slow ionization/recombination have to be included, as is shown in this paper) show that this contribution actually cause the lines to be in emission all the time even though the average line emission is weaker than is observed.

The improved physics in the simulations reported on here compared with the previous work (the inclusion of a corona and the absorption of coronal radiation in the helium continua) has improved the agreement; the line cores in the calcium H and K lines are not as dark as in the previous simulations, the rms variation in the Lyman continuum radiation temperature has changed from 256 K (Carlsson & Stein 1994, 1997b) to 97 K (compared to an observed value of about 40 K) and UV lines from neutral elements now seem to be in emission everywhere in the simulations.

The existing disagreements indicate that there is still physics missing from the simulations. For instance, we do not include line blanketing (especially the numerous iron lines), the CO molecule and the singly ionized magnesium atom (producing the Mg h and k lines). We use a crude approximation to partial redistribution of radiation in the Lyman lines and no partial redistribution at all in the calcium H and K lines. We do not include frequencies above 20 mHz in the driving piston. We do not include the effects of magnetic fields.

Assuming complete redistribution in the calcium lines probably leads to an overestimate of the cooling in these lines by up to a factor of two (Uitenbroek, 2002), partially compensating the omission of the magnesium h and k lines. High frequency waves may lead to heating of the mid chromosphere through dissipation in shocks and the neglect of high frequencies in the driving piston probably leads to a lower mean temperature.

Magnetic fields structure the chromosphere and corona (Aschwanden et al., 2001; Berger et al., 1999; Schrijver et al., 1999; Moses et al., 1994). They influence both energy transport (wave modes and quasi-static driving) and heating (nanoflares, resonant wave absorption). In the photosphere magnetic fields are concentrated into isolated flux tubes or loops, except in very high magnetic flux regions, leaving most of the photosphere nearly field free. The upper convection zone is the site of acoustic wave generation (Stein & Nordlund, 2001; Skartlien et al., 2000; Goldreich et al., 1994), while overshooting convective flows in the photosphere generate gravity waves. MHD tube waves are driven by convective motions acting on the localized magnetic fields (Musielak & Ulmschneider, 2001). The photospheric tubes and loops spread out with increasing height and decreasing gas pressure. When acoustic waves propagate into the the region where the magnetic pressure equals the gas pressure ($`\beta =1`$) and the magnetic field lines are curved, there is significant reflection and mode conversion into magneto-hydrodynamic waves propagating at approximately the Alfven speed (Rosenthal et al., 2002). These latter are less compressive and have a smaller amplitude for a given flux than the acoustic waves, because of their faster propagation speed, and hence are less visible. Where there is magnetic field extra heating occurs — continuum intensities and line emissions are substantially higher in magnetic network regions than in the internetwork. It is likely that the same processes will contribute to the heating in internetwork regions in the mid-chromosphere where the magnetic field spreads to form a ”magnetic canopy”.

In this paper we focus on the processes controlling the ionization of hydrogen in a dynamic atmosphere. We show that when chemical rates are slow, populations do not have time to reach their equilibrium values. Indeed, some of the disagreement with observations may be due to this non-equilibrium. The solar chromosphere will experience some additional heating besides that included in our model. However, as we show later, even in such hotter models the ionization/recombination rates are slow and the results presented here will be qualitatively correct.

## 3 Ionization and Excitation Processes

In the upper chromosphere and lower transition region the density is much too low to ensure LTE and the full rate equations have to be solved to calculate the hydrogen ionization. Furthermore, inspection of the rates involved shows that typical ionization/recombination time scales are much longer than dynamical time scales so that the ionization balance can be expected to be out of statistical equilibrium. It is thus necessary to take into account the advection and time-derivative terms in the equations. In this section we analyze the rate equations to find what processes dominate the hydrogen ionization balance. In the next section we analyze the time scales involved.

Figure 2 shows the hydrogen ionization fraction as a function of column mass ($`m_c`$) in the initial radiative equilibrium atmosphere. The ionization fraction rises from about $`10^5`$ at the classical temperature minimum at lg($`m_c`$)=$`1`$ (height of 0.5 Mm) to about 30% at the base of the transition region 1 Mm further up.

The processes that are important for the ionization balance are not the same in the transition region and in the chromosphere. In the solar chromosphere the hydrogen ionization is dominated by photoionization from the first excited state. The reason is that the much lower particle density of the first excited state ($`3\times 10^4\mathrm{cm}^3`$) than in the ground state ($`10^{11}\mathrm{cm}^3`$) is more than compensated by the higher mean intensity at the Balmer edge at 364 nm compared to that of the Lyman edge at 91 nm. At a radiation temperature of 5000 K the ratio is $`3\times 10^8`$. In addition, the atmosphere is optically thick to Lyman photons and optically thin to Balmer photons in all of the chromosphere. The result is that photoionization by Lyman photons plays an insignificant role in the hydrogen ionization throughout the chromosphere. Figure 3 shows the net rates between all the hydrogen energy levels at lg($`m_c`$)=$`4`$ which corresponds to a height of 1.2 Mm in the initial atmosphere. This picture is almost identical in the whole chromosphere from lg($`m_c`$)=$`1`$ up to the base of the transition region at lg($`m_c`$)=$`5.6`$. Ionization is primarily through a net rate from $`n`$=2 to the continuum (photoionization in the Balmer continuum) balanced by recombination to the higher levels cascading down to the $`n`$=2 level through bound-bound transitions with $`\mathrm{\Delta }n`$=1. The reservoir of electrons is the ground state. There is transfer of electrons between the ground state and first excited state by collisions. The rate of this transfer is very slow, because chromospheric temperatures are of order 1 eV and the energy jump to the first excited level is 10.2 eV.

In the very thin zone (2 km) where the ionization fraction goes from 30% to fully ionized the situation is different. Collisional ionization from the ground state dominates with net photorecombination in all continua and again bound-bound transitions with $`\mathrm{\Delta }n`$=1 dominating the return channel to the ground state (Fig. 4).

Figures 3-4 only show the net rates and do not answer the question what is driving the system and what rates are just adjusting to provide a closed loop. One way to investigate this aspect is to see how the system responds to perturbations. We have therefore solved the equations of statistical equilibrium repeatedly perturbing the rates one-by-one. Increasing the photoionization cross-section in the Balmer continuum increases the ionization while the opposite is true for the bound-free transitions from the higher levels, consistent with the rate picture where the ionization is through photoionization in the Balmer continuum and recombination to the higher levels. Even though the H-$`\alpha `$ ($`n`$=3 $``$ $`n`$=2) transition has the largest net rate in the chromosphere changing this rate has no effect on the ionization balance. The same is true for all other bound-bound transitions except for Ly$`\alpha `$. Increasing the Ly$`\alpha `$ rate increases the ionization slightly just below the transition region where absorption in Ly$`\alpha `$ of photons from the transition region and corona populates level 2.

Up to now we have studied only the static atmosphere. Hydrogen ionization in a dynamic chromosphere is illustrated in Fig. 5. In the left column the temperature is shown as a function of column mass as a dashed line with the values shown on the right hand axes. The same column shows the hydrogen ionization fraction on a logarithmic scale as a solid line with values on the left axes. The right column shows the net rates from the bound hydrogen levels $`n`$=1-5 to the continuum with positive values indicating net ionization, negative values recombination. Time progresses from the top panel downwards.

At a time 1600 seconds from the start of the simulation a shock is starting to form around lg($`m_c`$)= $`3.5`$ (height of 1 Mm). The rates are in statistical equilibrium (the sum of the rates is close to zero) with a balance between Balmer photoionization and recombination to the levels $`n`$=2-5. 20 seconds later (mid panels) the shock has formed. Balmer photoionization is greatly increased in the shock, because the population of the first excited state has greatly increased due to collisional excitation from the ground state . (The ratio of collisional excitation to de-excitation rates depends exponentially on temperature, $`C_{12}/C_{21}=g_2/g_1e^{\mathrm{\Delta }E/kT}`$.) The Ly$`\alpha `$ transition is found to be in detailed balance and so just bounces electrons up and down. The photoionization rate per atom is constant because it depends on the radiation field which is set deeper down in the atmosphere. Balmer photoionization dominates the rates into the continuum (dotted curve), but it is not balanced by recombination (because of slow rates) resulting in a net increase of protons (sum of rates positive and almost equal to the Balmer photoionization). As the shock progresses farther upwards, recombination dominates over photoionization behind the shock. The zone where hydrogen goes from 1% to 40% ionization is about 600 km thick.

At the top of the chromosphere hydrogen is between 30 and 40% ionized. The further ionization up to 100% ionization takes place in the zone where the temperature increases from 10000 K to 25000 K. The ionization rates are here dominated by collisional ionization directly from the ground-state while the Lyman continuum actually provides some net recombination. The thickness of this zone depends on the temperature gradient and is thus very model dependent — in the plane parallel calculations this zone is very thin, on the order of 2 km. Also in this zone we see rates that are out of statistical equilibrium with net ionization when shocks pass followed by net recombination.

## 4 Time Scales

The crucial result of our investigation is that the time scale for changes in the ionization of hydrogen in the chromosphere is long, because the rates for ionization and recombination under conditions typical of the solar chromosphere are small. Consequently the hydrogen populations do not adjust to the local conditions.

The relation of time scales and transition rates can be understood most readily for the case of a 2-level atom. We summarize the well known results here. The equations for the populations of levels 1 and 2 are

$`{\displaystyle \frac{Dn_1}{Dt}}`$ $`=`$ $`n_2P_{21}n_1P_{12}`$ (1)
$`n_1+n_2`$ $`=`$ $`n_H,`$ (2)
where $`n_i`$ is the population of level $`i`$ and $`P_{ij}`$ is the transition rate per atom (radiative $`R_{ij}`$ plus collisional $`C_{ij}`$) from level $`i`$ to level $`j`$, and $`n_H`$ is the total hydrogen number density. The solution for the population of level $`1`$ is

$$n_1(t)=n_1(\mathrm{})+\left(n_1(0)n_1(\mathrm{})\right)e^{t(P_{21}+P_{12})},$$
(3)
where $`n_1(\mathrm{})`$ is the equilibrium population achieved at infinite time,

$$n_1(\mathrm{})=\frac{n_HP_{21}}{P_{21}+P_{12}}.$$
(4)
Thus there is only one time scale for for the approach to ionization equilibrium

$$\tau _{relax}=1/\left(P_{21}+P_{12}\right).$$
(5)
It is wrong to talk about separate time scales, one for ionization and one for recombination, even though there are distinct ionization and recombination rates. Also, this result of an exponential approach to a final equilibrium state with a given time scale assumes that the ionization and recombination rates per atom are constant in time. In the real case where the rates per atom are themselves evolving (because of changes in the radiation field and electron density) the concept of a time scale is not well defined, although it is still useful.

When the rate of upward radiative transitions balances the rate of downward radiative transitions, $`n_1R_{12}=n_2R_{21}`$, the equations can be simplified by assuming that this holds exactly, which is called detailed balance. In this case, $`n_iP_{ij}n_iC_{ij}`$ and the solution for the population of level one is

$$n_1(t)=n_1(\mathrm{})+\left(n_1(0)n_1(\mathrm{})\right)e^{t(C_{21}+C_{12})},$$
(6)
where now the final equilibrium population is

$$n_1(\mathrm{})=\frac{n_HC_{21}}{C_{21}+C_{12}},$$
(7)
and the relaxation time scale is

$$\tau _{relax}=1/\left(C_{21}+C_{12}\right).$$
(8)
This assumption of detail balance thus implies that the radiative rates per atom change with time as the populations relax to their equilibrium state, so that

$$\frac{R_{12}}{R_{21}}=\frac{n_2}{n_1}=\frac{n_2(\mathrm{})\left(n_1(0)n_1(\mathrm{})\right)e^{t(C_{21}+C_{12})}}{n_1(\mathrm{})+\left(n_1(0)n_1(\mathrm{})\right)e^{t(C_{21}+C_{12})}}.$$
(9)
Another possible simplification is to express the radiative rates in terms of the net radiative bracket (NRB),

$$n_2R_{21}n_1R_{12}=n_2R_{21}\left(1\frac{n_1}{n_2}\frac{R_{12}}{R_{21}}\right)=n_2R_{21}NRB.$$
(10)
In this case the equation for the level 1 population becomes

$$\frac{Dn_1}{Dt}=n_2\left(C_{21}+R_{21}NRB\right)n_1C_{12}.$$
(11)
If it is assumed that both the rates per atom and the NRB are constant, then the solution for the population is

$$n_1(t)=n_1(\mathrm{})+\left(n_1(0)n_1(\mathrm{})\right)e^{t(C_{21}+C_{12}+R_{21}NRB)},$$
(12)
where now the final equilibrium population is

$$n_1(\mathrm{})=n_H\frac{C_{21}+R_{21}NRB}{C_{21}+C_{12}+R_{21}NRB},$$
(13)
and the time scale to approach equilibrium is

$$\tau _{relax}=1/\left(C_{21}+C_{12}+R_{21}NRB\right).$$
(14)
If the net radiative transition rate is small compared to the upward and downward rates individually, then the NRB will be very small, so the radiative rates will make only a small change in the relaxation time scale given by the collisional rates. The radiative transitions can either increase (if NRB $`<`$ 0) or decrease (if NRB $`>`$ 0) the relaxation time scale.

We have calculated the relaxation time scale (at each height and each time step) from our numerical simulation. We proceeded in the following way: The value of the hydrodynamic variables were taken from a given time step and kept constant in time for the relaxation time scale calculation. The population densities consistent with this state of the atmosphere were calculated by solving the equations of statistical equilibrium. This defined the initial population density in the ionized state, $`n_p(0)`$. This atmosphere was then perturbed by increasing the temperature by 1% throughout. The populations consistent with this perturbed state was calculated from the equations of statistical equilibrium giving the asymptotic solution, $`n_p(\mathrm{})`$. The time evolution from the initial state of the number of protons at a given height was also calculated using the full rate equations, defining $`n_p(t)`$. The numerical solution was cast in a 2-level form (see above) and the relaxation time scale was calculated from a least-squares linear fit to

$$\mathrm{ln}\left[\frac{n_p(t)n_p(\mathrm{})}{n_p(0)n_p(\mathrm{})}\right]=\frac{t}{\tau _{relax}}.$$
(15)
The relaxation time scale for hydrogen ionization/recombination, as found from the numerical simulation, from the photosphere to the transition region, is shown in Fig. 6 (thick solid line). The time scale increases outward from the $`1`$ sec. in the photosphere to $`10^5`$ sec in the mid chromosphere and then decreases to $`10^2`$ sec. at the base of the transition region.

The eigenvalues and eigenvectors of the rate matrix, $`P_{ij}`$ with $`P_{ii}=_jP_{ij}`$ (where $`P_{ij}`$ is the transition rate per atom from level $`i`$ to level $`j`$), help to clarify the processes controlling the time scale. The relaxation process can be represented by the equation

$$𝐧=\underset{i=0,5}{}c_i𝐯_𝐢e^{\lambda _it},$$
(16)
where $`𝐧`$ is the vector of the level populations,

$$𝐧=\left(\begin{array}{c}n_1\\ n_2\\ n_3\\ n_4\\ n_5\\ n_6\end{array}\right).$$
(17)
Here level 1 is the ground state and level 6 is the continuum. $`𝐯_𝐢`$ is the eigenvector corresponding to the $`i^{th}`$ eigenvalue $`\lambda _i`$. The coefficients $`c_i`$ depend on the initial conditions. There is a zero eigenvalue whose eigenvector is the equilibrium state. The other eigenvalues are all negative, since the populations are relaxing toward their equilibrium value.

The numerically determined relaxation time scale is compared with the time scale from the eigenvalue calculation, which is the inverse of the smallest (in absolute value) non-zero eigenvalue, in Fig. 6. We have calculated the eigenvalues using several different assumptions about the rate matrix. Note, first of all, that when all the processes are included in the rate matrix the time scale obtained from the eigenvalues (dotted line) is orders of magnitude smaller than found in the numerical solution (thick solid line). The time scale obtained from the rate matrix with only collisional rates (dot dash line) has the same general pattern as the numerical time scale, but is generally slightly larger. When radiative rates are added, but with all the Lyman transitions assumed to be in detailed balance, the time scale from the eigenvalues of the rate matrix (thin solid line) reproduce the numerical time scale up to the peak in the mid-chromosphere. This indicates that in this region the large Lyman radiative rates are in fact changing with time so as to maintain nearly detailed balance as the populations change. The general behavior of the relaxation time is thus controlled by the collisional processes. The time scale increases from the photosphere to the mid chromosphere, where the electron density has a minimum approximately as the inverse of the collisional recombination rate ($`n_e^2`$). From the mid-chromosphere to the transition region the time scale decreases, first because the electron density (thick dashed line) increases as the increasing ionization fraction of hydrogen more than offsets the outwardly decreasing overall density, and second because the temperature starts to increase due to absorption of coronal photons. The upward collisional rates are exponentially sensitive to the temperature. Radiative transitions generally reduce the time scale by a factor of about 4 from the all collision case. The wiggles in the relaxation time above its peak in the mid-chromosphere are nearly reproduced if the net Ly$`\alpha `$ radiative transition is included as a constant net radiative bracket (dashed line). Where Ly$`\alpha `$ photons begin to leak through the atmosphere they increase the relaxation time because in that region they produce a net increase in the population of the second level while collisions produce a net depopulation of the second level. Hence, the Ly$`\alpha `$ net radiative bracket is negative and as seen in the two level case this increases the relaxation time.

The role of the various transitions in the equilibration process can be determined from the eigenvectors and eigenvalues appearing in eqn. 16. We study the level in the atmosphere at column mass density $`m_c=10^4`$ in the initial atmosphere, and treat the Lyman transitions in detailed balance except for Ly$`\alpha `$ which is treated via its net radiative bracket. The slowest process, with eigenvalue $`2.2\times 10^6\mathrm{s}^1`$, has the eigenvector

$$𝐯_\mathrm{𝟏}=\left(\begin{array}{c}0.7\\ 2\times 10^7\\ 1\times 10^9\\ 4\times 10^{10}\\ 3\times 10^{10}\\ 0.7\end{array}\right).$$
(18)
which thus represents the relaxation between the ground state and the continuum. The next slowest process, which is 10 orders of magnitude faster, with eigenvalue $`2\times 10^4\mathrm{s}^1`$, has the eigenvector

$$𝐯_\mathrm{𝟐}=\left(\begin{array}{c}2\times 10^4\\ 0.7\\ 4\times 10^3\\ 1\times 10^3\\ 9\times 10^4\\ 0.7\end{array}\right).$$
(19)
and is relaxation between the second level and the continuum. The fastest relaxation processes, all with eigenvalues of order $`1\times 10^7\mathrm{s}^1`$, are between level 2 and levels 3, 4 and 5. In the chromosphere, unlike the transition region, ground state relaxation does not take place via direct level 1 $``$ continuum transitions, but rather through a level 1 $``$ 2 transition followed by level 2 $``$ higher levels and continuum transitions (Fig. 3). Thus the process of hydrogen ionization in the chromosphere can be thought of as a very rapid (small fraction of a second) equilibration of the 2nd and all higher levels with the continuum, together with a slow leakage of electrons from the ground state to the 2nd level, or visa versa, depending on whether hydrogen is ionizing or recombining.

## 5 Consequences of Slow Rates

The ionization structure of a dynamic atmosphere is very different from that of a static atmosphere. The ionization/relaxation time scale decreases dramatically (to about 10 - $`10^3`$ sec.) in chromospheric shocks, where the temperature increases significantly and to a lesser extent in the elevated temperature tail of the shocks (Fig. 7). As a shock propagates upward it strengthens and the time scale becomes shorter. However, the time scale is still too long for the ionization to reach its equilibrium value in the shock and the peak ionization occurs behind the shock front (Fig. 7). The slow ionization results in more of the shock energy going into raising the temperature rather than ionizing the gas (Carlsson & Stein, 1992).

The primary processes controlling the ionization in shocks is the same as for the static atmosphere: photoionization from the Balmer level and photorecombination to the higher levels, with slow leakage between the ground state and the 2nd level. The high temperature (and density) in the shocks increases the collisional and radiative leakage rates between the ground state and the 2nd level which leads to a smaller relaxation time. Balmer photoionization has the greatest increase, both within the shock and in the post shock tail (Fig. 5). In the shock itself, the sum of all the ionization/recombination rates is approximately equal to the Balmer photoionization rate alone. That is, almost all the ionization is occurring as Balmer photoionization and there is little balancing recombination. Thus the ionization is not in statistical equilibrium and there is net ionizing in the shock. In the mid-chromosphere there is net ionization even in the post shock tail, though at a slower rate. When the shocks reach the upper chromosphere, increased recombination to the third and higher levels leads to net recombination in the post shock tail, but at a much slower rate than the initial ionization in the shock (Fig. 5, t= 1650 s). Because the time scale is short in the shock where hydrogen is ionizing, due to the high temperature, but longer in the post shock region, where the temperature decreases and hydrogen is recombining, the steady state level of hydrogen ionization tends toward the value approximating the peak ionization in the shocks (Fig. 8).

In the photosphere (e.g. initial height of 0.4 Mm), the dynamical ionization fraction follows the statistical equilibrium one with a delay of 20-40 seconds and the ionization fraction is about $`10^5`$ (Fig. 8). By the low chromosphere (initial height of 0.6 Mm) the rates become too slow to keep up with the dynamic variations of density and temperature and the ionization fraction shows much smaller variations than in equilibrium. The ionization fraction increases with time from $`2\times 10^5`$ to a steady state value of $`10^4`$. In the mid-chromosphere (initial height 1.0 Mm) the rates are so slow that the ionization variation does not follow the dynamics at all. There is a slow, secular increase in the steady state ionization fraction from $`10^4`$ to $`10^2`$. At the same values of the hydrodynamic variables the statistical equilibrium ionization varies by six orders of magnitude. The behavior slightly higher (at 1.4 Mm) is similar. The slow ionization/recombination rates thus cause the ionization fraction to vary much less than in equilibrium and the mean ionization fraction represents the conditions at the peaks of the shocks and is substantially higher than the mean equilibrium value.

The small amplitude of the electron density variations compared with their equilibrium values and the fact that the mean electron density samples the peaks of the shocks rather than the mean conditions have several consequences for the proper interpretation of chromospheric diagnostics.

The intensity of collisionally excited lines will vary much less than the hydrodynamic variations would imply. An analysis based on equilibrium values of the electron density would give too low amplitudes for the atmospheric variations even if the effects of departures from LTE are taken fully into account.

Classical static models of the solar chromosphere are based on temporal and spatial averages of intensity — either continuum intensities shortward of the silicon edge at 152 nm (like in the models by Avrett and co-workers, e.g., Vernazza et al. (1981), Maltby et al. (1986), Fontenla et al. (1993)) or line profiles of resonance lines from ionized calcium and magnesium. The mean is taken in the ultraviolet part of the spectrum where the temperature dependence of the Planck function is more exponential than linear. Although the source function is quite decoupled from the Planck function and has a smaller amplitude than the local Planck function, the mean source function preferentially samples the high temperatures in shocks. This was dramatically illustrated in Carlsson & Stein (1995) where the best semi-empirical fit to the mean intensities in the simulation showed a chromospheric temperature rise while the mean temperature showed no increase.

At long wavelengths the Planck function has a linear dependency on temperature and one avoids the effect of averaging a non-linear function (at ultraviolet wavelengths) that gives an exaggerated temperature increase with height. The temporal average of the intensity at mm and sub-mm wavelengths should give the mean temperature as function of height provided the formation is in LTE. However, the electron density sampling of shock conditions at long wavelengths has a similar effect as temperature sampling of shocks at short wavelengths and an analysis based on equilibrium values will give an exaggerated temperature increase with height. The current simulations do, in fact, match observations made at long wavelengths remarkably well (Loukitcheva et al., 2001). The conclusion from the simulations is that you can not construct a model of the mean atmosphere from temporally averaged intensities.

The result that the hydrogen ionization/recombination is slow is robust. As detailed in section 2 the present simulations do not contain several important physical ingredients for the proper modelling of the upper chromosphere. However, the shocks propagating through the chromosphere span all reasonable chromospheric conditions and the time scale for ionization/recombination is always longer than the dynamical time scale. Even the VALIIIC semi-empirical model atmosphere with a temperature rise already at 0.5 Mm height and a rather high temperature throughout the chromosphere gives an ionization/recombination time scale around 103 seconds (Fig.7).

## 6 Conclusion: The Dynamic Chromosphere

The solar chromosphere is a dynamic region. This is obvious both from basic properties of the solar atmosphere and from observations. First, the solar photosphere is continually perturbed by convection. The large drop in density between the photosphere and the base of the transition region (14 scale heights) means that any disturbance in the photosphere that propagates upward will grow in amplitude with height in order to conserve energy unless it is strongly dissipated. Thus the chromosphere will experience large amplitude motions driven from the photosphere. Second, chromospheric spectral lines such as CaII H and K and Mg h and k and the ultraviolet continuum show significant temporal and spatial variability on a wide range of periods and spatial scales.

Observations of such dynamic regions can not be properly analyzed based on static model atmospheres. There are four robust results from the analysis of this paper that exemplify this basic truth and are independent of the details of the specific model: First, the hydrogen ionization and recombination rates under any reasonable chromospheric conditions are slow compared to the rate of dynamical changes there. As a result, the ionization/recombination time scale is longer than the dynamical time scale in the chromosphere. Second, the fluctuations in the hydrogen ionization are smaller than the variation of its statistical equilibrium value calculated from the instantaneous conditions. Because the ionization/recombination time scale is longer than the dynamical time scale, there is never enough time for the ionization state to reach its equilibrium value. Third, the hydrogen ionization is greater than the statistical equilibrium value for the mean atmosphere. Because the ionization and recombination rates increase with increasing temperature and density, ionization in shocks is more rapid than recombination behind them. Therefore, the ionization state tends to represent the higher temperature of the shocks. Finally, the mean value of a dynamic property is not the same as that property evaluated for the mean atmosphere. Where adjustments are not instantaneous, a property depends on the history of the atmosphere as well as its instantaneous state. In addition, observed properties are in general non-linear functions of the state of the atmosphere. For both reasons $`<f(x)>f(<x>)`$, for a property $`f`$ of a state $`x`$.

## Acknowledgments

RFS was supported by NASA grant NAG 5-9563 and NSF grant AST 9819799, for which he is grateful. This work was supported by the Norwegian Research Council’s grant 121076/420, “Modeling of Astrophysical Plasmas”.