Charlotte Froese Fischer—Her Work and Her Impact

Abstract

Charlotte Froese Fischer has been at the forefront of research in atomic structure theory for over 60 years. She has developed many of the methods currently used by researchers and has written associated computer programs which have been published and hence made accessible to the research community. Throughout her career, she has consistently encouraged and mentored young scientists, enabling them to embark on independent careers of their own. This article provides an overview of the methods and codes she has developed, some large-scale calculations she has undertaken, and some insight into the impact she has had on young scientists, and the leadership she continues to show as she reaches her 90th birthday.

1. Foundations

Charlotte Froese was born in the Ukraine, but soon after her birth, her family moved, via Germany, to Canada. She grew up near to Vancouver, and undertook her undergraduate and masters studies, in Mathematics along with Chemistry and some Physics, at the University of British Columbia (UBC) in Vancouver [1]. From there, she moved to the University of Cambridge, England, to undertake studies for her PhD under the supervision of Douglas Hartree. This interaction provided Charlotte (and therefore us) with a link to the very beginnings of atomic structure calculations. These had begun with the foundational paper of Schrödinger (1926) [2], which also dealt with the very simple case of hydrogenic ions. This was extended by Hylleraas (1928) [3] to the case of the ground state of helium, in which he incorporated the interelectronic distance $r_{12}$ explicitly in a form of the wave function which also contained variational parameters, determined by minimising the total energy. The use of interelectronic coordinates, while possible for very simple atomic systems, did not lend itself to extension to the calculations for many-electron atoms and ions. Instead, Hartree (1927, 1928) [4,5] proposed a method which was capable of application to any number of electrons in an atom or ion. In this, the N-electron wave function is represented by a product of N one-electron functions, or orbitals, with the motion of an individual electron determined by a single orbital.

$$\tilde{\psi }(1,2,..,N)=u_1\left(1\right)u_2\left(2\right)...u_N\left(N\right)$$

(1)

He used physical arguments to demonstrate that the orbital functions $u_i\left(i\right)$ should satisfy the non-linear equations

$$\left(-\frac{1}{2}{\nabla }_i^2-\frac{Z}{r_i}\right)u_i\left(r_i\right)+\sum _{j\ne i}{u}_i\left(r_i\right)\int \frac{|u_j\left(r_j\right){|}^2}{r_{ij}}d{\tau }_j={\epsilon }_i{u}_i\left(r_i\right)$$

(2)

He chose to take a spherical average of the final term of the potential so that the angular dependence of the orbital functions took the form of single spherical harmonics, $Y_l^m$, and also the radial parts of the orbital functions were assumed to be independent of m. The number of unknown functions was thereby substantially reduced. Nevertheless, the non-linearity of the equations meant that they had to be solved iteratively, with the aim of achieving self-consistency.

Hartree’s simple product form of the wave function did not satisfy the anti-symmetry requirement for the wave function. Fock (1930) [6] therefore wrote the wave function as an anti-symmetrised product of orbitals, represented by a determinant:

$$\Psi \left(LS\right)=\frac{1}{\sqrt{N!}}\left|\begin{array}{ccc}{\varphi }_1\left(1\right)& \cdots & {\varphi }_N\left(1\right)\\ ⋮& \ddots & ⋮\\ {\varphi }_1\left(N\right)& \cdots & {\varphi }_N\left(N\right)\end{array}\right|$$

(3)

in which the one-electron functions (spin-orbitals) now incorporated a function describing the spin of the electron:

$${\varphi }_i(r_i,m_{si})=u_i\left(r_i\right)\left[\alpha \left(i\right)\mathrm{or}\beta \left(i\right)\right]$$

(4)

with $m_s$ = $\pm \frac{1}{2}$ respectively.

This led to more elaborate equations—the Hartree-Fock equations—particularly because of the additional term in the potential, representing the possible exchange of the indistinguishable electrons:

$$\begin{array}{c}\hfill \left(-\frac{1}{2}{\nabla }_i^2-\frac{Z}{r_i}\right)u_i\left(r_i\right)+\sum _{j\ne i}[u_i\left(r_i\right)\int \frac{|u_j\left(r_j\right){|}^2}{r_{ij}}d{\tau }_j-\delta (m_{s_i},m_{s_j})u_j\left(r_i\right)\int \frac{u_j^{*}\left(r_j\right)u_i\left(r_j\right)}{r_{ij}}d{\tau }_j]={\epsilon }_i{u}_i\left(r_i\right)\end{array}$$

(5)

Comparing the potential terms in Equations (2) and (5), that in Equation (2) has the same form as the first of the two terms in Equation (5). It is referred to as the direct term, for it comprises the direct interaction between electron i and the field arising from all the other electrons. Equation (5) contains an additional term, arising from the antisymmetric nature of the wave function. Compared with the first term, it can be seen that the function $u_i$ has been interchanged with one of the $u_j$, so that it appears inside the integral in the second term. Because of this exchanging of the placing of these orbital functions, this second term is referred to as the exchange term. Its physical interpretation is that it models the indistinguishability of electrons, a feature missing from Hartree’s formalism. In the case of Hartree’s Equation (2), the potential is local, in that the equations can be expressed as

$$\left(-\frac{1}{2}{\nabla }_i^2-\frac{Z}{r_i}\right)u_i\left(r_i\right)+V\left(r_i\right)u_i\left(r_i\right)={\epsilon }_i{u}_i\left(r_i\right)$$

where V, representing the potential, does not involve $u_i$. By contrast, such a multiplicative potential cannot be formed for Equation (5), since $u_i$ appears inside the integral in the second term. Therefore the potential in this case is non-local.

In both cases, the equations are non-linear in the orbital functions, so must be solved iteratively. Typically, at any iteration stage, the opening forms of the orbital functions are used to evaluate the integrals, and so the integro-differential equations become linear differential equations. For (2), each equation then involves just a single orbital function, whereas for (5), the exchange term means that the equations are coupled, involving all the orbital functions. In both cases, the orbitals comprising the solution at each stage of iteration become the input orbitals for the next stage, although in practice the orbitals obtained at earlier stages might also be used in the formation of the input orbitals, in order to speed up the convergence towards self-consistency.

Hartree set to work on solving, again self-consistently, these more elaborate equations, and although Fock’s formalism did not require it, he chose to represent the angular dependence of the space part of the orbital functions again by single spherical harmonics, and the radial functions as being independent of $m_l$ and $m_s$. The Hartree-Fock equations then became equations for the radial parts of the orbital functions.

