Targeted optimization in small-scale atomic structure calculations: application to Au I

The lack of reliable atomic data can be a severe limitation in astrophysical modelling, in particular of events such as kilonovae that require information on all neutron-capture elements across a wide range of ionization stages. Notably, the presence of non-orthonormalities between electron orbitals representing configurations that are close in energy can introduce significant inaccuracies in computed energies and transition probabilities. Here, we propose an explicit targeted optimization (TO) method that can effectively circumvent this concern while retaining an orthonormal orbital basis set. We illustrate this method within the framework of small-scale atomic structure models of Au I, using the Grasp2018 multiconfigurational Dirac–Hartree–Fock atomic structure code. By comparing to conventional optimization schemes we show how a TO approach improves the energy level positioning and ordering. TO also leads to better agreement with experimental data for the strongest E1 transitions. This illustrates how small-scale models can be significantly improved with minor computational costs if orbital non-orthonormalities are considered carefully. These results should prove useful to multi-element atomic structure calculations in, for example, astrophysical opacity applications involving neutron-capture elements.


Introduction
Recent advancements in astronomical spectral analysis have introduced new demands on atomic physics.A good example is the advent of the James Webb Space Telescope operating in the infrared spectral regime (Böker et al., 2022) where transition data for most elements and their ions are sparse.Another important example is the observation of the first-ever kilonova following the neutron star merger (NSM) event GW170817 detected by LIGO/Virgo (LIGO Scientific Collaboration and Virgo Collaboration et al., 2017).This ground-breaking discovery revived the debate on the origin of the neutron-capture elements, suggesting that these mergers are major sources of r-process nucleosynthesis in the universe (e.g.Cowan et al., 2021).
Reliable atomic data are increasingly required for accurate astrophysical modelling.For instance, when trying to identify individual spectral features, one needs complete atomic data and transition parameters with a high accuracy for each species considered.On the other hand, lower-accuracy data may suffice for modelling broad-band energy distributions and opacities; but this still requires completeness as well as some degree of accuracy, since too large uncertainties can prevent a detailed chemical analysis of the studied astrophysical material (Kasen et al., 2017).
Unfortunately, the current databases are both incomplete and of limited accuracy when it comes to some of the heavier r-process elements, in particular in the infrared spectral region (e.g.Smartt et al., 2017).This lack of information on atomic energy levels and processes is partly due to the complexities involved in carrying out accurate atomic structure calculations for many of these elements.
A significant challenge present in many atomic systems is the presence of multiple configurations with similar energies built from the same set of subshells.This is particularly true in the lanthanide group of elements where the ground and first excited configurations have closely aligned 4f, 5d, 6s and 6p orbital energies, implying that different subshell occupations lead to varying screening effects, and thus ideally requiring a non-orthonormal representation of the orbitals used to build the various configurations.
However, current state-of-the-art atomic physics codes focused on atomic structure, such as the configurationinteraction-based codes Fac (Gu, 2008) and Hullac (Bar-Shalom et al., 2001), assume a fully orthonormal orbital basis set.The current version of the multiconfigurational Dirac-Hartree-Fock (MCDHF) code Grasp ‡ (Froese Fischer et al., 2019;Jönsson et al., 2023a) has, since the Grasp2k release (Jönsson et al., 2007), allowed for some degree of non-orthonormality through biorthogonal transformations between even and odd parity orbital sets (Malmqvist, 1986), i.e. we obtain a different set of orbitals for each parity.Nevertheless, performing separate Grasp calculations of two competing configurations of the same parity may reveal significant non-orthonormalities between the orbitals optimized on the two different configurations.These non-orthonormalities, together with the possible strong configuration interaction among these states, can cause inaccurate expectation values, such as energies and transition rates, as well as wrong level ordering.The latter may be a particular issue for obtaining a grossly realistic representation of e.g.resonance lines and thermodynamics in plasma modelling (e.g. in kilonova modelling; Pognan et al. 2023).
In this work, within the framework of the latest development version of the Grasp code (Grasp2018), we discuss a method for conducting atomic structure calculations that address the issue of orbital non-orthonormalities, while also ensuring that the wavefunctions remain compact enough for efficient computations of complete spectroscopic, largescale opacity, or electron-scattering calculations that can be applied to a wide range of elements.The method, referred to as "explicit targeted optimization (TO)", consists of circumventing the orthonormality condition of Grasp by carefully optimizing the orbitals on targeted eigenstates.We also investigate semi-empirical rescalings of diagonal matrix elements to further improve the quality of the wavefunction, as recently discussed in the context of the Grasp code by Li et al. (2023).
We apply our methods to neutral gold, which is a representative system for illustrating some of the challenges of the lanthanides, while having the advantage of having a simpler atomic structure.Although the atom is a nominal one-electron system with a 5d 10 6s ground configuration, gold is a heavy, highly relativistic element that displays competing states similar to the lanthanides.This makes neutral gold a challenging system for theoretical computations (Gaarde et al., 1994).
On the experimental side, several works on the spectrum of gold have been reported.The measurements of energy levels and lines made by Platt & Sawyer (1941) were updated by the works of Ehrhardt & Davis (1971) and Brown & Ginter (1978).The latter two are the references used by the current version of NIST database (Kramida et al., 2022), which will be used in this paper to compare our theoretical results to experimental values.George et al. (1988) also revised the energy values of a few levels.More recently, Civis et al. (2011) measured several high-lying levels and revised a few of the NIST values.The availability of experimental energy levels and transition probabilities is another good reason to use gold to test our theoretical methods.
On the theoretical side, however, calculations for gold in the literature are limited in number, or focused on only a few states or transitions (see Migdalek 2020 and references therein).Safronova & Johnson (2004) present a set of calculations for the 5d 10 nl Rydberg series of Au I using third-order many-body perturbation theory.They obtained good agreement with experimental data as compiled by the NIST database.However, discrepancies were found for some 5d 10 nl states that were attributed to important interactions with the omitted 5d 9 nln ′ l ′ states.More recently, McCann et al. (2022) calculated energies for the lowest 26 levels of Au I using the Grasp 0 (Parpia et al., 1996) and Fac (Gu, 2008) codes, where they obtained an average difference of 8.6% between their Grasp 0 energies and experimental data.
In the following, we present how our targeted optimization calculations with Grasp2018 are done and how we apply them to the case of gold in Sect.2-4.In Sect.5, we discuss the results obtained using our method, by presenting the calculated energy levels as well as transition probabilities for Au I, and comparing them with available Au I data in the literature.Finally, we summarize our findings and discuss possible future work in Sect.6.

