Exploring spin symmetry-breaking effects for static field ionization of atoms: Is there an analog to the Coulson–Fischer point in bond dissociation?

Discussions on the usefulness of symmetry-broken approximate solutions are a familiar topic in time-independent quantum chemistry.^1–6 The issue is encapsulated in what Löwdin has called the “symmetry dilemma:”⁷ a more flexible trial wavefunction that does not necessarily preserve all the symmetries of the exact Hamiltonian might lead to better energetics (i.e., lower energy on account of the variational principle) at the cost of losing good quantum numbers that characterize the state of a given system. Within the single determinant Hartree–Fock (HF) model, classification of these symmetry-broken solutions is based on group theory considerations,⁸ but, most commonly, in electronic structure, one uses a terminology that reflects constraints imposed on the orbitals that comprise the single determinant.⁹ For instance, requiring that both α and β spin–orbitals share a common set of spatial functions (i.e., the electrons are paired whenever possible) leads to the well-known restricted closed-shell (M_s = 0) Hartree–Fock (RHF) and the more general (M_s ≠ 0) restricted open-shell HF (ROHF) models. RHF and ROHF are both eigenstates of total spin, $Ŝ2$⁠, and its z component, $Ŝz$⁠. Lifting this spin pairing constraint, such that α or β spin–orbitals can have different spatial functions, gives us the unrestricted Hartree–Fock (UHF) model whose wavefunction is no longer an eigenfunction of $Ŝ2$⁠. Further symmetry lowering by abolishing the notion of separate sets of α and β orbitals leads to the generalized Hartree–Fock (GHF) wavefunction, which is not an eigenstate of either $Ŝ2$ or $Ŝz$⁠. Furthermore, number symmetry can also be relaxed following Hartree–Fock Bogoliubov theory,^10,11 which yields a state that is not an eigenstate of the particle number operator (⁠$N̂$⁠).

Within a finite basis, the HF energy, E_HF(θ), is a function of orbital rotation parameters, θ, that mix occupied and virtual orbitals.^12–14 A solution of the Hartree–Fock equations zeros the orbital rotation gradient such that ∇_θE_HF = 0 (i.e., it is guaranteed to be a stationary point). However, a solution is not necessarily a minimum, and an analysis of the orbital Hessian, $EHFθθ$⁠, is required to characterize the nature of a stationary point. If all the eigenvalues of the Hessian are greater than zero, we have found a solution that is a local minimum and is said to be stable. The question of which set of orbital parameters to include in the evaluation of the Hessian arises:¹⁵ if the HF solution is stable within the manifold determined by certain symmetry constraints (e.g., the RHF constraints), it is said that such a solution is internally stable. Lifting some or all these symmetry constraints leads to characterization of the stationary point in manifolds of higher dimension than the one it was originally optimized in (e.g., characterizing an RHF solution in the space of UHF variations). This corresponds to an analysis of the external stability of the given solution. It might be useful to note that stability analysis is closely connected with time-dependent density functional theory (TDDFT) and time-dependent Hartree-Fock (TDHF) linear response equations,^6,16 which are routinely available in standard quantum chemistry packages.

For common problems in ground state electronic structure, internal stability analysis avoids convergence to spurious saddle points and excited states. On the other hand, external stability analysis sometimes reveals interesting physical insights into the nature of electron correlation. For instance, when single bonds are stretched beyond the so-called Coulson–Fischer (CF) point,¹⁷ the lowest triplet excited state (T₁) starts to mix with the singlet ground state (S₀) in a process known as the “triplet instability,”^18,19 which leads to spin symmetry-breaking. The resulting spin polarized UHF state has lower energy than the RHF state but is no longer an eigenstate of $Ŝ2$⁠: the state is said to be spin-contaminated due to its mixed singlet–triplet character. The RHF solution that pairs electrons in order to keep a well-defined $Ŝ2$ eigenstate fails to provide a qualitatively correct description of the dissociation process. The RHF pairing constraint, which fixes natural occupation numbers at 0 or 2 even at the dissociation limit (when one would expect the two orbitals to be singly occupied), leads to a state that spuriously preserves ionic character, which consequently fails to approach the correct asymptotic limit at complete dissociation. On the other hand, UHF is able to successfully provide a qualitatively accurate description of the ground state potential energy surface for single bond dissociation.¹³ Similar considerations apply to molecules that are singlet diradicaloid^20–22 in character: because two orbitals have occupation numbers significantly different from two and zero, RHF cannot be qualitatively correct, while UHF exhibits spin contamination due to mixing of S₀ and T₁.