For the next 25 years or so Hartree, alone (e.g., Hartree 1934 [7]) or with his father (Hartree and Hartree 1938) [8] and with other co-workers (e.g., Hartree, Hartree and Swirles 1939 [9]), determined wave functions for a variety of atoms and ions, with the radial functions being in tabular numerical form. It was a laborious undertaking, demanding careful book-keeping of the set of numerical functions, even though Hartree was later able to enlist some of the first computers that became available, and indeed which he was instrumental in designing [10].

2. Hartree-Fock Calculations

It was into this scientific environment that Charlotte came in the mid-1950s. Her first paper with Hartree concerned the solution of the Hartree-Fock (HF) equations for Ne IV and Ne V, (Froese and Hartree 1957 [11]). It was important to understand in detail the steps which had to be taken, in a systematic manner, for self-consistency to be achieved, in the solution of the HF equations. Charlotte’s mathematical background equipped her well for the tasks of studying the equations themselves and for programming the EDSAC computer in Cambridge. Sadly, it was to be her only paper with Hartree, for he died early in 1958, shortly after his definitive book on the subject (Hartree 1957 [12]) was published. His untimely death was a significant loss to the atomic structure community and of course to Charlotte herself. But, picking up his mantle, she proceeded to publish further papers during the same year ([13,14]), in which she completed the numerical solutions of the HF equations for a number of ions, including some of astrophysical importance. Over the next few years, further studies of the HF equations followed, but it is interesting to note her focus on aspects of the mathematical or numerical solution of the equations, for example on applications to high Z-values ([15,16,17]). By then, Charlotte had returned to UBC in Vancouver, where a new computer had just been installed—the first at UBC. She was then able to follow up the approaches she had developed with Hartree during her time in Cambridge, by investigating the numerical solution of the HF equations with emphasis on how electronic computers could be used most effectively, as well as considering the ways in which accuracy could be assured [18]. She was a pioneer in the use of electronic computers in atomic structure calculations, and was becoming a world leader in the field—a position she still holds!

The 1960s proved to be very significant for Charlotte. In 1964, she was awarded an Alfred P. Sloan Fellowship, a prestigious award given to those early-career researchers who show exceptional promise. She was the first woman to be given that award and accolade. Clearly, the foresight of the awarders has been fully vindicated!

Further HF calculations followed in those years, for a range of atoms and ions, of wave functions and energies, and other atomic properties such as oscillator strengths and hyperfine structure. In addition, issues affecting the rate of convergence of the iterative process in the solution of the HF equations as well as other mathematical aspects of the HF equations were studied and resolved.

Above all, she met and married Patrick Fischer, and in 1968 they moved across Canada to the computer science department of the University of Waterloo, Ontario. Since then, she has published under the name of Charlotte Froese Fischer.

3. Extensions of HF

This detailed study of the mathematics of the solution of the HF equations, as well as her development of an algorithmic approach, stood her in good stead for taking forward some of the many ways of extending the HF method in order to achieve greater agreement with experimental results, in energy differences, but also in the calculation of other properties such as oscillator strengths, hyperfine structure or isotope shifts.

Hartree’s original method as well as the HF method can be expressed in terms of a variational approach to calculating energies (Slater 1930 [19]). Consequently, HF energies are reasonably accurate, at least for isolated and low-lying states, since the errors in the energies are of second-order for first-order changes in the orbitals. However, in the case of other operators (for example the dipole operator encountered in the calculation of oscillator strengths), the errors in their matrix elements are of first order.

A powerful and now widely used means of improving on HF is through the use of wave functions which, in $LS$ angular momentum coupling, are of configuration interaction format:

$$\Psi \left(LS\right)=\sum _{i=1}^M{a}_i{\Phi }_i\left({\alpha }_{i}LS\right)$$

(6)

where ${\alpha }_i$ represents the coupling of the angular momenta of the orbitals in each of the configuration state functions (CSFs) ${\Phi }_i$, and typically ${\Phi }_1$ is chosen to be the HF wave function. For any choice of the form of the CSFs, the optimal values of the mixing coefficients $a_i$ are eigenvector components of the Hamiltonian matrix whose typical element is $<{\Phi }_i\left|H\right|{\Phi }_j>$, with the corresponding set of eigenvalues providing the calculated energy values. Froese (1964) [20] demonstrated the insufficiency of the HF process in her study of multiplet transitions in Fe XV and Fe XVI. The Fe XVI ion, is essentially hydrogenic with a single electron outside completely closed shells, and the HF method leads to reasonably accurate oscillator strengths. By contrast, for some transitions in Fe XV, HF is not adequate. Froese found that considerable improvement could be achieved by the inclusion of configuration interaction, particularly for ${}^1$P${}^{\mathrm{o}}$−${}^1$D transitions. Even if the lower state 3s3p ${}^1$P${}^{\mathrm{o}}$ is represented by just the HF configuration, there are two possible upper levels—3s3d and 3p${}^2$, and each state needs to be represented by a linear combination of the two ${}^1$D configurations, and both configurations have large mixing coefficients in each state. The radial functions of the orbitals were generated in HF calculations on different states, and then used in a configuration interaction, or more precisely superposition of configurations, calculation. Hartree, Hartree and Swirles (1939) [9] had incorporated configuration interaction much earlier, albeit with the simpler system of oxygen ions. Froese incorporated the mixing coefficients in the Hartree-Fock equations, thus providing an early example of the multi-configuration Hartree-Fock (MCHF) calculation.

In the MCHF process, the variational method is used to generate the MCHF equations, using a trial wave function in the form of equation (i6), and by setting to zero the first-order change in the energy expression $<\Psi \left|H\right|\Psi >$ subject to the orthonormality of the orbital functions. The MCHF equations are similar to the HF equations in that they are the equations which determine the orbital functions, but additionally incorporate the optimisation of the mixing coefficients, $a_i$ of Equation (6). The MCHF method includes the additional process of the diagonalisation of the Hamiltonian matrix, thus providing the optimal values of the mixing coefficients for the current forms of the orbital functions.