Theoretical background
2.1.1.Grasp The theoretical investigations in this work are mainly performed with the Grasp2018 suite of codes (Froese Fischer et al., 2019;Jönsson et al., 2023b,a) making use of the MCDHF and relativistic configuration interaction (RCI) methods.A Grasp eigenstate is represented by an atomic state function (ASF; Ψ) that is defined by an expansion of configuration state functions (CSFs; Φ): where c i are the mixing coefficients, J and M J are the angular quantum moments, P is parity and γ i represents the occupied subshells of the CSF with their complete angular coupling tree information.The Γ label of the ASF is usually taken from the dominant CSF component for the specific eigenstate solution.The CSF basis is constructed through substitutions -singles (S), doubles (D), triples (T) etc.from a set of reference configurations within a specified set of active orbitals.Expanding the CSF basis essentially takes care of electron correlation effects.A multi-reference can also be adopted to include higher-order correlation effects in e.g. a SD model.
In the MCDHF scheme of Grasp, radial orbitals are obtained through a Dirac-Coulomb self-consistent-field (SCF) procedure where each electron moves in the field of the remaining electrons, and the radial functions are varied to minimize a weighted linear combination of a selected set of eigenenergies in the so-called extended optimal level (EOL) scheme.The SCF procedure is continued until the orbitals have converged.Based on this orbital basis set, subsequent RCI calculations may be employed to add relativistic corrections beyond the Dirac-Coulomb approximation, in the form of the Breit interaction and quantum electrodynamical effects, to obtain the final wavefunctions for the targeted states (Jönsson et al., 2023b).
It is well-known that the structure of neutral atoms are challenging to compute accurately, due to the relatively large screening effects among the outer electrons in comparison to the nuclear potential, i.e. large electron correlation effects.Moreover, when doing multiconfigurational calculations, balancing the correlation corrections among the targeted eigenstates is challenging, and the wavefunction representing each individual state may be less accurate than in a single-configuration calculation.Due to these difficulties faced with neutrals and multiconfigurational calculations, convergence issues often occur with the MCDHF SCF procedure in Grasp.To overcome this, we used the DBSR-HF code (Zatsarinny & Froese Fischer, 2016) to obtain better initial estimates of the orbitals.

Fac
To complement the atomic structure investigation carried out with Grasp, we also perform RCI calculations with the Fac code (Gu, 2008).While there are many similarities between the two approaches, e.g. in how the atomic eigenfunctions are represented in terms of j jcoupled CSF's (Eq. 1) and the use of a Dirac-Coulomb Hamiltonian at the core of both codes, a fundamental difference between the two implementations is the representation of the electron-electron interactions and how the set of radial orbitals is optimized.While Grasp is using a MCDHF-EOL scheme to optimize the orbitals, Fac instead employs the computationally more efficient and more practical Dirac-Fock-Slater method.In this approach, a common, local central potential is optimized on an average configuration that is explicitly specified or determined from a set of selected configurations.This means that, in Fac, all orbitals are defined by a common central potential, while in Grasp, the MCDHF approach allows for a larger degree of freedom with the sacrifice of a more complex computational procedure.As mentioned in the introduction, Grasp also makes use of biorthogonal transformation techniques to allow for separate treatment of the even and odd parity states, i.e. including some level of non-orthogonality among the orbitals, while Fac is restricted to use a common orbital set for all states.See the Fac git repository for further information and documentation of the code §.