There are surely other situations in which similar symmetry dilemma arises, but where its consequences have not been so thoroughly explored as bond dissociation processes and singlet diradicaloids. One example lies in the interaction of atoms and molecules with strong electric fields.^23–31 Strong field chemistry and physics are rapidly developing fields because of the rich array of new highly non-linear phenomena that emerge. For example, a strong oscillating field will nearly ionize bound electrons in one direction in a first half-cycle followed by reattachment and near ionization in the opposite direction in the second half-cycle. This leads to high-harmonic generation (HHG),^32–36 the emission of radiation by the driven bound system at frequencies that are many times that of the applied radiation. HHG and related phenomena are building blocks for the new field of attosecond science.^37–40 Therefore, the use of quantum chemistry methods to model strong field phenomena is also attracting increasing interest.^41–46 However, to date, there has been no systematic exploration of the role of symmetry-breaking in the mean-field HF method in this context, to our knowledge. This work represents a first step in this direction, though the issue of inadequate field-free Hartree–Fock reference for real-time methods in strong field environments has been briefly mentioned in some previous works.^47–50 

By definition, strong fields are those whose scale approaches the strength of internal electric fields experienced by the valence electrons of atoms and molecules (e.g., 1 a.u. for the H atom). Therefore, the molecule and the field must be considered as a combined system, rather than treating the field as a perturbation. If such fields are static, then strong field ionization becomes possible and some electrons may be unbound. Yet because the molecule-field system is treated as a whole, the symmetry dilemma may arise for at least some values of the applied field. In this sense, the magnitude of the field is a control variable similar to the degree of bond-stretching in the dissociation of a closed-shell stable molecule. This paper explores the role of symmetry-breaking in HF solutions for atom-field systems where a static field is scanned across values that can strip one or more electrons from the atom. This investigation is interesting in its own right and may also help to set the stage for a subsequent analysis in the context of time-dependent fields. Our paper is arranged as follows: in Sec. II, we discuss our approach to approximate the description of continuum-like states using ghost basis functions. An analysis of the HF spin symmetry-broken solutions for the static field ionization of helium and neon is presented in Secs. III and IV, respectively, and in Sec. V, we assess our main conclusions while presenting an outlook of future work.

The first issue that arises when applying an electric field to an atomic or molecular system is that, strictly speaking, bound states are no longer supported by the combined potential arising from the system and the field.^51,52 These states are now resonances, and a suitable discretization for both localized and continuum states is needed to describe them. From a theoretical perspective, the scattering community has developed a wide range of tools to treat these resonance states, such as exterior complex scaling (ECS),^53–60 complex absorbing potentials (CAPs),^61–69 and different grids as schemes to discretize space.^70–76 On the other hand, the quantum chemistry community has long advocated for the use of atomic orbital (AO) expansions based on Gaussian functions as an efficient way to compute ground and excited state properties of molecular systems. Since our goal is to study static field ionization exploring the common toolbox of quantum chemistry, we shall resort to the usual finite AO basis set treatment. What are the consequences of such choice? The most concerning one is that ionization losses are hindered since electronic density will be constrained near the atom/molecule. Two possible ways⁷⁷ of circumventing this issue are the use of highly diffuse basis sets and discretizing the relevant part of the space by adding ghost atoms/functions that could be populated by electron density as ionization takes place. We have opted for this second option and added a series of ghost functions in a line passing through the atom and in the direction of the applied field. This allows us to capture some of the effects of ionization as electron density can escape the system and populate the discretized space (Fig. 1). This setup is inspired by previous works that aimed to study the real-time electron dynamics of systems in the context of strong fields and high-harmonic generation (HHG), where the addition of diffuse and ghost functions to the basis set is important for a good description of the Rydberg and unbound states of the system, respectively.⁷⁷ 

Given this setup for our problem, we also need to interpret what the exact results would look like in this finite basis set approach. We start our analysis by noticing that, assuming an adequate discretization of the space by the ghost centers, we can separate the problem/potential into two parts: one associated with the molecular/atomic system and another related to the artificial box that is effectively created by the extent of the ghost centers. Introductory quantum mechanics teaches us that, to first order, the energy levels of a charged particle in a box subjected to a uniform static electric field decreases linearly with the strength of the field.⁷⁸ Moreover, increasing the field strength affects the molecular potential by suppressing the barrier to ionization for higher excited states (especially Rydberg states) and increasing the tunneling probability for low-lying states. In this sense, higher excited states easily become continuum-like states and “dive” down in energy as the field becomes stronger. At certain field strengths, these continuum-like states start to interact with the low-lying bound states and an avoided crossing is observed: the former bound state acquires continuum-like character by localizing around the edge of the artificial box, with an energy that decreases linearly as the field strength continues to increase. Note that the applied electric field turns all the states of the system into resonances with finite lifetime,⁷⁸ and a proper discussion on how to obtain DC Stark lifetimes for several molecular systems using quantum chemistry methods has already been presented elsewhere.^43,44,79