The key difference between MCHF and the superposition of configurations method is that in MCHF, the orbital functions are obtained as the solutions of the MCHF equations, whereas generally the term configuration interaction implies that the orbital functions are pre-determined separately, though with the same the diagonalisation of the Hamiltonian matrix to provide the CI mixing coefficients. For example, in an MCHF calculation of a light element such as Be I, the MCHF orbitals as well as the CI mixing coefficients would be determined directly by the two-stage process of solving the MCHF equations (to a pre-determined level of consistency) followed by the diagonalisation of the Hamiltonian matrix. The CI or superposition of configurations process might first fix those orbital functions which are `occupied’ in the HF approximation, with further orbitals, added from a variety of sources, to be used in the other CSFs. A bringing together of these somewhat different processes is sometimes undertaken for larger atoms, when, for example, the orbitals for some of the subshells might be fixed from simpler calculations, such as those of their HF forms or from an MCHF calculation using just the CSFs with substantial mixing coefficients, while the additional orbitals introduced for the other CSFs are treated as unknowns, to be determined by solving the MCHF equations. As a consequence, generally the MCHF wave functions are more accurate than CI wave functions, because of use of the variational principle in setting up the MCHF equations, though the difference can be fairly small.

A major application of the use of MCHF or CI wave functions is in the calculation of oscillator strengths of transitions for atoms or ions. The oscillator strength, or f-value, for electric dipole transitions is a dimensionless quantity, and therefore is usually evaluated in atomic units, and for an N-electron atom or ion can be expressed, in velocity form, as

$$f_v=\frac{2}{3}\frac{1}{g_i\Delta E}{\left|〈{\psi }_j\left|\sum _{k=1}^N{\nabla }_k\right|{\psi }_i〉\right|}^2$$

(7)

where $\Delta E$, given in atomic units, is the transition energy (the energy difference between the two states) and $g_i$ is the g-value of the energetically lower state, which means that in $LS$ coupling $g_i=(2L_i+1)(2S_i+1)$. The emission transition rate $A^{ji}$, sometime called the transition probability, is related to the absorption oscillator strength $f^{ij}$ by, in the case of electric dipole transitions,

$$f^{ij}=1.4997\times {10}^{-16}{\lambda }^2\frac{g_j}{g_i}{A}^{ji}$$

with A expressed in units of s${}^{-1}$, and $\lambda$ (in Å) is the wavelength of the transition. The equivalent length form of the oscillator strength is

$$f_l=\frac{2}{3}\frac{1}{g_i}\Delta E{\left|〈{\psi }_j\left|\sum _{k=1}^N{r}_k\right|{\psi }_i〉\right|}^2$$

(8)

Oscillator strengths appear in astrophysical modelling, for example in the determination of element abundances in stellar atmospheres, in the combination $log\left(gf\right)$ so oscillator strengths are frequently published as $gf$-values, with g as the g-value of the lower state/level of the transition.

When the wave functions of the two states involved are exact eigenfunctions of the Hamiltonian, so that the Hamiltonian commutes with each $r_k$, the length and velocity forms give the same result. For HF, the potential is non-local, so these conditions are not satisfied. Hence in most cases, the calculated length and velocity forms do not agree for HF wave functions. However, when a very simple local model potential is used, with the two wave functions each being exact eigenfunctions of this simple Hamiltonian, then length and velocity forms do agree, but the common value is not necessarily correct. So, while it is necessary for length and velocity forms to agree, their doing so is not a guarantee of accuracy. Rather, accuracy is achieved by studying the convergence of the two forms as the numbers of CSFs in expansions (6) of the CI or MCHF wave functions are increased. It is sometimes argued that the length value is the more reliable of the two, and so only the length is calculated or provided. It is indeed often the case that the length value is the more stable as wave functions are extended, as the velocity form is more affected by the degree of electron correlation included in the calculation. But to omit the velocity value is to remove one measure of accuracy. If the two forms converge at least closely to a common value, one can have confidence in the accuracy of that value. It is this convergence process that Charlotte has endeavoured to achieve through her calculations.

4. Some Illustrative Examples

Configuration interaction calculations were undertaken for the some challenging transitions in Al I (Froese 1965) [21] and for Si II (Froese and Underhill 1966) [22]. This work was extended [23] to include a consideration of the very different 3d radial functions when optimised on the two ${}^2$D states separately, and the effect of imposing conditions which ensure the orthogonality of the two ${}^2$D states. The results for Si II are shown in

Table 1. CI calculations for Si II.

CI Coefficients for ${}^2$D States of Si II
Configurations	3s3p${}^2$	3s${}^2$3d	3s${}^2$3d (orthog)
3s3p${}^2$	0.7908	−0.5629	−0.5263
3s${}^2$3d	0.5994	0.8016	0.8251
3p${}^2$(${}^1$S)3d	0.1146	0.1921	0.1977
3p${}^2$(${}^1$D)3d		0.0118
3p${}^2$(${}^3$P)3d		−0.0603	−0.0563
3s3d${}^2$	0.0467
$gf$-values from the 3s${}^2$3p ${}^2$P${}^{\mathrm{o}}$ ground state
Froese Fischer (1968) [23]	0.103		6.22
Froese Fischer (1981) [24]	0.006		6.83
Hibbert et al. (1992) [25]	0.011		6.69

.

In the first part of the table, it is clear that the two main ${}^2$D configurations have strong components in each state. In an extension of [22], Froese Fischer [23] found state the mean radius of the optimal MCHF 3d function when optimised on the 3s3p${}^2$ state was 3.215, whereas when optimised on the 3s${}^2$3d state, the mean radius of 3d was 5.412. When the orthogonality constraint is imposed, some changes in the mixing coefficients occur, but the strong CI mixing in both states is generally maintained. The strong CI mixing has a pronounced effect on the oscillator strengths (shown in Table 1 as $gf$ values, with g = (2L + 1)(2S + 1) for the state with lower energy = 6). The oscillator strength between the ground state and the 3s${}^2$3d state is enhanced by the CI mixing, whereas that in the transition to the 3s3p${}^2$ state, strong CI cancellation occurs, resulting in an abnormally small $gf$ value. In a much more extensive calculation of these transitions, involving many configurations, including those representing core polarisation [24], and the use of non-orthogonal orbitals to circumvent the problem of different 3d functions being optimal for the two ${}^2$D states, the cancellation is stronger still. That very small oscillator strength agrees well with an independent calculation of my own [25]. This transition is of astrophysical significance in the study of the abundance of Si in the interstellar medium, Shull et al. (1981) [26]. They determined, from observations, a $gf$ value of around 0.033, not quite in agreement with the most recent calculations but, given the extent of the cancellation, it can be seen that calculated values and results derived from observation are fairly close.

One important feature of Charlotte’s MCHF calculations is the demonstration of convergence of results as the CI expansions of the wave functions are increased in size. While there is a variational principle which ensures that a longer expansion leads to a monotonic lowering of energy for a particular wave function, there is no such guarantee of monotonic improvement either for energy differences or for oscillator strengths. One way of assessing the accuracy of results is to see the way in which these quantities change as wave function expansions are extended in a systematic manner. Such a systematic analysis is provided by Tong et al. (1995) [27] for low-lying quartet transitions in neutral nitrogen. The basic way of providing a systematic analysis is first by constructing an active set of correlation orbitals usually characterised by specifying the maximum n-value and allowing all possible orbitals up to that maximum value. Ideally, the active space should consist of all possible CSFs which can be formed using the active set of orbitals, but this would be a prohibitively large calculation. Instead, Tong et al. [27] undertook a development of a sequence of increasingly complex models. in which the active set is systematically increased.

The most challenging transition is 2s${}^2$2p${}^3$${}^4$S${}^{\mathrm{o}}$− 2s2p${}^4$${}^4$P, because of quite severe CI cancellation in the transition integral.

Some results are shown in

Table 2. Oscillator strengths for the 2s${}^2$2p${}^3$${}^4$S${}^{\mathrm{o}}$− 2s2p${}^4$${}^4$P transition in N I [27].

Model Number	Configuration Complexes ${}^{†}$	$\mathsf{\Delta }\mathit{E}$	${\mathit{f}}_{\mathit{l}}$	${\mathit{f}}_{\mathit{v}}$
1	{2}${}^3$ {2,3,...,6}${}^2$	87,271	0.3513	0.3628
2	{2}${}^2$ {2,3}${}^1$ {2,3,...,6}${}^2$	88,524	0.0533	0.0563
3	{2}${}^1$ {2,3}${}^2$ {2,3,...,6}${}^2$	88,375	0.0667	0.0693
3+	{2}${}^1$ {2,3}${}^2$ {2,3,...,7}${}^2$	88,356	0.0658	0.0687
Within model 3	{2}${}^1$ {2,3}${}^2$ {2,3}${}^2$	89,760	0.3163	0.4062
	{2}${}^1$ {2,3}${}^2$ {2,3,...,4}${}^2$	89,324	0.1108	0.1171
	{2}${}^1$ {2,3}${}^2$ {2,3,...,5}${}^2$	88,446	0.0701	0.0717
	{2}${}^1$ {2,3}${}^2$ {2,3,...,6}${}^2$	88,375	0.0667	0.0693
	{2}${}^1$ {2,3}${}^2$ {2,3,...,7}${}^2$	88,356	0.0658	0.0687
Exp. [28] Average over J		88,132

† : The notation indicates the orbital occupancy of the 5 outer electrons; for example {2}${}^2$ {2,3}${}^1$ {2,3,...,6}${}^2$ means two orbitals with n = 2, plus one from n= 2 or 3; plus two with any n value from 2 through to 6; additionally, all CSFs contain 1s${}^2$.

, demonstrating ways in which the convergence of results as the wave functions are increased in complexity can be studied. One of the features of this transition is that the 2s2p${}^4$${}^4$P state interacts very strongly with two other ${}^4$P states—2s${}^2$2p${}^2$3s and 2s${}^2$2p${}^2$3d—so that all three states must be treated in an equivalent manner. In the first part of the table, the increase in complexity is achieved by extending the type of correlation effects included. Thus, model 1 includes single and double replacements of n = 2 orbitals by an increasing number (up to n = 6) of correlation orbitals (the active space) in the reference set comprising 2s${}^2$2p${}^3$ and 2s2p${}^4$. Correlation effects are therefore not included in an equivalent way in the three ${}^4$P states. This defect is corrected in model 2, which includes configurations 2s${}^2$2p${}^2$3l in the reference set. Model 3 allows for CSFs with only one orbital from the {2s,2p} set. In this analysis, a major improvement in the oscillator strength is achieved through model 2. It is worth noting that even in model 1, the length and velocity forms of the oscillator strength are in fairly good agreement (they differ only by about 3%), but the common value is clearly incorrect. In the second half of the table, we see how Tong et al. [27] undertook a systematic analysis of oscillator strength and transition energy values as the size of the active set is increased, within model 3. At each stage, defined by the value of the largest n-value of the orbitals included in that stage, all orbitals are computed using the MCHF program. Only when the n = 5 orbitals does any sense of convergence appear. The change between n = 6 and n = 7, the latter denoted by model 3+ in the first part of the table, is quite small. This analysis is characteristic of Charlotte Froese Fischer’s approach towards achieving confidence in the accuracy of the results of her calculations.

5. Computer Programs for Atomic Structure

The publication of her HF code [18] was but the first example of her adoption of the principle that if computer codes were of wide applicability, they should be made available to other users and not simply retained for personal use. There was clearly a need to create a computer program for the MCHF method, similar to that of the HF method [18], in which the processes were again automated once the CSFs and the initial estimates for the radial functions of the orbitals were selected, and which could readily be extended to any number of CSFs, or to any atomic system. During the 1960s, Charlotte developed such a code. The energy functional was based on Equation (6), from which the MCHF equations were derived using the variational principle. As with the HF process, these equations were coupled integro-differential equations, but now the mixing coefficients $a_i$ were included. Again, the solution comprised an iterative process. At each iteration, the integrals in the MCHF equations were calculated from the previous iteration of the radial functions, and also the mixing coefficients were taken from the previous iteration; the resulting (coupled) differential equations were then solved to give a new set of radial functions, and with these new functions, the new mixing coefficients were obtained by diagonalising the Hamiltonian matrix, although in early versions of the MCHF code only one eigenvalue/eigenvector was allowed, with the eigenvector components computed iteratively. Thus the orbital radial functions were obtained directly from these MCHF equations, rather than from separate HF solutions, followed by a superposition of configurations process, as had been used in earlier calculations. Hence the MCHF orbitals themselves incorporated the effects of configuration interaction.

In 1969 a new journal, Computer Physics Communications (CPC), was launched, with my colleague in Belfast, Phil Burke, as Principal Editor and Charlotte as one of the subject editors for atomic structure. This journal provided the vehicle for the dissemination of the MCHF code (Froese Fischer 1969 [29]), updated a little later with the name MCHF72 (Froese Fischer 1972 [30]), and again later still (Froese Fischer 1978 [31]). This journal provided assurance for both authors and users of the published codes. Submissions to the journal required a description of the methods underlying the codes, as well as details of how the codes were constructed—for example, the subroutines, procedures, modules from which the code was built and how they linked together. In addition, sample input data had to be provided, as well as the corresponding output and a statement of which type of computer had been used in obtaining the output, and the size of the program. The codes and the papers describing them were then rigorously refereed before acceptance for publication. Users could then compare the output from the test data when they ran the code on their own machine against the output submitted with the code. They could then be reasonably well assured that the code of interest was what they required, and that it worked as expected. The authors had much more to do in preparing the paper and code for publication than is customary for a standard research paper, but having done so, they were much less likely to be asked for a copy of the code, or questioned about its reliability than if they had distributed the code privately. The underlying assumption was (and is) that the authors’ codes become widely available for others to use.

There are benefits to authors as well as safeguards for users. For example, their codes become known and used by a wider community, so that due recognition of the authors’ work is accorded. Unexpected or unforeseen `bugs’ might be eliminated at the journal’s checking stage. Requests for copies of the codes can be directed to the journal.