Spectroscopic CI models
We compare the atomic structure results obtained with the explicit targeted optimization method using Grasp (described in Sect.4) to more conventional spectroscopic CI models (SCI) with Grasp (MCDHF + RCI) and Fac (RCI), i.e. models with a CI space limited to include only those CSF's that represent physical states.SCI atomic structure models are of particular current interest as they are the predominant foundation in the atomic data used for computations of atomic processes and, ultimately, multi-element modelling of r-process dominated kilonovae plasmas (e.g.Kasen et al., 2017;Fontes et al., 2020;Tanaka et al., 2020;Pognan et al., 2023).
In all cases considered in this work, we focus on the first four even and two odd excited configurations (see Table 3).In the case of Grasp, the orbitals are optimized on the even and odd states in two separate MCDHF calculations.For Fac, typical calculations use a central potential optimized on the ground configuration as a starting point.For the present case of Au I, however, the central potential was instead optimized on the average configuration formed by 5d 10 6s and 5d 9 6s 2 for a better energy level structure in comparison to experimental energy levels.

The computational challenges for gold
As mentioned in the introduction, non-orthonormalities between the orbitals can be an important aspect of the computation of heavy r-process elements, such as the lanthanides.The electronic valence shell structure of neutral lanthanides can be summarised as 4 f q 1 5d q 2 6s q 3 6p q 4 , where the q's represent the electron occupation of each configuration.Depending on the subshell occupation, the shape of the 4f, 5d, 6s, and 6p orbitals may vary significantly.Therefore, in a small-scale spectroscopic calculation of a wide range of energy eigenstates of a given lanthanide system, say, it would often be a much too crude approximation to assume that e.g. a single 4f orbital is representative of all the 4f orbitals in the various configurations, and a non-orthonormal orbital set should ideally be used.
Moreover, the states belonging to different competing configurations are sometimes nearly degenerate and thus potentially highly mixed.The mixing coefficients for the configuration interactions and the relative energy positions of the different interacting configuration groups are therefore sensitive to the representation of their respective orbitals and can be challenging to compute correctly with an orthonormal orbital set (e.g.Carlsson, 1988).
Their complex shell structure and large number of levels make the lanthanides a difficult test-bench for model development.A simpler case with similar challenges in terms of non-orthonormalities between orbitals can be found in the copper-group elements.These elements do not have open f-shells and thus have considerably fewer levels per configuration.Their atoms have nd 10 n ′ l Rydberg series accompanied by a number of bound states belonging to the configurations nd 9 n ′ s 2 and nd 9 n ′ s n ′ p, with n = 3 to 5 and n'= 4 to 6 across the homologous group.Similarly for the lanthanides, these configurations involve the same s, p, and d orbitals with different electron occupations.Moreover, the nd 9 states can interact and perturb the nd 10 states (Carlsson, 1988;Gaarde et al., 1994).
In our work, we considered neutral gold, which is the third element in the copper-group, as a test case for the explicit targeted optimization method.Gold is an interesting element from an astrophysical point of view, in the sense that its cosmic origin is unknown to this day (e.g.Kobayashi et al., 2020).It lies close to the third rprocess peak (e.g.Roederer et al., 2022), and like most of the lanthanides, it might be produced in NSM's (e.g.Gillanders et al., 2021).
The electronic structure of the gold atom can be seen in the partial LS -term energy-level diagram in Fig. 1, which can help to better understand the possible interactions between the different states.For example, we pointed out above that the 5d 9 6s 2 and 5d 9 6s6p may interact with states in the 5d 10 nl Rydberg series.However, we can see from Fig. 1 that the 5d 9 6s 2 2 D 3/2,5/2 terms, which can only interact with the 5d 10 6d 2 D 3/2,5/2 states, are situated far below those.The levels belonging to the 5d 9 6s6p configuration, on the other hand, interlap with the 5d 10 np and 5d 10 n f series, and since they have the same parity, a strong interaction between these states is expected, perturbing the Rydberg series.
On top of the perturbed series, another computational difficulty is the shape of the radial orbitals.In Fig. 2, it can be seen how the shape of the 5d, 6s and 6p orbitals vary from states with ten 5d-electrons to states with nine 5d-electrons.To further illustrate the differences between them, the expectation values of their radii are represented with a dashed vertical line, and the values are given in Table 2.The 6s and 6p valence orbitals are more extended when optimized for the 5d 10 states, since they feel more screening effects due to the additional tenth 5d electron.This variation in optimal distribution of the radial one-electron functions is most severe for the 6p orbital, optimized for either the 5d 10 6p or 5d 9 6s6p states, as also shown by the difference in the expectation values of the radius of 6p in Table 2.The computed orbitals that are ideal for representing these various configurations are therefore non-orthonormal.
The combination of strong configuration mixing and non-orthonormalities between the orbitals has a big impact on the computed energy levels of Au I.In Table 1 (fifth column) and on Fig. 3, we show the energies resulting from the conventional Grasp SCI model, described in Sect.2.2.We can see that when the orbitals are forced to be orthonormal (in the sense that the different configurations are computed together using a common orbital set), the level ordering between the 5d 10 6p and 5d 9 6s6p states is wrong.
Moreover, to illustrate the sensitivity of the calculated energies to the way that the orbitals are optimized, we also Table 1.The first few energy levels of Au I calculated using three different orbital optimization schemes, and experimental values from the NIST database (in cm −1 )."Opt.5d 9 + 5d 10 " indicates that the orbitals are optimized on both the 5d 9 and the 5d 10 states.The other 2 columns indicate whether they are optimized on the 5d 9 or 5d 10 states.The root-mean-square deviation is calculated by taking the difference between the calculated and experimental energies.