Figure 2 illustrates this process for a He atom (described by the functions contained in the aug-cc-pVTZ basis set^80,81) in an artificial box comprised of 18 hydrogen STO-3G⁸² s-type ghost centers spaced by 0.5 Å spanning a range of 4.5 $Å$ around the atom (which should be adequate for an unambiguous determination of the state of ionization of the system, given that this distance is much larger than the atomic radius of He). We shall focus on the behavior of the ground state S₀ and the first singlet excited state S₁ as a function of field strength. Initially, in the absence of the field, these states are well separated in energy and there is no mixing between them. As the field strength increases, we observe the characteristic quadratic polarization effect on the energy (which can be described by perturbation theory).⁸³ Higher lying states are more polarizable, so this effect becomes more evident. These high energy excited states are also more easily ionized, and the transition between the quadratic behavior characteristics of a bound state to a linear dependence of the energy on field strength happens at lower fields. At F = 0.175 a.u., we observe an avoided crossing between the first singlet state (which has been previously ionized) and the unionized ground state. This leads to a change in character of the ground state after the avoided crossing: initially unionized, with both electrons localized around the atom, the ground state becomes singly ionized, with one electron detaching from the atom and localizing at the edge of the box. The second ionization is characterized by yet another avoided crossing between the singly ionized S₀ and doubly ionized S₁ at a field of F = 0.36 a.u. Between these two field strengths, we note that the singly ionized T₁ state is degenerate with the ionized ground state S₀. It is worth noting that the potential energy curves presented here as a function of field strength are very similar to stabilization graphs proposed by Simons⁸⁴ and widely used to describe temporary ions and other atomic/molecular resonances.⁸⁵ There, information about the resonance’s lifetime can be obtained by an analysis of the avoided crossing that arises by changing the variational parameters associated with the basis set. Here, however, the basis set is fixed and the energy is plotted as a function of varying field strength.

Our goal is to understand how different approximations, such as HF and MP2, can recover this exact finite basis behavior. We shall put emphasis on the effects of allowing for spin polarization (i.e., difference between using restricted and unrestricted orbitals) to describe field ionization as an analog to what happens in bond dissociation of molecules. All calculations were performed with the Q-Chem 5.3 package⁸⁶ following the same protocol: we used the aug-cc-pVTZ basis set for the atomic center and hydrogenic STO-3G s functions for the ghost centers. We also performed internal stability analysis in order to ensure that the solution within the domain of each HF class (RHF, UHF, and GHF) was the lowest possible.

We start our discussion about the relation between spin polarization and static field ionization by analyzing the behavior of the Hartree–Fock (HF) solution in a simple, yet instructive, minimum basis model of a two electron system. This model is comprised of a single basis function localized around the atom (represented by $A$ in Table I) and a single function (represented as $C$ in Table I) placed at a distance d from the atom in the direction of the applied field to represent the discretized continuum and that can support the flux of ionized electrons once the field is switched on. It should be noted that this is analogous to the toy model used to describe the bond-stretching of H₂ in the minimum basis.^6,13 Here, we have also assumed that the distance d is large enough such that there is no overlap between the orbitals representing the atom and the discretized continuum. Moreover, we assumed that the ghost function is diffuse enough that we can neglect its kinetic energy (⁠$〈C|T̂|C〉≈0$⁠), but not too diffuse relative to the distance d. We can then construct the HF orbitals as linear combinations of these two functions.

In this toy model, the RHF solution is controlled by a single variational parameter θ that mixes the atomic basis function with the continuum-like ghost function (Table I). Analyzing the stationary and stability conditions of the energy obtained as the expectation value of the Hamiltonian with the RHF determinant leads to two limiting cases. The solution characterized by the double occupancy of the atom-centered orbital (θ = 0) is found to be stable when the product of the applied field strength (F) and the distance between the basis functions (d) is less than IP₁, the first ionization potential (IP) of the system (i.e., Fd < IP₁). On the other hand, the $θ=π2$ solution, which corresponds to two electrons paired up on the ghost function that supports the ionized flux, is stable when Fd > IP₂, where IP₂ is the second ionization potential of the system. At intermediate field strengths, IP₁ < Fd < IP₂, the stable solution is given by a state that has mixed neutral and doubly ionized character. It is interesting to note that RHF only allows for the combined ionization of the electron pair and the intermediate state associated with single ionization is never reached. Hence, we see that restricted HF (RHF) cannot recover the exact behavior shown in Fig. 2 for the ground state, as it only connects the neutral solution (with the electrons paired up on the atom) and the doubly ionized state.