One additional benefit I have found, both as a user and authors, is that users and authors might and sometimes do work together on further code developments. It was this aspect that led me to begin working actively with Charlotte Froese Fischer. The MCHF code aims to solve the integro-differential equations for the radial functions of the orbitals. However, the Hamiltonian matrix elements in the energy functional can be expressed as a weighted sum of radial integrals. The integrals over the angular and spin coordinates can be achieved exactly and these data form part of the input to the code. My own code, WEIGHTS, also published in the first volume of CPC (Hibbert 1970 [32]), together with a follow-up the following year (Hibbert 1971 [33]), provided the weighting coefficients for the two-electron part of the Hamiltonian, while a later code (Hibbert 1974 [34]) provided the same data for the one-electron operators. This led to a strong collaboration between Charlotte and myself, which has continued from time to time ever since.

In the early years, it was not self-evident that the new journal would be a success. Some sceptics doubted whether authors would be willing to spend the time necessary for a thorough description of their codes, or to subject them to the rigorous refereeing scrutiny that was required. Others wondered how willing users would be to familiarise themselves with the detailed working of the codes. Their fears were unfounded, as is evidenced by the fact that as the 50th anniversary of the journal arrives, the volume count is almost at 250. This success is due in no small measure to the enthusiasm and determination of Phil Burke, Charlotte Froese Fischer and others who were well established researchers in their fields and whose involvement gave great credence to the value of the journal and to the quality of the papers and programs it would be publishing: their names alone gave the new journal considerable credibility. Over the next 50 years or so, Charlotte was to publish around 40 programs or procedures in CPC, and throughout that time, she maintained full support for the principles on which the journal was developed, and of course her codes were very robust and reliable.