Level index
Conf.
Term J E(Opt.5d 9 +5d 10 ) E(Opt.5d 10 ) E(Opt.5d 9 ) E(Exp.) show energies from the SCI model but where the orbitals are optimized either on only the 5d 9 or only on the 5d 10 states separately.The results are shown in Table 1 and Fig. 3.We observe that when optimized on the 5d 10 states, the levels belonging to the 5d 10 nl configurations get pushed down in energy, thus the 5d 9 states appear higher in the energy level structure as compared to experiment.The opposite occurs when we optimize them on the 5d 9 states.The MCDHF-SCI calculation with orbitals optimized on all states (5d 9 +5d 10 ) seems to be an in-between situation among the other two extremes, but still gives the wrong energy ordering compared to the experimental energy levels.It should be pointed out that separate Dirac-Focktype calculations with Grasp, for each configuration, would not be optimal either, since this would omit important configuration interactions.One could improve the SCI energies with a large-scale correlation model.The most dominant correlation con- tributors to the calculated energies are often the configurations belonging to the Layzer complex (Froese Fischer, 1997).By adding these to the reference set of configurations, a multi-reference is created that captures interactions arising from single and double substitutions, as well as the dominant higher-order interactions from triple and quadruple substitutions in the correlation model.However, for systems such as Au I, this implies a multi-reference with configurations such as 5d 8 5 f 2 nl and 5d 7 5 f 2 nl 2 that represent the important d 2 → f 2 substitution or substitutions to 5g (although in practice these are likely less important).Large-scale multiconfigurational calculations with such a multireference would quickly become computationally very expensive, especially when the final wavefunctions are to be used to compute other properties, e.g.radiative and electron-impact rates, or opacity-oriented calculations across a wide range of elements and their ions.
In the next section, we address the issue of finding a balance between computational cost, achieved accuracy and applicability, through an explicit targeted optimization approach.

Explicit targeted optimization
We showed in the previous section how Au I demonstrates clear non-orthonormal effects in the orbital basis.We also remind ourselves from the introduction that Grasp and most other current, state-of-the-art atomic structure codes are limited to using orthonormal orbital basis sets.Given this constraint, to correct for the differences between the optimal shape of some of the orbitals belonging to the 5d 10 nl and 5d 9 nl n ′ l ′ states in Au I, and without expanding the model to a large-scale correlation calculation, we investigate a socalled explicit targeted optimization method approach.We focus our attention on low-energy configurations of Au I.
To construct a TO model for Au I, we start by noting that the 5d, 6s, and 6p orbitals show nonorthonormal effects.To remedy this while retaining a compact correlation model, we add a few carefully chosen configurations to the CI space, built from a few new orbitals that are specifically tailored to address the nonorthonormality issues.
In more detail, we start from orbitals generated in a 5d 9 6s nl MCDHF calculation (with nl being 6s or 6p for the even and odd parity calculations, respectively).Then, to represent the remaining 5d 10 nl states (which, as we learned in the previous section, ideally should be described by more extended 5d orbitals), we allow for some radial variation in the subshell by expanding the original configuration as 5d 10 nl = 5d 9 (5d + n ′ d)nl = 5d 10 nl, 5d 9 nl n ′ d and adding these two configurations into the configuration basis, where the n ′ d orbital is a new, "correction orbital" with the same symmetry as the original 5d.For the present range of configurations included in the SCI model, the next available orbital to use for this purpose in the even calculation is 7d, and 6d for the odd.
Furthermore, to correct also for the larger non-orthonormalities present in 6s and 6p orbitals, we expand the configurations as 5d 9 (5d + n ′ d)(nl + n ′ l), where n ′ l are correction orbitals with s or p symmetries for the even and odd calculation, respectively.Expanding this final expression results in a set of new configurations, see Table 3, corresponding to a slightly enlarged CSF space.Using this basis, the final important step is then to optimize the correction orbitals on the 5d 10 nl states only, to force them to correct for the non-orthonormalities with the orbitals from the initial 5d 9 6s nl calculation.Since the correction orbitals are optimised within the active set containing all the other orbitals, the orthonormality is automatically ensured during the MCDHF optimisation procedure through the use of Lagrange multipliers, so that an orthonormal set of orbitals is obtained both for the even and odd parities separately.
The procedure can be summarized as follows: (i) For each parity, identify non-orthonormalities and group configurations accordingly.
(vi) Perform final RCI calculation to obtain all ASF eigenstates.
The scheme outlined here for Au I is readily generalized to any atomic system with similar challenges.After identifying non-orthonormalites, the configurations are grouped with those that have approximately orthonormal ideal orbitals.An important aspect to pay attention to in the construction of a TO model, is which configuration (or configuration group) to start from.In the case of Au I, we started from the 5d 9 states, but one could also perform step (ii) on the 5d 10 nl group.This, however, has the disadvantage of producing correction configurations with 5d 8 subshells (due to the expansion of 5d 9 6s nl, see step (iii)), resulting in a substantially more complex CSF space.Moreover, this will potentially cause a larger CSF space if a correlation calculation based on S and D substitutions from the correction configurations is adopted, due to the presence of configurations with 5d 6 subshells.
The difference between the TO approach and a more traditional MCDHF model lies in the initial observation of present orbital non-orthonormalities, and based on this, how the CSF basis is constructed and the configuration-specific orbital optimization scheme.This procedure of explicitly adding configurations follows what is sometimes called a "configuration-driven approach", in contrast to the more commonly applied "orbital-driven approach" (e.g.Jönsson et al., 2023b), where the CSF basis is constructed through substitutions within an active set of orbitals.Each column represents a different optimization scheme:"Opt.5d 9 + 5d 10 " is a spectroscopic calculation where the orbitals are optimized for all states, while the "Opt.5d 9 or 5d 10 " calculations are done by optimizing the orbitals either on the 5d 9 states or 5d 10 states.Finally, "TO" represents the results obtained with the targeted optimization method.
Another important distinction of the TO method is how we selectively optimize the orbitals on the states of different configurations to coerce them to correct for the non-orthonormalities that we know are present.This is in contrast to conventional approaches where the orbitals usually are optimized on all states.