The analysis of the unrestricted case is more intricate and interesting. Two independent variational parameters (θ₁, θ₂) control the mixing of the basis functions for the independent sets of α and β orbitals (Table I). The stationary condition on the energy expectation value for the UHF determinant leads to four different solutions for the (θ₁, θ₂) pairs. The solutions (θ₁, θ₂) = (0, 0) and $(θ1,θ2)=(π2,π2)$ reduce to RHF, corresponding to double occupation of the atomic and ghost orbitals, respectively. The former case is stable when Fd < IP₁, whereas the stability for the latter is achieved when Fd > IP₂. The difference in the UHF case is the possibility of achieving intermediate stable solutions $(θ1,θ2)=(0,π2)$ and $(θ1,θ2)=(π2,0)$ for the intermediate range of field strengths IP₁ < Fd < IP₂. These extra UHF solutions are characterized by having one electron occupying the atomic orbital and the other one localized on the ghost function. Therefore, we see that UHF is able to characterize the singly ionized state. This ability, however, comes at a cost: in an analogy to the common problem of bond dissociation at the UHF level, these single determinant unrestricted solutions for the intermediate range of field strengths are composed of an equal mixing of singlet and triplet configurations that results in a spin-contaminated state with ⟨S²⟩ = 1.¹³ Nonetheless, UHF is qualitatively able to describe the static field ionization behavior expected from our analysis of the exact results (Fig. 2), showing the three distinct regimes corresponding to the three possible levels of ionization. In addition, the large separation between $A$ and $C$ leads to a vanishing singlet–triplet gap, so the contamination has no direct impact on the energetics. Finally, we note that the 1:1 quantitative mapping between the position of the “kinks” in the energy curve as a function of field strength and ionization potentials (IPs) is only possible if the conditions outlined in Table I are satisfied. The model can be easily modified to account for the electrostatic repulsion between $A$ and $C$ and the kinetic energy of the populated ghost function (⁠$〈C|T̂|C〉$⁠), leading to a better numerical agreement for the IPs. Alternatively, a decomposition of the total energy into atomic, continuum, and interaction contributions could lead to better quantitative agreement for the IPs. Qualitatively, however, the main features do not depend on the specific nature of the ghost functions.

We move onto an analysis of the performance of other wavefunction-based methods on a larger basis set to describe the static field ionization of He. We should point out that our emphasis is, once again, on the difference between how different flavors of self-consistent field (SCF) solutions can capture ionization effects in this discretized continuum model. In this case, He is represented by the larger aug-cc-pVTZ basis and the continuum is represented by a track of 18 1s hydrogenic STO-3G ghost basis functions spanning the range between −4.5 and 4.5 Å (with an equal spacing of 0.5 Å between each ghost center) along the direction of the field. Figure 3(a) summarizes the results in this larger basis set.

Just as expected from the analytical model, RHF cannot qualitatively describe the correct behavior for single-electron ionization due to the constraint that the two electrons are paired. For field strengths smaller than 0.175 a.u., we observe the neutral (two electrons localized around the atom) ground state, and for field strengths greater than 0.36 a.u., we observe the doubly ionized state of the atom. We note a smooth transition between the neutral atom and its atom dication state for intermediate field strengths. This can also be observed when we analyze how the total electric dipole [Fig. 3(b)] of the system and the charge of the He atom [Fig. 3(c)] change as a function of the applied field. UHF, on the other hand, captures the essential features expected from the exact solution: the energy curve shows three distinct regions separated by first derivative discontinuities. These kinks are associated with spin polarization effects as illustrated by the transition between ⟨S²⟩ = 0 and ⟨S²⟩ = 1. The sharpness of this transition (i.e., how fast the spin polarization occurs) is highly dependent on the size of the basis set: for smaller basis sets, ⟨S²⟩ seems to continuously change from 0 to 1, which indicates another remarkable similarity between the analysis presented in this work and the nature of the Coulson–Fischer point for bond dissociation. For larger basis sets, the spin polarization transitions appear to be much more sudden and discontinuous. We next consider the behavior of electric dipole and He charge for the UHF solution [Figs. 3(b) and 3(c)]. We do, however, note that the numerical magnitude of the change in the dipole moment is directly proportional to the size of the box created by the ghost functions. Both observables indicate that UHF qualitatively captures the three distinct regimes expected for the ionization of helium: the polarization of the electron pair on the neutral atom, the singly ionized state, and the doubly ionized state. Again, the ability of UHF to describe this singly ionized state comes at the cost of spin polarization as indicated by ⟨S²⟩ = 1 for field strengths in the range that generates one unit of positive charge on helium.