6. Extensions and Enhancements of MCHF—Non-Relativistic Treatment

We consider here the extensions to the MCHF process, but still in the context of calculations in $LS$ angular momentum coupling. We will report on a relativistic approach in the following section.

6.1. Non-Orthogonal Orbitals

In using the MCHF process, it is customary to express all the atomic states in terms of a common set of orbitals which are orthogonal to each other. But this is a restriction, which for elements with only a small number of electrons outside closed shells can sometimes be overcome with additional configuration interaction (essentially by extending the set of orbitals), but in other cases even very extensive additional CI would be required and possibly the restriction cannot fully be overcome even then. It is therefore necessary to consider the possibility of using orbitals which, for the same l-value, are not orthogonal to each other (orbitals with different l-values are of course mutually orthogonal).

Some exemplars of where non-orthogonal orbitals are useful include:

• He: 1s${}^2$ and 1s2p—radially, the 1s function of the 1s2p state is close to being hydrogenic whereas the 1s function for the ground state resembles a screened hydrogenic function.
• Be: [1s${}^2$]2s${}^2$${}^1$S, 2s2p ${}^3$P${}^{\mathrm{o}}$ and ${}^1$P${}^{\mathrm{o}}$—the 2s functions differ somewhat from state to state, but the more significant feature is that the mean radii of the 2p functions in the two excited states differ by around a factor of two.
• Al-sequence: we have already noted in Table 1 that the optimal 3d function in the two lowest ${}^2$D states is very state dependent. A more appropriate CI expansion would have configurations of the form:
  ${}^2$D: 3s${}_1^2$3d; 3p${}_1^2$(${}^1$S)3d, 3p${}_2^2$(${}^1$D,${}^3$P)3d${}_1$, 3s${}_2$3p${}_3^2$
  where the same $nl$ orbital but with different subscripts need not be mutually orthogonal.
• The 3d orbital in open d-shells, for example in the iron group elements, can be very term-dependent even for an individual ion.

The lifting of the restriction that the same orbital set, comprising mutually orthogonal functions, be used for all states can lead to substantial improvements in the accuracy of the results, and/or much shorter CI expansions to achieve comparable effects. For example, even in the light element neutral nitrogen, Robinson and Hibbert (1997) [35] found that just a few configurations could achieve for quartet transitions an accuracy as good as, and in some cases much better—comparing length and velocity values—than could be obtained with a much larger calculation using orthogonal orbitals (Hibbert et al. 1985) [36]. The difficulty was obtaining agreement between the length and velocity forms, and although the calculations of Robinson and Hibbert were not definitive, they did achieve comparable agreement with Tong et al. (1994) [27], who found that it was necessary to use of some thousands of CSFs. Some results are shown in Table 2.

To be able to undertake CI or more specifically MCHF calculations in the framework of non-orthogonal orbitals, I renewed my collaboration with Charlotte as we worked out how to modify the codes which undertake the angular and spin integrals, in order to incorporate the possibility of using non-orthogonal orbitals. This resulted in work which was published in CPC (Hibbert et al. 1988) [37], work which also began my collaboration with Michel Godefroid (Brussels). The extent of the non-orthogonality was limited: essentially we considered various pair correlations in any atomic state and allowed non-orthogonality of the orbitals in different pairs.