TO versus other optimization schemes
The energy values resulting from the explicit targeted optimization method are given in the sixth column of Table 4. Compared to the more conventional MCDHF-RCI calculations of SCI type with the orbitals optimized together on all the 5d 10 and 5d 9 states (fifth column of Table 1), we can see that the TO optimization corrects the wrong relative position of interacting states belonging to 5d 10 6p and 5d 9 6s6p, and improves the separation between the 5d 10 6p 2 P 1/2 and 2 P 3/2 levels.Moreover, we observe a significant improvement in the even parity 5d 10 nl states.The relative energy differences ((E exp.− E calc.)/E exp. ) drops from 25.5% to 9.7% for the Table 4.The first few energy levels of Au I calculated using the explicit targeted optimization method (with TO) and without the TO method, as well as with Fac, literature values from McCann et al. (2022), and experimental values from the NIST database (in cm −1 ).The last two lines are theoretically predicted energy levels for which no experimental could be found in the literature.The root-mean-square deviation is calculated by taking the difference between the calculated and experimental energies.
Table 5. Theoretical LS coupling composition of the states of Au I calculated with the explicit targeted optimization method, with level indices in the same ordering as in Table 1 and Table 4.For clarity, we show only the three biggest contributors to the wavefunction.

Level index State
LS -composition 1 5d 10 6s 2 S 1/2 0.95 + 0.02 5d 9 7d 2 S 2 5d 9 6s 2 2 D 5/2 1.00 3 5d 9 6s 2 2 D 3/2 1.00 4 5d 10 6p 2 P 0.78 + 0.15 5d 9 6s 6p 2 F • + 0.04 5d 9 6s 6p 4 F • 5d 10 7s 2 S 1/2 level and from ∼ 25% to ∼ 10% for the two 5d 10 6d levels.This shows the impact of the correction configurations on the 5d 10 series.Overall, the energy results show that, when compared to the other orbital optimization schemes (i.e.optimized for 5d 10 + 5d 9 states simultaneously, or for 5d 10 or 5d 9 ; see Sect.3), the TO method gives the best agreement with the experimental energies, as can also be seen in Fig. 3. Nevertheless, despite the overall improvement achieved by the TO method, we note that the first two excited levels 5d 9 6s 2 2 D 5/2 and 2 D 3/2 appear too high in energy compared to the experimental values and has large percentage differences with 44.7% and 15.7%, respectively.Remaining inaccuracies in the TO excitation energies are likely due to uneven energy contributions from the correction CSF's.For instance, since the spin-angular symmetry of the 5d 10 8s correction configuration is equal to the 5d 10 6s ground configuration, the eigenstates of the latter may be lowered more than the states with other symmetries, such as the ones belonging to 5d 9 6s 2 .Problems like this can be remedied by expanding the small-scale TO model and converging correlation contributions.
We also note that the calculated 5d 9 6s6p energy levels are far too extended in energy, i.e. the configuration is too wide in energy compared to experiment.A possible explanation of such situations can be found in that low-spin states of a given configuration typically have higher excitation energies, and that these tend to require more electron correlation corrections than the lower-lying high-spin states, for which the electrons are forced to spatially avoid each other due to the Pauli principle.This often causes small-scale atomic structure models, i.e. models where correlation contributions are not converged, to predict configurations with too wide-spread excitation energy structures.
The improvement in the energy level structure from applying the TO method could potentially be attributed to the electron correlation contributions when we explicitly include correction configurations and thus slightly expand the CSF basis (Eq.1).To exclude the effect of correlation in the comparison, we also compare our results with a calculation performed within the exact same CSF basis (i.e. the same spectroscopic and correction configurations), but without any targeted optimization in obtaining the orbitals.Instead, we first calculate the 5d 10 and 5d 9 states simultaneously, and then add the correction configurations and optimize the correction orbitals for the spectroscopic states; just as one would do for correlation orbitals.The results from this calculation, labelled "Without TO", are shown in the fifth column of Table 4. Comparing these results with the conventional SCI model "Opt.5d 9 +5d 10 ", we can see that the impact of the competition between the 5d 10 6p and 5d 9 6s6p states is less severe than in the SCI model, but the relative position of the levels is still improved significantly with TO.The representation of the even 5d 10 nl states is also better when we do a targeted optimization, as well as the general agreement with the experimental values.This indicates once again that the optimization scheme of the orbitals has a significant impact on the calculated energies, and how a careful hands-on approach such as TO may improve the results significantly in small-scale calculations.
To further display the effect of the TO method on the energy structure as a whole, we also estimated the rootmean-square deviation for the energies obtained with the different optimizations schemes (see bottom of Table 1 and  Table 4).We find that the deviation from the experimental energies is the smallest for the TO method compared to the other methods.
However, it should also be noted that Au I represent an atomic system that is neither well-represented in either jjor LSJ-coupling.This is particularly true for the 5d 9 6s6p states, as can be seen in Table 5.The strong mixing for some states complicates their labelling and therefore also the level-by-level matching with the experimental levels.Detailed analysis of the dominant LS terms has been done to identify the levels, but some possible differences between the theoretical and experimental dominant LS terms have been reported in the captions of Table 1 and Table 4.These differences could be due to the theoretical model we use and the representation of the wavefunction, or due to the labelling of the experimental levels.Nevertheless, one should be cautious when considering the differences between the theoretical and experimental energies, since this could in some cases be caused by the high mixing in the 5d 9 6s6p states.
5.1.2.Comparison of TO to Fac SCI and other works To further test the TO results, we also compare to Fac SCI calculations (described in Sect.2.2).The resulting energies are reported in Table 4.The obtained excitation structure is comparable to the TO results in terms of relative energy differences to experimental results (an absolute average of 11.5% with TO versus 13.7% with Fac).However, the relative energy positions of the 5d 10 6p 2 P 3/2 and 5d 9 6s6p 4 P 5/2 states remain an issue with the Fac SCI model.
Interestingly, the Fac SCI model gives energies that are in better agreement with experimental values than energies obtained with the corresponding Grasp SCI "Opt.5d 9 + 5d 10 " model (with an absolute average difference of 16.2%).One should note that the Fac model was tuned to obtain the optimal central potential yielding the most accurate energy level structure, which partly explains the good agreement with TO results and experiments.In contrast, the TO model relies only on observing the nonorthonormalities through separate smaller calculations on the various configurations, and then taking these differences into account via state-selective optimization of correction orbitals.
We also compared our calculated energies with the most recent theoretical calculation of these levels in the literature, done by McCann et al. (2022).The authors reported energy levels for the 5d 10 6s, 5d 9 6s 2 , 5d 10 6d, 5d 10 6p, and 5d 9 6s6p configurations, calculated using a modified version of the Grasp 0 code (Parpia et al., 1996).They included correlation by explicitly adding important even and odd parity configurations, resulting in a multireference of a total of 14 configurations, giving rise to 2202 CSFs.Their choice of configurations was done with the motivation to achieve a comparatively good representation of the system while keeping the calculation compact enough for a subsequent R-matrix electron-impact calculation.
It is important to point out that the size of the calculation, i.e. the number of CSFs available to the wavefunction representation, is an important parameter to take into account.Indeed, a larger CSF expansion typically captures more correlation contributions.However, we need to find a balance between the accurate representation of the wavefunction and the size of the calculation.This is particularly difficult for heavy systems such as gold, where the number of CSFs can quickly grow to unmanageable proportions, as was discussed in the end of Sect. 3. Considering that, while the Grasp 0 calculation of McCann et al. (2022) is relatively compact for being a modern Grasp calculation, the TO model achieves a comparable accuracy in the energy level structure with an order of magnitude smaller basis size (with 146 CSFs in total).The TO model also has the benefit of being a systematic and fully ab-initio approach, based only on theoretical observations.Moreover, we can see the issue of competing levels belonging to 5d 10 6p and 5d 9 6s6p in their data set, which is successfully addressed with the TO method.