Moreover, we briefly analyze the MP2 (Fig. 4) performance for strong field ionization of He using restricted (RMP2) and unrestricted (UMP2) orbitals. As expected, since UHF is a better reference to describe all possible charge regimes for our systems, UMP2 recovers all the qualitative features of the exact solution. By contrast, RMP2 seems to oscillate above and below the exact solution: for field strengths closer to the ionization points, RMP2 energies are, as expected, above the exact values and seem, at first, to follow the same pattern as the RHF curve. However, RMP2 drops below the exact curve for fields in the mid-range for the singly ionized state. This indicates that RMP2 is over-correcting for the bad reference provided by the restricted orbitals in this regime. Nonetheless, this behavior is not as catastrophic as in the RMP2 potential energy curve for the dissociation of H₂, where incipient orbital degeneracies cause a divergent PES as the bond length increases.¹³ Although not explored in the present work, this failure could be a fertile ground for an investigation of the efficiency of orbital-optimized MP2 (OO-MP2) schemes with regularization⁸⁷ to cheaply but accurately recapture some of the correlation effects associated with field ionization.

Our final analysis for the He atom is based on an attempt to remove most of the electric field and spurious discretized continuum contributions to the total energy, accounting only for the purely electronic energy of the atom as a function of the field strength. Within the dipole approximation, we defined an internal energy operator

The electronic energy of the combined atom–continuum system can then be calculated via

Equation (2) includes three main components: (i) the electronic energy of the isolated atom (⁠$Eelecatom$⁠), (ii) the electrostatic interaction between the discretized continuum and the He atom (E_elst), and (iii) the electronic energy of the populated discretized continuum (⁠$Eelecghost$⁠).

In the limit of two non-interacting subsystems (the isolated atom and the discretized continuum represented by the ghost functions), the super-system wavefunction is given by the product of the two individual wavefunctions and the super-system density matrix (D) can be written as a direct sum of the atomic density matrix (D^A) and a density matrix associated with the ghost functions (D^G). For our setup, however, the definition of the atomic subsystem is somewhat fuzzy since, in the complete basis set limit for the atom, the most diffuse basis functions will overlap with the nearest ghost centers added to represent the continuum electrons. This situation could be ameliorated by increasing the distance between this nearest ghost center and the atom or by noting that, after ionization, electrons will tend to occupy the ghost functions closer to the edge of the artificial box used for the discretization of the continuum states. This allows us to have a clear definition of the ghost subsystem, and we can then define the atomic contribution as the difference between the total energy and the energy contributions from the ghost basis functions.

By following this approach, D^A and D^G are defined as blocks of the super-system density matrix associated with the electron density represented by basis functions localized around the atom and with the density represented by basis functions centered in the track of ghost centers, respectively. The electrostatic interaction can then be calculated as

In Eq. (3), we separated the contribution of the attractive interaction between the electrons populating a ghost function and the He nucleus (which is accounted by the ghost function block of the nuclear attraction one-electron integrals, V^G) and the repulsive screening between the continuum-like electron density (D^G) and the remaining electron density localized around He (D^A).

The electronic energy of the populated discretized continuum is expressed as a trace over its subsystem density matrix, D^G,

where T^G is the kinetic component of the one-electron integrals for the ghost basis functions and F^G is a Fock matrix in which we replaced the usual one-electron integrals by only its kinetic component (since for the discretized continuum subsystem, there is no nuclear attraction).

Equation (4) is clearly correct in the non-overlapping limit. In the presence of overlap, Eq. (4) corresponds to the embedded mean-field theory^88,89 decomposition of the total energy of two composite systems (in our case, He and ghost centers) into the energies of each subsystem and their interaction energy. It is worth noting that we have not explicitly included other energy contributions arising from the mixed atom–ghost density contributions or atom–ghost exchange effects in Eq. (2). Implicitly, these contributions are accounted for in the definition of $Eelecatom$ [Eq. (5)] and are small due to large separation between the atom and the ghost function on the edge of the track.

