Method of accounting for higher harmonics in the calculation of the electronic structure of the states of a linear molecule excited by high-intensity X-rays

Modern development of high-intensity and high-resolution X-ray technology allows detailed studies of the multiphoton absorption and scattering of X-ray photons by deep and subvalent shells of molecular systems on a wide energy range. The interpretation of experimental data requires the improvement of computational methods for obtaining excited and ionized electron states of molecular systems with one or several vacancies. The specificity of solving these problems requires the use of molecular orbitals obtained in one-center representation. Slow convergence of one-center expansions is a significant disadvantage of this approach; it affects the accuracy of the calculation of spectroscopic quantities. We offer a method of including higher harmonics in onecenter expansion of a molecular orbital with the use of wave functions of electrons of deep shells of a ligand (off-center atom of a molecule). The method allows one to take into account correctly electron density of a linear molecule near the ligand when describing vacancies created in a molecular core leading to radial relaxation of electron density. The analysis of the parameters of one-center expansion of the ligand functions depending on ligand’s charge is performed. The criteria for the inclusion of higher harmonics of one-center decomposition of the ligand functions in the molecular orbital are determined. The efficiency of the method is demonstrated by the example of diatomic molecules HF, LiF, and BF by estimating energy characteristics of their ground and ionized states.


Introduction
The advent of compact and tunable X-ray sources (X-ray free electron lasers (XFEL) [1][2][3], LWFA-based inverse Compton sources (LWFA) [4,5], Laser-plasma accelerators (LPA) [6], Highharmonic generation sources (HHG) [7]) opened up new opportunities for the study of the electronic structure of various many-electron systems, including simple molecules and biomolecules, nanoparticles and nanocrystals [8].It becomes possible to study in detail the processes of inelastic scattering, multiple absorption and ionization in gas-phase molecular systems with high resolution.
Studies of the process of ionization / excitation of molecular systems [9][10][11] carried out so far show that the creation of a vacancy in the molecular core can lead to a significant redistribution of the electron charge and has a many-electron character, which has certain peculiarities when the system is ionized by high-intensity X-rays [12][13][14].When a molecule absorbs x-ray photons, innershell multiple ionization induces fragmentation dynamics [15].Chemical bonds are weakened, and electrons and holes rearrange themselves before the molecule breaks apart [16][17][18].In this case, the improvement of calculation methods of highly excited and ionized states of molecule is crucial for understanding the interactions between a molecule and the intense X-ray pulses.
One of the most well-known and widely used methods for calculating the energy characteristics of molecular systems is the linear combination of atomic orbitals (MO LCAO) method [19].Within the framework of the present problem, this method has limitations, since it is mainly used to obtain ground and low-excited states of molecules.When solving the problem of interaction of an x-ray photon with a molecular system for the calculation of the observed values, the use of a one-center (OC) method is more effective [20,21].The OC method is based on the representation of the molecular orbital (MO) in the form of expansion over spherical harmonics with respect to a single point in space usually located at the origin of the molecular system.The use of this approach is necessary when calculating the observed spectroscopic characteristics for the processes of absorption and/or scattering of x-ray photons, it allows one to significantly simplify the calculations of the energy characteristics of molecules [22,23].In the OC approach, the mathematical apparatus developed for atomic calculations is easily generalized to the case of molecular calculations for one-and two-particle matrix elements which form the basis for most of the calculated physical characteristics of polyatomic systems [24,25].In the course of its development, the OC method was implemented through solving a secular equation using basic Hartree-Fock functions obtained in the field of the central atom of the molecular system [13,14,26], and through solving a system of coupled integro-differential equations [27,28].The principal disadvantage of the OC representation of MO is the slow convergence of the expansions in spherical harmonics.This leads either to an increase in computational problems when using huge basis sets for solving secular equations, or to the growth of the order of a system of coupled integro-differential equations.In this case, the inevitable limitation of the number of members of the OC expansion leads to a significant problem of decreasing of the accuracy in the calculation of the observed values, especially the energy characteristics.Attempts to solve this problem in [29,30], were based on an increase in the participation rate of the terminal partial harmonic of the OC expansion in order to compensate the lack of MO's electron density caused by the limitation of the OC series.However, such a "technique" does not allow correct calculation of the energy of the Coulomb interaction which expectedly leads to the loss of accuracy in the calculated energy characteristics of the molecular system.
In this paper, we present a theoretical basis for taking into account the higher harmonics of the OC expansion of MO's of highly excited and ionized states of molecular systems ensuring accurate calculation of the highly-coupled electron density of molecular ligands in polyatomic systems.The method is universal; it can be used both to solve secular equations with the construction of MO on certain basis set, and to solve a system of coupled integro-differential equations.A computational routine for obtaining ligand wave functions is presented.The characteristics of the ligand wave functions are analyzed with respect to the charge of its core, and the criteria for their inclusion in the MO are formulated.The efficiency of the developed method is demonstrated by the example of calculating the energy characteristics of a number of diatomic molecules: HF, LiF, and BF.