Transition properties
Despite several theoretical and experimental works on the transitions of Au I (Jannitti et al., 1979;Hannaford et al., 1981;Gaarde et al., 1994;Migdalek & Garmulewicz, 2000;Safronova & Johnson, 2004;Fivet et al., 2006;Zhang et al., 2018)), only 20 lines out of the 191 lines listed in the NIST database have an available transition rate to date.The quasi-totality of these lines lies in the ultraviolet and visible regions, and only very few lines have been observed in astrophysical objects.In the Sun, the most useful Au I line is the 312.28 nm (wavelength in air) between 5d 9 6s 2 2 D 5/2 and 5d 10 6p 2 P o 3/2 (Grevesse et al., 2015), with a precise transition rate measured by Hannaford et al. (1981) using laser-induced excitation.In metal-poor stars, the 267.59 nm line (5d 10 6s 2 S 1/2 and 5d 10 6p 2 P o 1/2 ) can be used as an abundance indicator (Cowan et al., 2002;Sneden et al., 2003;Roederer et al., 2012;Placco et al., 2015), while the 242.79 nm (5d 10 6s 2 S 1/2 and 5d 10 6p 2 P o 3/2 ) can sometimes also be detected (Roederer et al., 2022).Gillanders et al. (2021) have recently searched for a number of Au I lines in the AT2017gfo kilonova, placing upper limits on the abundance of gold in the ejecta.
In order to evaluate the impact of the TO method on the transition rates, we computed E1 transitions for Au I based on the three calculation models: TO, without TO, and MCDHF SCI with both 5d 9 and 5d 10 states.The resulting A-values and experimental A-values from NIST are given in Table 6.We selected transitions with A-values higher than 1.0 × 10 7 s −1 , since weak transitions are more challenging to compute.We also selected the ones that have a NIST accuracy classification of at least lower than 10%.This selection contains three stellar lines mentioned above, as well as one of the lines considered by Gillanders et al. (2021), namely the 751.07 nm (air).We only present the Avalues with the Babushkin gauge (or length form), which is the preferred gauge for transitions involving lowly excited states (Hibbert, 1974), while there is an indication that the Coulomb gauge is more suitable for transitions with highlying Rydberg states (Papoulia et al., 2019).
Table 6 also shows the difference in percentage between the calculated and experimental transition rates, to allow for a better comparison between the different optimization schemes.We once again observe that the optimization scheme of the orbital basis have a significant impact on the calculated properties, like the transition rates.When comparing to the experimental A-values, the TO method gives the best agreement among the three tested calculation models, improving the transition rates by up to 20% for some transitions.Notably, the experimental A-values of the two strongest transitions (at 242.8 and 267.6 nm), which are also the diagnostics for gold in stellar atmospheres, agree within the experimental uncertainty with the values from TO, whereas the other optimization schemes do not.The A-values of transitions involving the 5d 9 6s 2 2 D 3/2 and 2 D 5/2 states have the poorest agreement with NIST values, probably due to the important difference between the calculated and experimental energies.