In order to study transitions, it was also necessary to allow for non-orthogonality between the orbitals in the two states or levels in the transition. In particular, the two wave functions used in the calculation of the oscillator strength might be calculated separately, using different orbitals in each case. This situation was allowed for by the use of a bi-orthonormal transformation of the orbital functions by which the methods applicable to orthogonal orbitals can be used, before calculating the transition matrix elements. The general theory of this process was introduced by Olsen et al. (1995) [38], in which the application to the MCHF codes was developed by Godefroid.

Its power was amply demonstrated by a study of the 2s${}^2$2p ${}^2$P${}^{\mathrm{o}}$− 2s2p${}^2$${}^2$D transition in boron. An example, given by Olsen et al. (1995) [38], of the success of the use of non-orthogonal orbitals, when combined with the bi-orthonormal transformation of the orbitals is shown in

Table 3. Oscillator strengths for the 2s${}^2$2p ${}^2$P${}^{\mathrm{o}}$− 2s2p${}^2$${}^2$D transition in B I [38].

Active Set: max n	$\mathsf{\Delta }\mathit{E}$(cm${}^{-1}$)	${\mathit{gf}}_{\mathit{l}}$	${\mathit{gf}}_{\mathit{v}}$
3	53,197	0.6876	0.8156
4	48,720	0.2456	0.2606
5	48,440	0.2625	0.2695
6	48,125	0.2891	0.2866
7	48,051	0.2928	0.2900
7 ${}^E$	47,847	0.2916	0.2912
Other results
Method		${\mathit{gf}}_{\mathit{l}}$	${\mathit{gf}}_{\mathit{v}}$
MCHF [39]		0.243	0.274
Expt: LIF [40]		0.283 ± 0.020

E: Using a weighted average (over J) of the individual levels given in [28].

for the 2s${}^2$2p ${}^2$P${}^{\mathrm{o}}$ – 2s2p${}^2$${}^2$D transition in neutral boron.

Table 3 exemplifies the process of making systematic improvements to the calculations. Specifically, the concept of active space of orbitals is used, consisting of all possible orbitals with an increasing value of n. The reference set of configurations consisted of (1s${}^2$2s${}^2$2p, 1s${}^2$2p${}^3$; ${}^2$P${}^{\mathrm{o}}$), (1s${}^2$2s2p${}^2$, 1s${}^2$2s${}^2$3d; ${}^2$D). The CSF set included all which are obtained by single and double replacements of the orbitals occupied in the reference set by any from the active set. The n = 3 results allow only for a single 3d orbital, whereas the optimal 3d from the 1s${}^2$2s${}^2$3d differs substantially from the optimal 3d of 1s${}^2$2s2p${}^2$, and that difference is not included at that stage of the calculation. This is substantially rectified by n = 4. The trend from n = 4 to n = 7 shows a systematic improvement in the transition energy $\Delta E$ and the agreement between the length and velocity forms of the $gf$-values. The final results, labelled 7${}^E$, are obtained from the n = 7 results by using the experimental rather than calculated transition energy. We also give in Table 3 the MCHF results using orthogonal orbitals [39] and the experimental result of O’Brian and Lawler [40] who used the laser-induced fluorescence method. The improvement obtained by using non-orthogonal orbitals can be clearly seen, and the final calculation of $gf$ lies well inside the experimental error bars. (Incidentally, the MCHF result of [39] is an example of the velocity form being closer than the length form to the experimental or converged calculated value of the oscillator strength, a counter-example to the view that the length value is the better of the two).

6.2. Use of B-Splines

The MCHF method focuses on the best possible way of introducing short-range electron correlation (assuming $LS$ coupling) into the solution of the Schrödinger equation. It works well when the active electrons (those in the outermost shells) have rather similar mean radii, so that electron exchange is a significant effect. Equally, it works well for excited, but not too highly excited, states. But for studying Rydberg series in atoms and ions, where the outer electron has a much larger mean radius than all the other electrons, then electron exchange has a lower probability and representation of highly excited states is not so easy to accomplish with the customary MCHF approach. Indeed, Hartree’s original method, which ignores electron exchange, becomes an improving approximation. One of the difficulties of the use of orthogonal orbitals in MCHF calculations is that an MCHF orbital $nl$ with a high n-value contains oscillations in the radial function arising from the requirement of orthogonality to those with lower n, which the use of non-orthogonal orbitals would not introduce. To overcome this difficulty, Froese Fischer and co-workers (e.g., Brage and Froese Fischer 1994 [41]) have modified the MCHF method by using B-splines to represent the orbitals, and particularly those of the outermost electrons. The number of B-splines required can be adjusted to encompass the radial range even of the outermost orbital. Each B-spline is a relatively localised function, being non-zero over only part of the radial range covered by the orbitals, so each orbital function must be represented by a sum of appropriate B-splines. While the B-splines are mutually orthogonal, the orbitals represented by these sums are not necessarily orthogonal, so a non-orthogonal approach is necessary. This flexibility can lead to more accurate representations of the wave functions. This use of B-splines is part of a vast field of study well described by Bachau et al. (2001) [42]. An example of the use of B-splines in atomic structure is given by Brage and Froese Fischer (1994) [41]. They studied several Rydberg series in neutral calcium. They allowed for some limited non-orthogonality of the outer orbitals and added the effect of core polarisation by means of a model potential. A small selection of their results is included in

Table 4. Binding energies (cm${}^{-1}$) of 4snp ${}^1$P${}^{\mathrm{o}}$ states in Ca II [41].

Label	Exp [28]	MCHF+BS	MCHF
4s4p	25,654	25,472	24,689
4s5p	12,574	12,684	12,160
4s6p	7627	7638	7060
3d4p	5372	5269	4609
4s7p	3881	3799	3247
4s8p	2826	2786	2405
4s9p	2122	2102	1853
⋯
4s22p	271	271

.

It can be seen that the inclusion of the effects of core polarisation, together with the use of a B-spline representation of the outer orbitals, (MCHF+BS), gives a substantially improved agreement with experiment for the binding energies, compared with the conventional MCHF approach. In particular, at the upper end of the Rydberg series, the agreement is excellent.

There are distinct similarities between this approach and that adopted in R-matrix calculations [43,44], where the outer (free) electron is customarily represented by a linear combination of basis set of continuum orbitals. Most R-matrix calculations to date have required the continuum orbitals to be orthogonal to the orbitals describing the bound orbitals of the N-electron core, but Zatsarinny and Froese Fischer [45] have undertaken an R-matrix calculation of the photoionisation of Li using B-splines and non-orthogonal orbitals.