Method of calculating molecular orbitals taking into account higher harmonics
In the one-center approach, a single-particle wave function describing the behavior of an electron of a given MO is a expansion in spherical harmonics; it has the form [31]: Where    () =   ()  (, ) is l-symmetry partial harmonic of MO;  ≡ {, ,  } are spherical coordinates;  (  ) are quantum numbers characterizing one-electron states of the molecule with the point symmetry group  ∞ .The proposed approach to accounting for higher harmonics in a single-particle wave function is based on the nature of the action of the local ionic potential of the ligands included in the molecular system.With the growths of l of the partial harmonic    which is included in the OC expansion of MO, the effect of the ligand field on the radial part   () of the harmonic increases, and becomes decisive.Starting from a certain l = ζ the topology of the radial part   () of the partial harmonic    becomes unchanging.In this case, a single-particle MO can be defined by two components and have the form: Where Ω  are the participation factor of the wave function    of the L-th ligand of the molecule.
The first sum in (2) includes partial harmonics with orbital moments  ≤  , which are the solutions of one of the common methods for obtaining MOs in the OC approach for molecular systems [26][27][28].The behavior of these partial harmonics is largely determined by the molecular core of the system and the electron-electron interaction characterizing the chemical bonding in the system.The second sum in ( 2) is defined as a linear combination of the so-called ligand function.The functions    describe the electron density of deep electron shells strongly bound to the ligand of the molecule.In    , only higher harmonics of the OC expansion (l>ζ) are included.They represent OC expansions of the wave functions of electrons in the shells of the off-center ligand atom, and have the form: To determine the participation coefficients Ω   of the ligand function    in one-electron MO, the condition of equality of the radial parts of MO   () of the partial harmonic with l = ζ from the first and the second sum in ( 2) is used: The solution of the system of linear algebraic equations with respect to unknown participation coefficients Ω   based on relation (4) allows, considering the normalization procedure (2), to determine the desired function Ψ  for each MO of the molecular system.