Fine-tuning of energies
Theoretical energy levels and transition properties can potentially be further improved without increasing the size of the CSF basis by semi-empirically fine-tuning the diagonal elements of the Hamiltonian matrix to observed excitation energies.Adjustment of ab-initio diagonal energy parameters to experimental energies has been used extensively in CI calculations by Hibbert (e.g. 2003), where it was demonstrated that the fine-tuning process gives a better description of the mixing between different LS terms with a common J-value than ab-initio calculations.It has also been shown that fine-tuning can circumvent the problem of wrong relative positions of interacting states and perturbers within Rydberg series (Carlsson, 1988;Lundberg et al., 2001).
Fine-tuning has so far been applied to LSJ-coupled multiconfigurational Hartree-Fock (MCHF) and configuration interaction (CI) calculations.However, it was not applied to jj-coupled multiconfigurational Dirac-Hartree-Fock (MCDHF) or relativistic configuration interaction (RCI) calculations due to the typical presence of large off-diagonal Hamiltonian matrix elements.In their work, Li et al. (2023) investigated a method that can overcome this problem.In their proposed method, the Hamiltonian matrix is transformed from jj to LSJ coupling, applying fine-tuning in this representation, and then transforming back to jj-coupling again, which ultimately modifies both the diagonal and off-diagonal elements.This method has been implemented in the development version of Grasp through the new programs jj2lsj 2022 and rfinetune.
We applied this fine-tuning procedure to the Au I energy matrix obtained with the TO method and we obtained mitigated results.For the even parity levels, the computed energies rapidly converged to experimental values as the jj-LSJ-jj fine-tuning method was applied iteratively.However, this was not the case for the odd-parity levels, which diverged more and more at each iteration of the fine-tuning process.This could be explained by the fact that the 5d 9 6s6p odd parity configuration is subject to large off-diagonal matrix elements in either the jj or LSJ representation, and the best theoretical diagonal energies are therefore far from the experimental ones.A possible way forward could be to apply a non-linear fit of the eigenvalues by varying the diagonal matrix elements.