7. Inclusion of Relativity

The MCHF method is essentially non-relativistic. But in order to study both allowed and forbidden transitions among levels, it was necessary for Charlotte to incorporate fine-structure effects into the calculation. Two approaches were available: a multi-configurational formalism based on the Dirac equation rather than Schrödinger’s equation, or an approximation to this approach in which the non-relativistic Schrödinger Hamiltonian is augmented by the operators of the Breit-Pauli Hamiltonian. For light atoms and ions, the latter process is adequate for the accuracy required in many applications. For heavy elements as well as for highly ionised systems with fewer electrons, the fully relativistic approach based on the Dirac equation is normally required.

7.1. Breit-Pauli Calculations

Since many of the applications being considered were for the lighter elements, Froese Fischer chose the Breit-Pauli approach. Her calculations followed the following pattern.

• The radial functions were optimised in an $LS$ MCHF calculation.
• The angular and spin integrals of the relativistic operators were evaluated using the Racah algebra analysis given by Glass and Hibbert (1978) [46] and are input to the MCHF+BP code.
• Then the full Breit-Pauli Hamiltonian matrix was diagonalised to give the $LSJ$ wave functions. These wave functions take the form

  $$\Psi \left(J\right)=\sum _{i=1}^M{a}_i{\Phi }_i\left({\alpha }_i{L}_i{S}_{i}J\right)$$

  (9)

  so that CSFs with different $L_i$ and $S_i$ can be combined to a common total J.

Initially, Charlotte studied fine structure separations of low-lying term energies and of the forbidden magnetic dipole and electric quadrupole transitions between them. For example, she calculated the splitting between the ${}^2$P${}_{0.5}^{\mathrm{o}}$ and ${}^2$P${}_{1.5}^{\mathrm{o}}$ of boron-like ions [47]. Some results are shown in

Table 5. Term splitting (cm${}^{-1}$) of the ground state in B-like ions.

Method	B	N${}^{2+}$	Ne${}^{5+}$	Si${}^{9+}$	Fe${}^{21+}$
MCHF+BP [47]	15.0	170.0	1292	6961	119,175
MCDHF ${}^{\mathrm{a}}$	15.7	172.4	1298	6968	118,177
Exp ${}^{\mathrm{b}}$	15.3	174.5	1307	6990	118,255

${}^{\mathrm{a}}$ Huang et al. (1982) [48]; ${}^{\mathrm{b}}$ NIST [28].

and compared with the fully relativistic treatment.

It can be seen that the Breit-Pauli approximation gives results which are very close to the fully relativistic results of [48] near the neutral end of the sequence but, as expected, diverge as Z increases, whereas the fully relativistic calculations are consistently close to the experimental values.

Similar calculations followed for the ground terms of other fairly light elements: C([49]), N([50]) and O([51]), and these papers also included the rates of forbidden transitions between the ground term levels. Around the same time, electric dipole transitions were studied, using the same MCHF+Breit-Pauli formalism [52], and this work permitted the study of intercombination lines. Of particular interest was the calculation of the 2s${}^2$${}^1$S${}_0$− 2s2p ${}^3$P${}_1^{\mathrm{o}}$ line in C III [53], as displayed in

Table 6. MCHF+BP transition rates of the 2s${}^2$${}^1$S${}_0$− 2s2p ${}^3$P${}_1^{\mathrm{o}}$ line in C III [53].

Degree of Correlation	Active Space	$\mathsf{\Delta }\mathit{E}$ (cm${}^{-1}$)	${}^1$P${}_1^{\mathbf{o}}$−${}^3$P${}_1^{\mathbf{o}}$	${}^3$P${}_2^{\mathbf{o}}$−${}^3$P${}_0^{\mathbf{o}}$	A (s${}^{-1}$)
Val ${}^a$	n = 3	52,746	51,592	76.79	89.3
	n = 6	52,733	50,684	77.08	95.6
+CP ${}^b$	n = 3	52,640	50,567	78.06	97.6
	n = 6	52,520	50,098	80.18	105.7
+CC ${}^c$	n = 3	52,362	50,948	77.33	91.9
	n = 6	52,343	50,230	79.53	103.1
CIV3 [54]		52,369	50,325	78.9	103.8
MCDHF [55]	n = 8	52,384	50,098	79.86	102.72
Experiment	[28]	52,391	49,961	80.05
	[56]				121 ± 7
	[57]				102.94 ± 0.14

${}^a$ Valence correlation only; ${}^b$ a + core polarisation; ${}^c$ b + partial core correlation.

.

As in other calculations undertaken by Charlotte, a systematic development of the results can be seen, both in the sense of growing complexity in the type of correlation included and, with each type, the variation in results as the active space is enlarged. The final result differs from experiment by about 20%, but is in complete agreement with another, independent calculation [54]. A later calculation [55] gave a transition rate of 102.72, which when extrapolated to take into account the slight inaccuracies in the ab initio energy separations resulted in a recommended value of 102.87, lying within the very narrow error bars of the recent heavy-ion storage ring experiment [57]. This calculation was updated by Froese Fischer and Gaigalas [58] to yield a transition rate of 103.0 s${}^{-1}$, with an estimated uncertainty of 0.4 s${}^{-1}$.

Many further calculations followed, for different ions, and this work culminated in the extensive tabulations of energy levels and electric dipole oscillator strengths by Froese Fischer and Tachiev (2004) [59] for the first row elements and their ions, and by Froese Fischer et al. (2006) [60] for the second row elements and their ions. These important compendia provide a set of accurate atomic data for transitions between a substantial number of levels of these elements. The data are characteristic of Charlotte’s work: they are undertaken in a consistent manner, and consider all the main correlation and relativistic effects appropriate for these elements. In the discussion in the papers, there is a strong attention to detail and where possible, comparison is made with experiment, especially the energy levels. They demonstrate why her work is considered to be first-rate and reliable for other researchers, including astrophysical modellers, to use with confidence.

7.2. Other Atomic Properties

Once MCHF+BP wave functions are available, it is possible to determine atomic properties other than transition rates. For example, Charlotte undertook calculations of isotope shifts and hyperfine interactions [61], photoionisation [62], and autoionisation [63].

7.3. Further MCHF-Based Computer Packages