Finally, the quantity of interest $ΔEelecatom$ is given in terms of the atomic subsystem electronic energy at a given field strength

relative to the corresponding quantity at zero field,

Based on our previous analysis, this approach should be more suitable for the UHF solution in which the separable density assumption holds to a greater extent as it can properly account for ionization. Figure 5 shows that the kinks in the total UHF energy [Fig. 3(a)] and the dipole and charge discontinuities for the UHF solution [Figs. 3(b) and 3(c)] are, once again, intrinsically related to ionization events, as we can directly recover information about the first and second ionization potentials (represented by the dashed lines in Fig. 5) of He by removing the spurious field and ghost function contributions to the total energy. Finally, the analysis based on the decomposition of the density matrix into its atomic and ghost contributions allows us to obtain more information about the electronic structure of each subsystem: an eigenvalue decomposition of the combined atom–ghost, atomic, and ghost density matrices indicates the number of electrons for each subsystem, as illustrated in Fig. 6.

In Sec. III B, we thoroughly discussed how spin polarization is essential to qualitatively describe ionization by static fields in our model in which the continuum was discretized by a track of ghost functions placed along the direction of the applied electric field. At least for the He atom, we concluded that UHF properly captures, at the expense of spin contamination, the features associated with ionization within this limited model: we observe kinks in the potential energy curves as a function of field strength that are accompanied by discontinuities in the total electric dipole moment and the charge on the He atom. Now, we move onto investigate if this simplest flavor of spin polarization (i.e., allowing α and β orbitals to have different spatial parts) is enough to capture the features of atomic ionization in a heavier atom. Figure 7(a) shows the RHF and UHF energy for the Ne atom. Once again, we see characteristic spin polarization and kinks in the UHF potential energy surface that are associated with ionization. We performed the same previous analysis for the electric dipole moment [Fig. 7(b)] and Mulliken charges on the Ne atom [Fig. 7(c)], and they both agree with our previous discussion: at the expense of losing information about the total spin of the system (i.e., spin contamination and our wavefunction not being an eigenstate of the $Ŝ2$ operator), we can now ionize an odd number of electrons for Ne through UHF, whereas those states are not accessible through a RHF reference.

The electronic structure of Ne is more complicated than He, which also allows us to analyze other qualitative features of the ionized states of the system as the strength of the applied field increases. With this in mind, it is important to point out what one would expect for each of the Ne ionized states within this limited discretized continuum model that we used in our analysis. Table II summarizes simplified electronic configurations expected for Ne, Ne⁺, Ne²⁺, and Ne³⁺ and the track of ghost functions (with net charge 0, −1, −2, −3) that represent the continuum-like states. Even though a proper analysis of the nature of the real ionized states of the atom is hindered due to the single determinant nature of HF theory, one can note that the S₁ and S₂ states in Table II corresponds to one of the components of the physical (field-free) ²P state of the Ne cation and ³P state of the Ne dication, respectively. Our goal now is to see how we can obtain qualitative information about these physical states through Hartree–Fock theory within our limited model for ionization.

For the first ionization, we expect one unpaired electron at both the atom and the ghost functions, leading to a spin-contaminated UHF state in a similar fashion to what was described for the He case (Fig. 6). Figure 8 summarizes some of the features of the UHF for a field strength of F = 0.25 a.u., which lies within the region for single ionization. We can identify two unpaired electrons for the whole system [Fig. 8(a)]. Performing an eigenvalue analysis on the partitioned density matrices for the atom subsystem (D^A) and the continuum/ghost subsystem (D^G) to obtain the natural orbital occupation numbers (NOONs) reveals that, as expected, each subsystem has a single unpaired electron [Figs. 8(b) and 8(c)]: an α electron for the discretized continuum and an excess β electron for the atom [Fig. 8(d)], keeping M_s = 0 for the contaminated ground state.