Numerical procedure for obtaining one-center wave functions of ligands
The calculation of the functions    () is based on the possibility of representing the wave functions of electrons of shells of an off-center ligand atom in the form of a series in spherical harmonics relative to a selected center of the molecular coordinate system.In this case, the wave function of the ligand atom's electrons has the form [32]: Here r'≡ {r ', ϑ', φ '} are spherical coordinates in the atomic system associated with the molecular ligand, r≡ {r, ϑ, φ} are spherical coordinates in the molecular coordinate system, which are interconnected by the relation

Results and Discussion
As objects for testing the method of accounting for higher harmonics of OC decomposition of MO, diatomic molecules were chosen: hydrogen fluoride (HF), lithium monofluoride (LiF) и boron monofluoride (BF).The ligand nucleus charge changing with a constant step allows estimating ligand's effect on the parameters of the ligand functions, as well as determining the conditions for the inclusion of higher harmonics of the OC expansion in the MO of electronic states using the example of selected molecules.In addition, these molecules have relatively small number of electron shells; this allows us to find out easily the possibilities of the method developed in the work when calculating the energy characteristics of their ground and ionized states.
During the calculations, the molecular reference frame was centered on the F atom.To systematize one-electron states, the quantum numbers characteristic of the point symmetry group molecules  ∞ were used.When calculating the wave functions, the parameters of the ligands used in the calculation of the higher harmonics of the OC expansions of MO, and the equilibrium distances (R) between the ligand and the fluorine atom, i.e. interatomic distances [33] in the HF, LiF and BF molecules, are presented in Table 1.The calculation of ground and ionized states MOs of the molecules under consideration was performed by solving the system of coupled integro-differential equations in accordance with [27].The choice of this method allows avoiding the problems of selecting the optimal basis sets for calculating the electronic states of each of the molecule under consideration.
The ligand functions were calculated in accordance with the procedure presented in Section 3 of this work.The wave functions of ligand electrons (see Table 1) were obtained beforehand by solving the Hartree-Fock system of integro-differential equations using the software developed at the Department of Physics at Rostov State Transport University [13,34].To improve the accuracy of the energy characteristics of molecular states when calculating MO wave functions, the computational scheme with a constant step along radius r was replaced with the scheme with a variable step along r in accordance with the expression: The parameters δ and θ in (7) determine the logarithmic dependence of the r values of the new radial grid Λ, which provides an increase in the concentration of radial points in the vicinity of the origin, and thus provide necessary accuracy in the region of oscillations of the fluorine atom wave functions.The parameter ω in (7) provides an increase in the density of r points in the region where the ligands are located capturing the region of oscillation of wave functions in the vicinity of the ligand.The OC expansions of atomic wave functions in spherical harmonics with respect to molecular reference frame is performed for orbital moments l up to 100.
Figure 1 presents calculated one-electron energies (OEE) [35] of 1s states (panel a) and 2p states (panel b) of the ligand atoms H, Li, and B. Calculations are performed with different numbers of included OC expansion terms.For placing the curves in Fig. 1 closer together, calculated OEE are normalized to exact OEE obtained by the Hartree-Fock method for isolated ligand atoms on their own centers.For the 1s states (Figure 1a), it can be seen that with an increase in the charge of the ligand, the number of the OC expansion terms necessary for reaching the saturation levels of OEE increases.Indeed, to achieve 99.9% of the exact OEE for the case of the ligand H not more than 30 OC expansion terms are required, Li ligand requires not more than 80 OC terms, and for the B ligand not more than 90 OC expansion terms are needed.At the same time, the transition to non-hydrogen ligands causes a sharp increase (by about 3 times) of required OC expansion terms, then a rather weak further relative growth of this parameter is seen when going from the Li to the B ligand (only by about 11 %).The behavior of the 1s states with an increase in the nucleus charge of the ligand indicates the stabilization of the number of terms required for inclusion in the OC expansion.In this respect the 2p ligand states are different; the required number of the OC expansion terms grows steadily with the increase of ligand's nucleus charge.

Preprints
Energetically, 2p states lie above the 1s states which "shields" the effect of the ion field on them.As can be seen from Figure 1b, the saturation of calculated OEE is observed with a smaller number of the OC expansion terms.Thus, to achieve 99.9% of OEE for the case of the H ligand, not more than 5 OC terms are necessary, Li ligand requires not more than 11 OC terms, and for the B ligand, not more than 20 terms are needed.The observed pattern indicates a decrease in the role of higher harmonics of the SC decomposition and allows including only the ligand functions of s-and psymmetry, when calculating molecular system's ground, excited and ionized states.
Figure 2 shows the radial parts of the partial OC harmonics of the ligands' 1s states: H1s (panel a), Li1s (panel b), B1s (panel c).Shown are the partial harmonics with l up to 100.With the increase of ligands nucleus charge, a change in the topology of the radial parts of the partial OC harmonics is observed, accompanied by an increase in their symmetry.The electron density under the action of the nucleus charge field is concentrated around the nuclei of the ligands, reducing the occupied area of space along the radial direction from 10 a.u.for the H atom to ~ 1.5 a.u.for B (more than 6 times).The growth of the nucleus charge triggers an increase in the number of necessary OC expansion terms having an amplitude  1 () •  ≥ 0,1: from l ≤4 (for H) to l ≤ 25 (for B).Such behavior of the OC expansion spherical harmonics indicates a rapid transition to the profiles symmetrical with respect to the positions of ligands, due to an increase in the nucleus charge of the ligand.These harmonics differ from each other only by their maximum values; this makes it possible to significantly simplify the summations caused by the OC expansions of the ligand wave functions.Symmetric profiles of partial harmonics in the case of larger nuclear charges indicates the possibility of a significant reduction in the need for computing resources.The calculations of the integrals in the matrix elements with higher harmonics can be reduced to the algebraic summation of the components of the OC series, which differ from each other only by numeric coefficients.In our case, it is necessary to determine the asymptotic behavior of the partial harmonics against spherical harmonic orbital momentum quantum number.
Figure 3 shows the dependence of the maximum value of the radial part of the partial harmonic on its l quantum number for the case of ligands 1s states.The maximum values of the radial parts of the partial harmonics are those at r = R, the F-to-ligand equilibrium distances (see Table 1, Figure 2).For the imposition of curves in Figure 3, the normalization is performed to the maximum value of the radial part of the partial harmonic with l = 25, which is maximum in the range of l ≥ 25.The hyperbolically decreasing curves in Figure 3 do not change significantly with increasing the ligand nucleus charge keeping their topology.Drawing the approximating curves allows one to estimate the maximum value of the radial part of partial harmonics at the intervals of the OC expansion terms with l > 100.
The higher harmonics of the OC expansion of ligand functions thus characterized make it possible, when calculating the observed quantities, to introduce an effective and low-cost method of taking them into account and solve the problem of slow convergence by going from integration to algebraic summation for the harmonics with symmetrical profiles differing only by numerical coefficients.
The results of the calculation of total energies of the ground and ionized states of HF and LiF are presented in Table 2. Column (a) of Table 2 lists the total energies calculated by solving the system of coupled integro-differential equations [27] with a limited series of OC expansions, l ≤ 7. Column (b) of Table 2 presents total energies calculated by solving the system of coupled integro-differential equations [27] with accounting for higher harmonics of OC expansions presented in this work.In these calculations the harmonics with l up to 35 were calculated straightforwardly.Column (c) of Table 2 presents the total energies calculated within the MO LCAO approach using GAMESS software package [36].The following abbreviations are used in Table2  ≡ 1 2 2 2 3 2 1 4 (4 2 );  − ≡ 1 1 2 2 3 2 1 4 (4 2 );  − ≡ 1 2 2 1 3 2 1 4 (4 2 );  − ≡ 1 2 2 2 3 1 1 4 (4 2 ).Comparison of calculations ground state (GS) energies of molecules shown in columns (a) and (b) of Table 2 showed that the inclusion of higher harmonics of the OC expansions of ligands wave functions an additional contribution to the total energy proportional to the number of terms in OC expansions.The contribution of higher harmonics to the energy of the electronic configuration of the ground state of the molecule increases with increasing ligand nucleus charge.The changes in energy observed in column (b) of Table 2 are due mainly to the contribution of the 1s-states of the H and Li, the ligands in the HF and LiF molecules.The inclusion of higher harmonics of the OC expansions, all other conditions being equal, provides obtaining of ground state total energies of HF and LiF comparable in accuracy to the results of the MO LCAO method (see  2).When switching to molecules with non-hydrogen ligands, the inclusion in the calculation of energy characteristics of the higher harmonics of the OC expansions ensures accounting for the contribution from electron density strongly related to the molecular core; it significantly improves the results (see column (a), (b) and (c) Table 2 for LiF).The results of calculations for ionized states of 1σ -1 , 2σ -1 and 3σ -1 of HF and LiF are given in columns (a) and (b) of Table 2. Ionization of the deep 1σ-shell leads to an increase in the electron density near the nuclei of molecules due to radial relaxation, and partially compensates for the limited inclusion of the number of terms in OC expansion of MO.At the same time, a consistent increase in the number of included OC expansion terms will provide an increase in the accuracy of the calculated characteristics of molecular systems.Ionization of valence 2σ and 3σ shells of HF and LiF has a markedly smaller effect on the behavior of the electron density of the molecular core.In this connection, an expected small changes in the values of the total energies of the ionized states of the HF and LiF molecules are observed upon taking into account the higher harmonics of the OC expansions of ligands wave functions.3 presents the total energy calculated by the method of a system of coupled integro-differential equations [27] with a limited series of OC expansions with l ≤ 4. In column (b) of Table 3, total energies are calculated taking into account the higher harmonics of the ligands wave functions, up to l = 60.
The growth of the charge of the ligands of the molecular system significantly increases the influence of the higher harmonics of the OC decomposition on its energy characteristics.The data in columns (a) and (b) of Table 3 for these states show that with increasing ligand charge, the influence of higher harmonics of the OC expansions determined by the electron density strongly bound to the nuclei of the molecular system, increases significantly.Thus, the presented method of taking into account the higher harmonics of the OC decomposition in MO of excited and ionized electronic states becomes effective in the studies of molecules with non-hydrogen ligands.

Conclusions
This paper presents a method for taking into account the higher harmonics of the one-center decomposition of ligands functions in the calculation of molecular orbitals which removes the principal limitation of the one-center method, its slow convergence.The effect of the charge of the nucleus of the ligand on the parameters of the higher harmonics of one-center expansions is analyzed.It has been established that the growth of the charge of the ligand allows, in practical calculations, to restrict oneself to the ligand functions of only s-and p-symmetry.Switching to non-hydrogen ligands causes the "stabilization" of the number of one-center expansion terms necessary to achieve required accuracy of the calculation of the energy characteristics of molecular systems.The nature of the asymptotic behavior of partial harmonics of one-center expansions against their orbital momentum quantum number is determined.The results showed that taking into account the higher harmonics of the expansion of the wave functions of ligands electrons of describing the electron density near the nuclei of the molecular system is a necessary condition when performing calculations in the OC approach to obtain accurate energy characteristics of molecular systems excited and ionized by Xrays.

Figure 1 .
Figure 1.Dependence of calculated one-electron energies of the nl states of ligands normalized to respective energy obtained by own-center atomic calculation in the Hartree-Fock approximation as functions of the number of spherical harmonics included in the OC expansions: (a) 1s states of H, Li and B atoms, (b) 2p-states of H, Li and B atoms.

11 December 2018 Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 11 December 2018 doi:10.20944/preprints201812.0121.v1
being the distance from the ligand to the beginning of the molecular coordinate system;   (′) is the wave function of electrons of atomic shell nl of the ligand,    () is the nl electron wave the function in the molecular coordinate system,   is a normalization coefficient.

Table 1 .
Fluorine-ligand distances and ligand wave functions included in calculations

Table 2 .
Total energies of ground and ionized states of HF and LiF (in Ry).

Table 3 .
The total energy of isoelectronic configurations (in Ry).

Table 3
presents calculated total energies for HF and molecular ions of LiF and BF.Calculations are performed for the electronic configuration of the ground state of the HF molecule.Column (a) of Table