By the 1980s, CPC had begun to include papers on computational methods, as well as continuing to publish computer codes and their descriptions. Accordingly, in order to pull together in one place several of the developments in MCHF and BP procedures which had taken place during the 1980s, Charlotte requested that almost the whole of one issue of CPC was devoted to these developments. The details are summarised in

Table 7. MCHF and Associated Codes.

Authors	Short Title	Type of Paper
Froese Fischer [64]	The MCHF atomic structure package	Methods
Froese Fischer [65]	MCHF support libraries and utilities	Package
Froese Fischer, Liu [66]	Configuration-state lists	Package
Hibbert, Froese Fischer [67]	Angular integrals with non-orthogonal orbitals	Package
Froese Fischer [68]	General MCHF program	Package
Hibbert, Glass, Froese Fischer [69]	Angular integrals for Breit-Pauli Hamiltonian	Package
Froese Fischer [70]	General CI program	Package
Froese Fischer, Godefroid, Hibbert [71]	Angular integrals for transition operators	Package
Froese Fischer, Godefroid [72]	Programs for $LS$ and $LSJ$ transitions	Package

.

These codes constitute a comprehensive package of codes allowing users to undertake a wide range of MCHF calculations. Once again, Charlotte was demonstrating her continuing commitment to making codes available to other users, codes which had been thoroughly tested.

8. Fully Relativistic Codes

The use of the Breit-Pauli Hamiltonian works well for lighter ions, but for transitions in heavy ions or for heavy elements, the use of fully relativistic wave functions becomes more accurate. In a multiconfigurational context, Desclaux (1975) [73] published the multiconfiguration Dirac-Fock (MCDF) code, based on a Dirac rather than Schrödinger formalism. This nomenclature of MCDF lacks the acknowledgement of the importance of Hartree’s work, and Charlotte preferred the designation MCDHF, the `H’ denoting `Hartree’. An enhanced version, GRASP, was published by Dyall et al. (1987) [74]. One of the authors (Parpia) spent some time working with Charlotte and as a consequence, Charlotte began to work with Grant’s group, culminating in an updated version of GRASP: Parpia et al. (1996) [75,76], with a further update in 2007 [77]. Further enhancements followed and were published in CPC, including the calculation of other properties using relativistic wave functions, such as hyperfine structure [78] and isotope shifts [79].

Some of the enhancements involved the use of more efficient methods for the calculation of the angular momentum integrals in GRASP. In this, Charlotte made use of the opportunity to collaborate with the group in Vilnius, whose work she had long admired [80]. The following 20 years saw the publication of many calculations of atomic properties for medium and heavy ions. Often, these were the outcome of an on-going international research team—CompAS—comprising Charlotte and groups based in Lund, Malmö, Vilnius, Krakow, Gdansk and Brussels. The group meets periodically to review progress and plan for future collaborations. I was present for part of their discussions in Lund in 2018 and in Brussels in 2019. I found Charlotte thoroughly involved in those discussions, an amazing degree of engagement at the age of 90.

9. In Summary

Charlotte Froese Fischer has been at the forefront of atomic structure developments for over six decades. In her early years in research, she developed programs for the calculation of atomic properties using electronic computers, which were only just becoming available. In subsequent years, the initial Hartree-Fock methods were extended, characterised by

HF → MCHF → MCHF+BP → MCHF+BS → Non-orthogonal → MCDHF → GRASP.

Charlotte has consistently been keen to use the most up-to-date computer architecture available to her, and to develop new numerical methods to exploit such facilities. This work involved the study of the convergence of MCHF and MCDHF iterations, and the efficient calculation of the angular momentum integrals (which are one of the most time-consuming parts of the calculations). Additionally, efficient methods and codes for the diagonalisation of huge matrices were developed. For large-scale MCHF or MCDHF calculations, with many thousands of CSFs included in the CI expansion of the wave functions, the Hamiltonian matrices are very large, but only a relatively small number of eigenvalues/eigenvectors, those which are the lowest in energy, are needed for subsequent calculations. The iterative approach, initially developed by Davidson [81], was further developed and programmed for this purpose [82].

Charlotte was part of the original team of scientists which set up the journal Computer Physics Communications in 1969. The publication of her codes in this journal has amply demonstrated her commitment to the ethos of the journal, not least in ensuring that others can make direct use of her work in undertaking their own calculations. That commitment was one of the reasons why the journal was able to develop in its early years, because other scientists could see the value of publishing their programs in that new medium. The ongoing success of the journal is a tribute to the commitment of Charlotte and the other editors, and of course to the long-term vision of its first Principal Editor Phil Burke.

As well as her original research papers, Charlotte has also published a number of textbooks, suitable for graduate students, which explain the methods she has developed and provide a detailed explanation of how calculations are undertaken. Amongst these are a discussion of the general Hartree-Fock method [83], and more recently an explanation of the computational approach to solving the MCHF equations and the calculation of atomic properties [84]. These books facilitate other scientists, including those starting out in the field, in developing an expertise in undertaking atomic structure calculations for themselves.

Charlotte is the holder of a number of awards recognising her international reputation as a world leader in the field of atomic physics. The Alfred P. Sloan Fellowship, awarded in 1964, was an early recognition of her promise. She was elected Fellow of the American Physical Society in 1991 for her contribution to the discovery of the calcium negative ion and for her extensive and innovative researches. Her research standing was also recognised overseas. In 1995, she was elected Fellow of the Royal Physiographic Society in Lund, Sweden. In retirement, her research continues, and in 2004 she was elected a foreign member of the Lithuanian Academy of Sciences. Then as recently as 2015 she was awarded an honorary doctorate by the University of Malmö, while in 2016 she was elected a Fellow of the Royal Society of Canada, especially appropriate as she returned to Vancouver.

But it is not just for her scientific contributions to the field of atomic structure, immense though they are, that Charlotte should be applauded. She has also had a large influence on the careers of many scientists, particularly on younger colleagues. As well as demonstrating her openness in making her computer codes available to the entire community, and in providing support for users, she has been keen to work with young scientists and to be a mentor to them, through the provision of encouragement as well as guidance on the standards of research expected in the field. And she has continued to collaborate with those who were young scientists and who have been able to develop their careers and their activity in the field. As Charlotte reaches her 90^th birthday, her enthusiasm and commitment continue unabated.

Charlotte Froese Fischer is a remarkable woman, and I consider it a privilege to be able to count her as a colleague and friend.