The doubly ionized state in He simply corresponds to the pairing of the two electrons in one of the ghost functions, leaving behind just the atomic nucleus (Fig. 6, the curve indicated by a field strength of F = 0.5 a.u.). For Ne, however, we expect the second ionized electron to pair with the first one and localize around a ghost center, while Ne becomes a triplet state characterized by two unpaired electrons with the same spin, and consequently, total M_s = ±1 (see Table II). This state, however, cannot be achieved by following the UHF potential energy curve, which is constrained to the M_s = 0 subspace. Instead, the UHF solution for this doubly ionized state of Ne is characterized by a ⟨S²⟩ = 1 spin-contaminated state comprised of α electron and β unpaired electrons on the atom and no unpaired electrons for the track of ghost functions, as illustrated in Fig. 9. Some physical insights can be gained here by analyzing which states of the field-free Ne dication are included in the UHF M_s = 0 field dependent solution in the doubly ionized regime: the external field, breaking spherical symmetry, mixes the ¹S and ¹D states of Ne²⁺. Moreover, due to spin symmetry-breaking in UHF, we observe that one of the components (M_s = 0) of ³P is also included in the mixture. It should be noted that such a behavior also highlights one of the limitations of our approach to describe the ionized flux: by constraining electron density to escape the atom through a unique track of ghost functions, we observe a driving force towards pairing of the ionized electrons. This prevents neon from achieving the desired and expected triplet configuration. A possible, yet not explored in this work, way to minimize this issue would be to expand our description of the discretized continuum by adding more ghost functions into two separate parallel tracks (as illustrated by the gray crosses in Fig. 1), which would possibly allow for the delocalization of the ionized electrons, leaving the atom subsystem with M_s = ±1 and the ghost subsystem with M_s ∓ 1, while maintaining total M_s = 0 for the UHF solution.

Alternatively, since the main interest should be to capture as much physical insights for the atomic subsystem as possible within our limited model, it is possible to remove the M_s = 0 constraint in UHF by exploring the manifold of generalized Hartree–Fock (GHF) solutions. In GHF, all the exact spin symmetries of the system (namely, $Ŝ2$ and $Ŝz$⁠) can be broken by taking spin–orbitals that no longer belong to two independent sets of α and β orbitals but now are linear combinations of orbitals in both sets. This allows the spin quantization axis of individual electrons to rotate and not necessarily be parallel to each other in order to achieve the variationally lowest possible energy state. This extra flexibility does not mean that the GHF solution cannot be obtained through UHF though. A GHF wavefunction is said to be non-collinear (and sometimes referred to as a true GHF solution) if its spin cannot be quantized along some axis, and this illustrates a real UHF → GHF instability of the orbital Hessian.^90,91 This is the case for the double bond dissociation in CO₂^90,92 and for systems presenting spin-frustration.^93,94 If such a spin quantization axis exists, the GHF solution can be found by exploring UHF wavefunctions with different spin multiplicity or lying within different M_s manifolds. One can then find the difference between the number of up and down electrons by projecting the spin operators in this quantized axis.⁹⁰ In this context, we see that GHF can provide a convenient way to adiabatically follow the lowest UHF solution independently of M_s constraints.

Figure 10(a) illustrates the behavior of the different M_s UHF solutions compared to GHF for the strong field ionization of the neon atom within the race track of ghost functions to represent discretized continuum states that we have been using so far. Even though we cannot note significant differences between UHF and GHF in terms of their description of charges, dipoles, and the total number of unpaired electrons for the supersystem, the GHF potential energy curve follows the lowest UHF solution, highlighting its collinear nature. For instance, for field strengths lower than F = 0.15 a.u., we observe that the M_s = 0 UHF solution is the one with the lowest energy and that there is no difference between this solution and GHF, whereas the other M_s states are higher in energy. For electric fields ranging from F = 0.15 and F = 0.309 a.u. (single ionization), however, the M_s = 0 and M_s = 1 states become degenerate, characterizing the RHF → UHF instability previously discussed. For the regime between F = 0.309 and F = 0.60 a.u. (which corresponds to the doubly ionized regime), the M_s = 1 UHF solution becomes the lowest in energy, illustrating the transition between different M_s subspaces that is captured naturally through GHF. Our first UHF analysis was constrained within the M_s = 0 subspace, and therefore, the wrong qualitative behavior (one α and one β electron) was observed for the double ionized state of the atom. GHF, however, abolishes the M_s = 0 constraint, and we can access the relevant doubly ionized state lying in the M_s = ±1 manifold [Fig. 10(b)] and Ne achieves the expected configuration with two unpaired α electrons, which better resembles the M_s = ±1 component of the physical ³P field-free state of the Ne dication.

Finally, even though the plateau regions in Figs. 7(b) and 7(c) are well characterized by collinear GHF solutions that represent UHF wavefunctions with different M_s values, the transition between these ionization states might not present the same behavior. At the onset of ionization, we suspect that, similarly to what was observed for the RHF → UHF instability, one would need a non-collinear GHF solution to properly describe the transition between one ionized state and the other. However, due to the dependence of the sharpness of this transition on the size of the artificial box delimited by the ghost functions, we did not pursue further investigation of this question. Hence, we note that GHF was used merely as a tool to access the lowest UHF state possible without any initial constraints on the number of α and β electrons.