Conclusions
In this paper, we investigated the application of a method within the realm of MCDHF calculations using the Grasp2018 atomic structure code to improve theoretical excitation energies and transition rates through a small-scale calculation using a compact basis set, such as those employed in multi-element astrophysical opacity calculations.This was achieved by a targeted orbital optimization approach focusing on the perturbing configurations of neutral gold, in order to address the issue of non-orthonormality between the orbitals.
The efficacy of the TO method was assessed by comparing computed energies and transition rates with those obtained through other, more conventional, optimization schemes, as well as with different codes.We compared the results with calculations employing the same CSF basis definition but without targeted optimization of the orbitals, thus excluding the influence of added electron correlation.We found that, in comparison to the other optimization schemes we tried, the energies and transition rates obtained with TO showed the best agreement with the experimental values.The improvement in energies was particularly pronounced for the 5d 10 nl states.Additionally, the TO method effectively rectified the issue of incorrect relative energy positions in the odd states.We also carried out small-scale calculations with the Fac code and they agree well with Grasp when targeted optimization is used, but with TO still yielding a better energy level ordering than Fac.
Certain challenges persist, such as achieving a balance between the optimization of the 5d 10 nl and 5d 9 nln ′ l ′ series simultaneously, especially concerning the even states.We found that while the energy of the 5d 10 nl states improved, the 5d 9 6s 2 states were too high relative to the experimental values.Furthermore, for systems such as Au I that show high mixing between states, in both jj or LSJ representation, the labelling and identification of some levels pose difficulties, thereby complicating comparisons with experimental data as well as other theoretical works.
This work highlights the impact of the choice of orbital optimization scheme and the importance of incorporating non-orthonormal calculations in small-scale calculations that do not rely on a converged energy structure through a systematically expanded CSF basis.While codes such as the (Dirac) B-spline R-matrix, codes by Zatsarinny ((D)BSR; Zatsarinny, 2006;Zatsarinny & Bartschat, 2008) allow for fully non-orthonormal structure calculations, they are fundamentally scattering codes and computationally intractable for calculations of complete, spectrum calculations of heavy elements.Implementations of fully non-orthonormal atomic structure methods pose many numerical challenges, making the TO method an effective compromise.
Looking ahead, the TO method is promising for application to more complex systems, such as the lanthanide group of open f-shell elements.Consider, for instance, Gd I, which shares a similar structure to the copper-group elements that Au is part of, featuring configurations with seven and eight f-electrons and degenerate energy levels from different Rydberg series.We predict significant non-orthonormality effects in for example the 4f and 6s orbitals from the 4 f 7 and 4 f 8 series, which could be addressed through targeted optimization to obtain a more accurate average orbital for these two series.Furthermore, this method can be seen as a correction of small-scale structure models in e.g.multi-element efforts in astrophysical applications, or as a basis for large-scale calculations of single ions targeting high-accuracy atomic data.

Figure 2 .
Figure 2. The electron density representing the shape of the 5d, 6s and 6p orbitals optimized for different configurations, based on separate Dirac-Fock calculations.The dashed vertical lines represent the ⟨r⟩ expectation values.

Figure 3 .
Figure 3.Comparison between the calculated energy levels of Au I (first four columns) and the experimental energies from NIST (last column).Each column represents a different optimization scheme:"Opt.5d 9 + 5d 10 " is a spectroscopic calculation where the orbitals are optimized for all states, while the "Opt.5d 9 or 5d 10 " calculations are done by optimizing the orbitals either on the 5d 9 states or 5d 10 states.Finally, "TO" represents the results obtained with the targeted optimization method.
In our calculations, this term is highly mixed with the 4 F 3/2 term which is the dominant LSJ term b In our calculations, this term is highly mixed with the 4 P 1/2 term which is the dominant LSJ term c In our calculations, this term is highly mixed with the 2 P 3/2 term which is the dominant LSJ term d In our calculations, this term is highly mixed with the 4 F 5/2 term which is the dominant LSJ term Civis et al. (2011). (2011)

Table 2 .
Partial energy-level diagram for neutral gold.The black lines represent the even parity states while the blue ones represent the odd parity levels.The numbers next to the levels of Rydberg series represent the principal quantum number n.For the 5d 9 6s6p states, the fractions are the values of J.The energy values are experimental values from the NIST database.Expectation values of the radius ⟨r⟩, in atomic units, of the 5d − , 5d, 6s, 6p and 6p − orbitals optimized on the states belonging to the configurations in the first row, based on separate Dirac-Fock calculations.

Table 3 .
The list of the spectroscopic configurations and the correction configurations used in the targeted optimization calculation.

Table 6 .
Transition rates (s −1 ) calculated with different models, compared to experimental values from NIST, for selected E1 transitions of Au I.