In conclusion, the present work explored an analogy between the Coulson–Fischer points that have been extensively characterized for bond dissociation problems and static strong field ionization phenomena. Such an extension is, at first, hindered by the fact that ionization implies an outgoing electron flux that is poorly described by the atom-centered Gaussian basis sets commonly used in electronic structure packages.⁸⁶ Hence, our basic “race-track” model was introduced to account for this ionized flux by adding ghost basis functions spanning an appropriate region of space to allow for some characterization of the continuum-like states of the ionized system.

A preliminary analytical analysis of the ionization problem of He in a minimum basis with a single distant ghost function hinted that ionization events, which would be characterized by the localization of the electron density at one of the auxiliary ghost centers of the “race-track,” would not be well characterized by a restricted Hartree–Fock approach. In particular, the RHF solution could describe the neutral and the doubly ionized states of the system but would lack any information about the singly ionized state, given its electron-pairing constraint. This analysis also indicated the existence of a RHF → UHF instability for the singly ionized range of field strengths that was akin to the triplet instability that leads to spin polarization in the bond dissociation problem. The predictions from this toy model were validated by computations in a larger basis set. Moreover, the UHF potential energy surface as a function of field strength reveals kinks that, after an analysis of the charge population of the atom, the dipole moment of the system, and the expectation value of the $Ŝ2$ operator, indicate the existence of an analog of the Coulson–Fischer point for static strong field ionization. We have also made an effort to eliminate possible spurious contributions stemming from our limited model of the discretized space by devising a partitioning scheme based on taking traces of the RHF and UHF density matrices over the appropriate basis functions centered only around the atom. Such a scheme allowed us to, once again, determine the character of each state of the atom for different field strengths, as well as to recover information about the ionization energies of He. While RHF cannot represent the singly ionized state, UHF recovers the first and second ionization energies of He at the cost of breaking spin symmetry and being spin-contaminated.

We then extended our analysis to a more complicated system: the Ne atom. The main takeaways are still valid: by constraining the HF solution to be an eigenfunction of $Ŝ2$⁠, i.e., constraining ourselves to the RHF manifold of solutions, we cannot describe ionization of an odd number of electrons as the strength of the external field increases. By exploring the solutions contained in the spin polarized UHF manifold, on the other hand, we have a qualitatively correct description of the atomic charge associated with these ionized states. However, due to a strong drive to pair the electrons at the end of the auxiliary “race-track” of ghost functions and the constraint imposed on the value of $⟨Ŝz⟩$ for UHF solutions, the results do not correspond to the expected Hund’s rule configuration of Ne for some field strengths. Thus, we also analyzed the use of GHF solutions to obtain a better qualitative description of the static strong field ionization of Ne. We observed that, while maintaining the same behavior as UHF for the atomic charges, GHF can properly rotate spins at the cost of losing information about the expectation value of $Ŝz$ in order to give the proper high-spin description of the unpaired electrons on the atomic subsystem. Exploring different kinds of symmetry-broken solutions, such as allowing complex polarization in Hartree–Fock (cRHF, cUHF, and cGHF),^95–98 could also prove potentially useful to obtain a qualitatively correct description of strong field ionization of more complex systems.

Even though the results and analysis presented in the current work were obtained for the case of an applied external static electric field, it can be argued that some of the consequences of imposing spin constraints on the Hartree–Fock solutions could also lead to a qualitatively wrong description of the real-time dynamics of the system when a physically meaningful laser pulse is used. For the more common problem of bond-stretching in electronic structure, it has already been shown that using a qualitatively bad HF solution as the starting point for more elaborate methods that account for dynamic electron correlation could lead to failures in obtaining a good description of potential energy surface, for example. In this sense, work toward extending the analysis here to real-time time-dependent electronic structure is underway in our group. To achieve more physically meaningful results, a better approach to represent the continuum states, perhaps through a mixed Gaussian-plane wave basis set or complex absorbing potentials, is needed.^51,99

See the supplementary material for the raw data (energies, dipole, and charges) organized in XLXS format and additional plots (PDF).

This research was supported by the Director, Office of Science, Office of Basic Energy Sciences, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. L.A.C. would like to thank Dr. Josh Cantin and Juan Arias-Martinez for stimulating discussions. M.H.-G. is a part-owner of Q-Chem, which is the software platform in which the computations described here were carried out.

The data that support the findings of this study are available within the article and its supplementary material.