Self-Consistent Field Methods for Excited States in Strong Magnetic Fields: a Comparison between Energy-and Variance-Based Approaches

: Self-consistent ﬁ eld methods for excited states o ﬀ er an attractive low-cost route to study not only excitation energies but also properties of excited states. Here, we present the generalization of two self-consistent ﬁ eld methods, the maximum overlap method (MOM) and the σ - SCF method, to calculate excited states in strong magnetic ﬁ elds and investigate their stability and accuracy in this context. These methods use di ﬀ erent strategies to overcome the well-known variational collapse of energy-based optimizations to the lowest solution of a given symmetry. The MOM tackles this problem in the de ﬁ nition of the orbital occupations to constrain the self-consistent ﬁ eld procedure to converge on excited states, while the σ - SCF method is based on the minimization of the variance instead of the energy. To overcome the high computational cost of the variance minimization, we present a new implementation of the σ -SCF method with the resolution of identity approximation, allowing the use of large basis sets, which is an important requirement for calculations in strong magnetic ﬁ elds. The accuracy of these methods is assessed by comparison with the benchmark literature data for He, H 2 , and CH + . The results reveal severe limitations of the variance-based scheme, which become more acute in large basis sets. In particular, many states are not accessible using variance optimization. Detailed analysis shows that this is a general feature of variance optimization approaches due to the masking of local minima in the optimization. In contrast, the MOM shows promising performance for computing excited states under these conditions, yielding results consistent with available benchmark data for a diverse range of electronic states.


INTRODUCTION
Atoms and molecules may exhibit exotic chemistry in strong magnetic fields.−4 Since this range of magnetic fields cannot be investigated on Earth, quantum chemistry is an essential tool for understanding chemistry under these conditions.In this context, Lange et al. showed that a new molecular bonding mechanism, the perpendicular paramagnetic bonding in H 2 and He 2 , should exist in strong magnetic fields. 5However, computing properties under such conditions beyond the perturbative regime needs particular machinery.In addition to the modifications of the molecular Hamiltonian due to the external field, gauge-origin independence of energies and properties has to be ensured.This may be achieved by means of London atomic orbitals (LAOs) but leads to a computational cost that is higher than that of calculations with standard Gaussian-type basis sets. 6,7xcited states are of fundamental importance in chemistry, not only in the absence of magnetic fields where they are essential in many areas such as photochemistry but also in the presence of magnetic fields.From the theoretical point of view, their study has led to the development of a wide range of methods.−25 The computation of excited states in the presence of magnetic fields is essential for understanding and interpreting the spectra observed from stellar objects where such strong fields exist.−29 Studies including electron correlation have also been presented at the full CI (FCI) level 30−34 and by means of quantum Monte-Carlo by Jones, Ortiz, and Ceperley. 35More recently, linear response HF 36,37 and EOM-CCSD 38−40 have been extended to the strong-magnetic field regime.However, the complex algebra required by the use of LAOs makes the adaptation of these sophisticated methods often a complicated task, in addition to significantly increasing their computational cost.
In this context, the prospect of computing excited states through self-consistent field (SCF) approaches is attractive.Indeed, using a mean-field approach presents several advantages: (i) it provides a simple picture of excited states which may be interpreted in the same way as the ground-state (GS) ones, (ii) it yields a reference wave function to which correlated methods may be applied, (iii) unlike linear response-based methods, it could provide a correct description of many types of excitations (including, e.g., double excitations and charge transfer), 41,42 and (iv) the computational cost should be no more than that for the GS calculation with the same method.However, the main difficulty is to find a way to circumvent the variational collapse problem, since a variational optimization of the energy will generally lead to the lowest energy state of a given spin and spatial symmetry.
−45 The principle of the Δ-SCF method is to optimize a non-auf bau configuration by promoting an electron from an occupied orbital to a virtual orbital of a GS wave function.However, since the order of molecular orbitals (MOs) may change between SCF iterations, tracking the desired configuration is often difficult.The MOM attempts to overcome this issue by defining the occupation according to the orbitals that have the most overlap with those from the previous iteration.−56 Despite its success, the MOM suffers from deficiencies, among which is the collapse to other solutions during the SCF procedure, which cannot always be prevented, and other strategies have been employed to variationally optimize excited states through energy optimization-based schemes.Colle et al. 57 proposed to variationally determine excited states by constraining the wave function to be orthogonal to the GS.While this constraint is more complicated to implement, it has led to several methods either in HF 58−60 or DFT, 61−65 with which excited states may be calculated.−74 These include the possibility to overcome the variational collapse of the energy optimization with direct optimization techniques.In this context, Hait and Head-Gordon proposed an extension of direct orbital optimization to optimize excited states by minimizing the square of the orbital gradient to avoid variational collapse; this is termed the squared-gradient minimization (SGM) approach. 71Furthermore, the use of quasi-Newton optimization approaches in this context has been exploited by Levi et al. 73,74 In order to improve further on the Δ-SCF method, as the MOM does by tracking the orbital occupations, the statetargeted energy projection method enforces the computation of non-auf bau configurations by restricting the occupied-virtual orbital rotations by means of a level shift in the definition of the Fock matrix. 72−68 An alternative possibility to overcome the variational collapse of the energy optimization is to instead optimize a different measure of the wave function.Based on the minimization of the variance, 75,76 Ye et al. have recently introduced the σ-SCF method for computing either the GS or excited states. 69,70This method proceeds through two successive minimizations with (i) an energy-targeting minimization to target the desired state and to provide an initial guess for (ii) a variance minimization of the state.While the variance is a quantity commonly used in quantum Monte-Carlo methods, less attention has been paid to its exploitation elsewhere in electronic structure theory.
The present paper is organized as follows: after introducing the theoretical background of quantum chemistry in strong magnetic fields, Section 2 presents the MOM and σ-SCF method and discusses their implementation in strong magnetic fields.Numerical applications are presented in Section 3, and we compare the performance of both methods in the presence of a magnetic field.The helium atom, the hydrogen molecule, and the methylidinium ion have been chosen as prototypical examples since they are representative of the species that may exist in the vicinity of magnetic white dwarf stars, 1−4,77−79 and the results are discussed in comparison with higher levels of theory, such as FCI and EOM-CCSD calculations.Finally, in Section 4, we discuss the practical limitations of variance-based optimization and explain their origin.

THEORY
2.1.Quantum Chemistry in Strong Magnetic Fields.The electronic non-relativistic Hamiltonian for an N-electron system in the absence of a magnetic field is given by with the momentum operator p i = −i∇ i for the electron i with coordinates r i .Here, Z K and R K are the charge and position of the nucleus K, respectively.The effect of a magnetic field B is introduced by a modification of the kinetic energy operator with the vector potential A, defined such that B = ∇ × A, and with a spin-dependent term B•S as which can be re-written as Due to the dependence of the vector potential on the gauge origin R O in a uniform magnetic field, eigenfunctions of the Hamiltonian in eq 3 become dependent on the choice of R O , which is completely arbitrary.This behavior cannot be modeled using a finite basis set, except by explicitly including the gauge-origin dependence; this approach is adopted using LAOs, which consist of standard Gaussian-type functions multiplied by complex phase factors 6,7 These exhibit the correct behavior to the first order with respect to the magnetic field and ensure that calculated observables are gauge-origin-invariant.
2.2.MOM and Initial MOM.In order to variationally determine excited states, the MOM provides a way to track the desired non-auf bau configuration between SCF iterations, since the change in the order of MOs at each step may complicate this task. 41At the nth iteration, one may define the overlap matrix between the current MOs and those from the previous iteration as with C being the matrix of MO coefficients and S being the overlap matrix in the atomic orbital (AO) basis.Thus, defining the projection where i and j run over the M occupied MOs of the (n − 1)th and nth iterations, respectively, and p and q run over the N AOs, allowing the current set of MOs to be occupied by selecting those which have the greatest overlap with those from the previous iteration.Another possibility is to choose the initial-guess orbitals as the reference set of MOs This method, called as the initial MOM (IMOM), may enable the computation of some states often difficult to access with the MOM, such as doubly excited states. 42In order to reach convergence to the desired excited state, the starting orbitals should be as close as possible to those of this solution.Often, orbitals obtained from a GS calculation provide an adequate guess for MOM and IMOM calculations when occupied with the required non-auf bau configuration.This choice has been used throughout the current work.
Therefore, the implementation of the MOM and IMOM is very simple since it mainly consists of a small modification to the SCF process.Furthermore, this simplicity is of particular interest for our purpose since the generalization of these methods to the strong magnetic field framework is straightforward.
2.3.Variance Minimization and the σ-SCF Method.The variance functional provides an alternative measure for the accuracy of a trial wave function Φ. Expressing this normalized wave function in the basis of the exact eigenstates Ψ i of the Hamiltonian (9)   where c i corresponds to the expansion coefficients; its variance may be expanded in terms of the corresponding eigenvalues E i as From this equation, it is readily seen that if Φ is an eigenstate of the system, the variance is equal to zero, and otherwise, S[Φ] measures the departure of Φ from one of the eigenfunctions of the Hamiltonian.Therefore, the variance presents a local minimum for each eigenstate, and a selfconsistent optimization of this quantity should approximate not only the GS but also excited states.Thus, one may expect to compute a desired state if the initial guess is already located in its variational well.The σ-SCF method provides a way to locate the desired well by performing an energy-targeting minimization, which provides the guess density matrix for the subsequent variance minimization.The energy-targeting function is dependent on the user-defined parameter ω. Figure 1 illustrates the meaning of both quantities defined by eqs 8 and 11 for the example of H 2 in the u-3-21G basis, where the prefix 'u-' indicates that the basis set is uncontracted.The black curve corresponds to the energy-targeting minimization, which presents clear wells with respect to ω.Each of these wells is associated with a plateau in terms of variance (blue) and energy (orange).Thus, σ-SCF proceeds through two successive minimizations: (i) W[Φ](ω) (black circles) may be optimized with a given value of ω, providing a guess in the desired well.Once this quantity is optimized, (ii) the minimization of the variance converges to the minimum of the related energy-targeting well (blue).The optimized wave function allows us to calculate the energy of the state (orange).Thus, one may target a particular state by selecting an appropriate value of ω, which should be close to the electronic energy of that state.For further discussion about the physical meaning of W[Φ](ω) and ω, the reader should refer to ref 69.
Although variance minimization-based methods are less popular than their energy optimization counterparts, this approach has been studied in the 60's and 70's 80−83 and is commonly applied in quantum Monte-Carlo methods. 84The main difficulty associated with such methods arises from the square of the Hamiltonian occurring in eqs 8 and 11.Indeed, expressing eq 8 in an orthogonal AO basis gives  (12)   Here, D pq is the element of the density matrix between AOs p and q, F pq is an element of the Fock matrix, and Q pq is an element of the matrix defined as Q = I − D, with I being the identity matrix.Note that two electron integrals (TEIs) are given in an orthogonal AO basis.This expression clearly shows the introduction of 8-index summations over the TEIs.This results in an increase in the computational cost compared to an energy-based optimization such as HF or MOM calculations with a scaling of N ( ) 5 , where N is the number of basis functions, coupled with a large prefactor. 69So far, this has limited the applications of σ-SCF to atoms or diatomic molecules in small basis sets. 69,70However, the size of the basis sets used for atoms and molecules in strong magnetic fields is an important requirement since they must be sufficiently flexible to describe the deformation of the electronic structure due to the field.Furthermore, care should be taken in selecting basis sets for excited states since their electron densities may be more diffuse than those of ground states.The scaling of σ-SCF in addition to the basis set requirements makes the use of the RI approximation attractive; a reformulation of eq 12 using the RI approximation has been developed and implemented in the present work.For the sake of convenience, all working equations and the derivation of the RI-based implementation are presented in the Appendix, along with numerical illustrations of its effectiveness.These developments significantly reduce the computational cost, and the RI-version scales as , with M being the number of occupied orbitals and K being the number of auxiliary basis functions.It is worth noting that the σ-SCF equations are expressed in terms of the Fock matrix and the TEIs, and thus, the introduction of the terms in the Hamiltonian occurring from the magnetic field does not modify the σ-SCF equations since they result only in a modification of the one-electron operator.However, due to the use of LAOs, there is less permutational symmetry for the TEIs 38 and so less opportunity to accelerate calculations by taking advantage of this.

RESULTS AND DISCUSSION
The MOM, IMOM, and σ-SCF have been generalized for calculations in strong magnetic fields and using LAOs in our development platform QUEST, 85,86 where all three methods have been implemented using standard TEIs and also with the RI approximation (see the Appendix for numerical comparisons of RI and conventional results).Although the MOM/ IMOM may be used with DFT, in the present work, they are only applied to Hartree−Fock to enable a direct comparison between the MOM/IMOM and σ-SCF, since the latter is not directly applicable to DFT.All calculations have been performed using the unrestricted formalism, permitting a direct computation of different spin multiplicities.Magnetic field strengths are given in atomic units, where 1B 0 ≈ 2.35 × 10 5 T.
As a proof of concept, in Section 3.1, we initially study the two-electron example of the helium atom in the 6-311G basis under strong magnetic fields and compare the excited-state energies obtained from both the MOM and σ-SCF with the results of FCI calculations, computed using the BAGEL package, 87 and results of CIS calculations obtained using the time-dependent linear response implementation in QUEST.
Following this, in Section 3.2, the energy spectra of the hydrogen molecule were computed with the MOM/IMOM and σ-SCF at zero field and in the presence of magnetic fields parallel and perpendicular to the internuclear axis and compared with CIS spectra.These were computed in the uaug-cc-pVDZ and u-aug-cc-pVQZ basis sets; in general, uncontracted orbital basis sets are used to allow greater flexibility to respond to the effects of the magnetic field.Here, the RI approximation was employed with the AutoAux auxiliary basis, 88 which provides large flexible RI basis sets tailored automatically to each orbital basis used throughout.
Subsequently, the energies of selected states of the hydrogen molecule with a magnetic field parallel and perpendicular to the molecular axis were calculated using the MOM and compared with the results of EOM-CCSD calculations from ref 39 at magnetic field strengths from 0 to 1B 0 .A similar analysis is conducted for the methylidinium ion (CH + ) with a magnetic field applied parallel to the internuclear axis in Section 3.3.In both cases, the RI approximation was used with the u-aug-cc-pVQZ and u-aug-cc-pCVQZ orbital basis sets for hydrogen and carbon atoms, respectively, and the AutoAux auxiliary basis.The same molecular geometries as in ref 39 were used throughout.For CH + , a magnetic field perpendicular to the molecular axis was also studied with the MOM, and the results are compared to those of the FCI calculations from ref 40 with the same molecular geometry.Note that this study uses the more modest cc-pVDZ basis.In all cases, the Cartesian form of these basis sets was used to enable a comparison with refs 39 and 40.
3.1.Helium Atom.The example of the helium atom in the 6-311G basis corresponds to a system with six m s = 0 configurations and three configurations for each m s = ±1.While in the absence of a magnetic field, the m s = ±1 determinants of a given triplet state are degenerate, the spin Zeeman effect that arises by virtue of an applied magnetic field breaks this degeneracy and makes it necessary to consider all m s components.Figure 2 presents the energies of all m s = ±1 components computed with the MOM and σ-SCF, denoted by squares and crosses, respectively, as a function of magnetic field strength.Since the 6-311G basis set only features s-type orbitals for the helium atom, only the spin part of the Zeeman effect is relevant, causing the splitting of both m s components of each triplet state.Over the range of magnetic field strengths for which these were computed (from 0 to 1B 0 ), the diamagnetic term does not seem to be a significant contribution and the energies show a near-linear variation with the field strength.Both methods are in very good agreement, and energy differences between the MOM and σ-SCF calculations at each field strength are below the convergence criterion (all energy differences are below 10 −10 E h ).
Figure 3 (left panel) presents the energy of the m s = 0 components computed using both methods with respect to the magnetic field strength, and the same good qualitative agreement is observed.However, plotting the differences between the MOM and σ-SCF energies with respect to the field for each determinant in the right panel of Figure 3 reveals different features.The quantitative agreement between both approaches for the m s = 0 configurations is less impressive than that for the m s = ±1 components.The ground state (1s α 1s β ) and the highest excited state (3s α 3s β ) present the smallest energy difference, which is non-negligible (about 0.4 mE h ), while the energy given by both methods for the first excited state differs by about 4 mE h at zero field.The difference seems to be neither related to the nature of the excited configuration (open-or closed-shell configuration) nor the position of the states in the energy spectrum; those with the two largest differences are the first and fourth excited configurations.It may be worth noting that the three states presenting the largest differences include an occupied 2s orbital, while all differences between both methods decrease when the magnetic field strength increases.Finally, in contrast to ref 69, the MOM solutions do not necessarily appear lower in energy than the σ-  SCF ones, and indeed, it may be seen in Figure 3 that this is only the case for the ground and fourth excited states.
Table 1 presents the mean absolute error (MAE) between both SCF methods and FCI calculations over a range of magnetic field strengths from 0 to 1B 0 for all states.While the MOM and σ-SCF give a consistent description of the GS, their MAEs compared with FCI calculations are relatively large, and this quantifies how much the GS is stabilized by electron correlation compared to its mean-field description.In contrast, the triplet configurations are relatively well described by both SCF methods, with an MAE of about 0.2 mE h for the configurations occupying the 1s2s orbitals and about 3.0 mE h for the configurations occupying the 1s3s and 2s3s orbitals.It should be noted that the energies of the m s = 0 triplet components are obtained by simply removing the Zeeman contribution from the total energy of the m s = ±1 component.Interestingly, both closed-shell excited states have an MAE with the same magnitude as that for the triplet states.In addition, the largest MAEs occur for the open-shell singlet states, with values larger than 10 mE h .However, the energy of these states is obtained using the simple projection technique 45 = − ↑↓ ↑↑ which becomes in the presence of an external magnetic field due to the spin Zeeman energy contributions present in the m s = ±1 components.Evaluating open-shell singlet-state energies through this kind of approach is only appropriate when the singlet state and the corresponding triplet state have similar charge distributions.The same idea is also exploited to compute magnetic exchange couplings by means of broken symmetry determinants (i.e., singlet−triplet gaps in magnetic compounds). 89,90However, this approximation may be no longer valid for some excited states in which the orbital relaxation leads to significantly different total densities between the singlet and triplet states.For comparison, MAEs for CIS calculations are also presented in Table 1.We have confirmed that the errors relative to FCI do not vary substantially for each approach as a function of field and so may be reasonably summarized by the MAEs.The CIS results are presented using both the lowestenergy 1s β 2s β and 1s α 1s β SCF solutions as the initial determinant.The former provides a good reference for the components of the triplet 1s3s and 2s3s excited states, yielding MAEs essentially identical to those of the MOM and σ-SCF procedures, which optimize the orbitals in these configurations.In contrast, using the 1s α 1s β reference to calculate the 1s3s triplet-state components leads to MAEs that are substantially larger, indicating the importance of orbital relaxation effects.A similar effect is observed for the 1s2s triplet components using the same reference.As expected, the open-shell singlet states are not well described by CIS.Interestingly, in most cases, the CIS MAEs for these states are significantly higher than those

Journal of Chemical Theory and Computation
for the MOM and σ-SCF, indicating that orbital relaxation may significantly reduce errors for these states, despite the fact that these methods do not provide spin pure states, and the energies are obtained through the projection of eq 14.Nonetheless, substantial errors remain, reflecting the wellknown difficulty in describing open-shell singlet states with single-determinant wave functions.
3.2.Hydrogen Molecule.An important computational requirement for atoms and molecules in strong magnetic fields is the size of the basis set, which must be sufficiently flexible to represent the distortion of the orbitals due to a strong magnetic field. 38,91Such calculations then usually require the use of large and uncontracted basis sets, such as the u-aug-cc-pVQZ basis set used in this work.Previous calculations with σ-SCF have been limited to very small basis sets. 69,70In the Appendix, an implementation using the RI approximation is presented.This allows calculations to be performed at lower cost, with more appropriate basis sets.Before comparing with reference calculations computed with the EOM-CCSD method, the MOM and σ-SCF calculations are first compared in larger basis sets.
Figure 4 presents the energy spectrum computed using CIS, the MOM/IMOM, and σ-SCF for the m s = −1 components for the hydrogen molecule in the absence and in the presence of a magnetic field of 1B 0 parallel and perpendicular to the molecular axis.Results in the two basis sets are compared: first in the modest u-aug-cc-pVDZ basis (left) and second in the larger u-aug-cc-pVQZ basis (right), more appropriate for calculations in strong magnetic fields.In both cases, applying a magnetic field stabilizes the m s = −1 components due to spin Zeeman effects.Furthermore, it may be seen that the energy spectrum in the parallel case presents more states than at zero field; this observation may be attributed to the effect of the orbital Zeeman contribution, which removes the degeneracy between orbitals with angular momentum greater than zero.In the same way, the number of states is even larger in the perpendicular case since the symmetry of the system is further reduced by the application of a magnetic field in this orientation.
Considering first the spectra computed in the u-aug-cc-pVDZ basis, shown on the left of Figure 4, it can be seen that the values obtained from CIS, the MOM/IMOM, and σ-SCF in the absence of a field and in parallel and perpendicular fields are largely consistent between the methods.This confirms that most singly excited states in the CIS spectra can be obtained with the MOM/IMOM and σ-SCF in this basis set.Furthermore, some states are accessible with the MOM/ IMOM and σ-SCF, which do not appear in the CIS spectra, for example, those with energies of around −0.22 and −0.24 E h for H 2 in a parallel magnetic field; analysis reveals that these are multiply excited states, which cannot appear in CIS spectra.However, the results presented here show only those states readily obtained from the standard MOM/IMOM approach using the ground-state orbitals as the guess from which to construct the non-auf bau solutions.
Despite the generally good agreement, the correspondence between the MOM/IMOM and σ-SCF spectra and CIS is not exact; a relatively small number of states present in the CIS spectra were only accessible with either the MOM/IMOM or σ-SCF but not both.In addition, several σ-SCF calculations converge to unphysical solutions, highlighted in gray in Figure 4.The presence of these spurious solutions has been highlighted by Ye and Van Voorhis in ref 70, while Goodisman interpreted similar observations when using variance-based methods as a probable mixture of several states. 92These unphysical solutions in the σ-SCF spectra can be problematic when they lie close in energy to a physical solution, for example, as indicated by the gray line close to the ground-state solution in the absence of a magnetic field.In practical calculations, unless the value of ω is carefully scanned and the resulting solutions are analyzed in detail, these cases could easily be misassigned as one of the physical solutions.This problem can be exacerbated in the presence of a magnetic field since this may lower the symmetry of the system; for example, additional unphysical solutions are visible around the first excited state in the parallel field and around the solution at −0.4 E h in the perpendicular field.
The spectra obtained using the larger u-aug-cc-pVQZ, shown in the right panel of Figure 4, present a somewhat different picture.The comparison between the MOM/IMOM and CIS spectra is similar to that observed in the smaller basis, with the great majority of states in the CIS spectra accessible with the MOM/IMOM, while the MOM/IMOM spectra contain additional lines corresponding to multiple excitations.However, the σ-SCF spectra are distinctly more sparse than their CIS and MOM/IMOM counterparts.Indeed, it may be seen that not all states computed via the MOM/IMOM approach are accessible with σ-SCF.In the absence of a field, even the GS cannot be computed through the variance-based method, while several low-lying excited states are missing in the presence of the perpendicular magnetic field.It is worth noting that although the energy-targeting step usually correctly approaches the desired state, the subsequent variance minimization converges to a different state.The fact that some states are impossible to access through variance minimization in the larger basis set will be analyzed in more detail in Section 4. A particularly striking example is that of the ground state of triplet H 2 in the absence of a field.For Δ-SCF approaches, this would be considered as the most trivial case since it is the lowest of the given spin multiplicity.The relative sparseness of the spectrum obtained with σ-SCF in large basis sets in addition to the presence of spurious states makes its use difficult.Indeed, the comparison with higher levels of theory such as FCI, providing spin pure states, requires the definition of spin-purified energies, which are obtained through eq 14 in this work.In this projection technique, several determinants are necessary to describe open-shell singlet states, and frequently, one of them may be missing with σ-SCF.For this reason, more detailed comparisons with EOM-CCSD calculations are made only for the MOM calculations.
The lowest 1 Σ g + , 1 Π g + , 1 Σ u + , 3 Σ u + , and 3 Π u + states of the hydrogen molecule at the equilibrium geometries of the 1 Σ g + state 32 have been studied and compared to the EOM-CCSD reference calculations recently published by Hampe and Stopkowicz. 39Figures 5 and 6 present the comparison between EOM-CCSD 39 and MOM calculations with a magnetic field parallel and perpendicular to the bond axis, respectively.In order to make this comparison possible, the MOM open-shell singlet states 1 Σ u and 1 Π u are obtained by the projection defined in eq 14.For the parallel case, Figure 5 presents the total energy of the states as a function of the field strength for EOM-CCSD (top left) and the MOM (top right).In both cases, the energies of the 1 Σ g and 1 Σ u states increase with the magnetic field strength due to the diamagnetic term.While the

Journal of Chemical Theory and Computation
1 Π u state is also a singlet, the occupation of a π orbital results in a competition between the orbital Zeeman effect, stabilizing the energy, and the diamagnetic term, leading to an energy minimum around B ∥ = 0.2B 0 .For the triplet states 3 Σ u and 3 Π u , the energy is largely stabilized by the spin Zeeman contribution, and the diamagnetic term does not significantly counteract this effect over the range of magnetic field strength considered here.While the 3 Π u state is higher in energy than 3 Σ u at zero field, the orbital Zeeman effect due to the occupation of the π orbital in 3 Π u results in a crossing between both states at around B ∥ = 0.4B 0 for the MOM and EOM-CCSD, with the 3 Π u state becoming the ground state at higher fields.However, the crossing between the 1 Σ g and triplet states predicted by both methods is slightly different, occurring at B ∥ > 0.4B 0 for EOM-CCSD and at B ∥ < 0.4B 0 for the MOM.The same observation may be made for the crossing between 1 Σ u and 1 Π u , noting that 1 Σ u appears slightly higher than 3 Π u at zero field for the MOM, while this ordering is reversed for EOM-CCSD.Nevertheless, the overall picture is qualitatively and quantitatively well described by the MOM, and the similar profiles of the excitation energies from the 1 Σ g state presented in the bottom left and right panels of Figure 5 confirm this trend.Furthermore, it is worth noting that quantitative agreement could be improved by the inclusion of correlation in the MOM solutions.Since the MOM corresponds only to a modification of the orbital occupations, it may be applied to DFT in the same way as for HF, while post-HF treatments such as CC may be applied to the MOM solution as well.
The perpendicular case is presented in Figure 6, showing the energy of the states with respect to magnetic field strength for EOM-CCSD (top left) and the MOM (top right).Applying a magnetic field perpendicular to the molecular axis reduces the symmetry of the system, from C ∞h in the parallel case to C 2h .Then, the nature of the states changes compared to the parallel case, leading to important differences; this may be demonstrated with the application of a magnetic field to the 1 Σ u + state.Applying a parallel magnetic field to 1 Σ u + results in the 1 Σ u state, which is only affected by the destabilizing diamagnetic term, and its energy increases with the square of the magnetic field strength.In the perpendicular field, the 1 Σ u + state becomes the 1 B u state.In this reduced symmetry, a greater mixing of orbitals may occur; as a result, this state becomes stabilized by the orbital Zeeman effect in addition to the destabilization by the diamagnetic term.This results in the presence of a minimum at B ⊥ ≈ 0.1B 0 .The effect of a perpendicular magnetic field on these states has been already discussed by Hampe and Stopkowicz, and the reader should refer to refs 5, 39, and 93 for a deeper analysis.In our comparison, the MOM correctly represents the energy changes with respect to the magnetic field strength, and it may be seen that the crossings of the 1 A g state with the 3 B u and 3 A u states are qualitatively well represented.As already observed in the parallel case, unlike the MOM, at zero field, the 1 Σ u + state is lower than 3 Π u + with EOM-CCSD, resulting in a crossing between the 1 B u and 3 A u states, which does not occur with the MOM.However, as observed for the parallel case, the MOM provides results that are otherwise generally in very good agreement with those of the EOM-CCSD calculations, which may be further improved by the inclusion of correlation in the MOM solution with DFT.
3.3.Methylidinium Ion.The methylidinium ion (CH + ) is an interesting example since this molecule may exist in the

Journal of Chemical Theory and Computation
vicinity of white dwarf stars, 79 while it cannot be studied by FCI in a large basis set. Figure 7 shows the different results obtained from the EOM-CCSD method and MOM in u-augcc-pCVQZ for the lowest 1 Σ, 1 Δ, 1 Π, 3 Σ, 3 Π, and 3 Δ states in a magnetic field parallel to the molecular axis.Since analyzing the behavior of the different states yields similar conclusions to those drawn in the hydrogen molecule example, the present discussion is mainly focused on the comparison between the two methods.
Figure 7 presents the total energy of each state with respect to the field for EOM-CCSD (top left) and the MOM (top right).While the MOM and EOM-CCSD calculations mainly differ due to the absence of the correlation effects in the MOM, the evaluations of the triplet states are in very good qualitative agreement.Thus, the energies of the 3 Σ, 3 Π, and 3 Δ states as a function of field strength present very similar trends for the MOM and EOM-CCSD, while the EOM-CCSD energies are shifted by the correlation contribution.While this may also be said for the 1 Π state, the situation is different for the 1 Σ and 1 Δ states.Indeed, the 1 Σ state appears higher in energy with the MOM than with EOM-CCSD; this may be seen at zero field, where the 1 Σ state is almost degenerate with the 3 Π state for the MOM.The 1 Δ state presents the same variation in energy with field strength with the MOM and EOM-CCSD, while it appears significantly higher in energy compared to the 3 Σ state with the MOM.This point is confirmed by comparing the excitation energies from the 3 Π state for EOM-CCSD and the MOM, presented in the bottom left and right panels of Figure 7, respectively.The energy differences are slightly overestimated with the MOM for both the 1 Σ and 1 Δ states, confirming this trend.However, it cannot be argued that this result comes from the projection technique since both states are closed-shell; instead, it shows how the additional treatment of correlation could improve MOM calculations in strong magnetic fields.
As previously described, the electronic structure of molecules in strong magnetic fields can be highly sensitive to the orientation of the molecule with respect to the field, since the symmetry of the orbitals can be significantly lower in some orientations than in others.Figure 8 presents the total energy of the 1 Σ + /1 1 A′, 1 Π/1 1 A″, 1 Π/2 1 A′, and 1 Δ/3A′ states with respect to a magnetic field perpendicular to the molecular axis for FCI (top left) from ref 40 and the MOM (top right) in the cc-pVDZ basis set.The overall picture is qualitatively well represented by the MOM, especially for the lowest 1 Σ + /1 1 A′ state and the 1 Π/1 1 A″ state.For the 1 Π/2 1 A′ state, the behavior is slightly different with the MOM, particularly above 0.4B 0 , where the energy appears to change less with a further increase in field strength with the MOM than with FCI.This results in a crossing between 1 Π/2 1 A′ and 1 Π/1 1 A″ at about 0.9B 0 for FCI, while it occurs at slightly lower field strength for the MOM, at about 0.8B 0 .Finally, the 1 Δ/3A′ state computed with the MOM also presents a somewhat different shape, since the change in curvature at about 0.5B 0 in the FCI plot does not appear with the MOM.Despite these slight differences, the excitation energies from the 1 Σ + /1 1 A′ state presented in the bottom of Figure 8 show a consistent description with the MOM, despite the lack of correlation.We have also confirmed that MOM calculations with the substantially larger u-aug-cc-pCVQZ basis (not presented) yield results in close agreement with the top right panel of Figure 8.In particular, the crossing of the 1 Π/2 1 A′ and 1 Π/1 1 A″ states remains at ∼0.8B 0 and the 1 Δ/3A′ state also remains smooth around 0.5B 0 .

LIMITATIONS OF THE VARIANCE-BASED SCHEME
While the variance-based minimization has been shown to yield results consistent with those of energy-based optimization for some examples in Section 3.2, there are significant limitations that have been observed.Section 3 highlighted that the energies computed with the two methods may be different, while some states appear completely inaccessible through variance minimization.In addition, spurious solutions may be obtained using the variance minimization, which do not correspond to a physical solution. 70Although at first surprising, these limitations have been observed in the literature in different contexts both during the 60's 80,94 and recently. 84As discussed in Section 2, the variance function of eq 8 defines a sum of local minima, each of which has a minimum of zero.Each minimum corresponds to an eigenfunction of the Hamiltonian, and then, in contrast to the energy-based optimization that will have a (possibly degenerate) global minimum, all states are equally represented in terms of variance.However, in practical calculations, finite basis sets and single-determinant wave functions are used and may lead to minima with non-zero values, which vary according to how well the basis and wave function ansatz represent each particular state. 94.1.Difference between Energy-Based and Variance-Based Determinations.We now turn back to the example of the helium atom.In this example, the MOM and σ-SCF provide an equal description of the m s = ±1 determinants, while for some of the m s = 0 configurations, the difference in energy is significant.One of the properties of the variance is to describe how accurately expectation values are evaluated and provide upper and lower bounds to the eigenvalue of the ith state 94

⟨Φ| |Φ⟩ −
Here, one may readily see that smaller the variance is, the more accurate the wave function should be.Table 2 presents the energy difference between both methods with the associated variance for all configurations of the helium atom in the 6-311G basis at zero field.As one may see, the square root of the variance of the MOM and σ-SCF solutions is equal to zero for the m s = ±1 determinants, resulting in a similar evaluation of the energy with both approaches.For the m s = 0 configurations, [Φ] S is greater than zero, and the energy differences between the MOM and σ-SCF solutions are nonnegligible.Even if further trends are difficult to discern, one may conclude that since the m s = ±1 states seem to be well defined in this basis set, the variance-based and energy-based methods are mutually consistent.For the m s = 0 configurations, large values of [Φ] S yield a wider interval, in which the approximate energy of the ith eigenstate may be contained according to eq 15; there is then no reason for both approaches to give the same description.
4.2.Problem of the Masked States.In the example of the hydrogen molecule in Section 3, we commented that some states are inaccessible through variance-based minimization, which often converges to other solutions.It was observed that this can happen even for the ground state of a given symmetry, which would usually be considered easy to access via Δ-SCF approaches.Recently, Filippi et al. 84 reported the same problem in the context of QMC and pointed out that "this finding is unexpected, especially, considering that variance  minimization has been the method of choice in QMC for decades and is still routinely used".Noting that the minima may have non-zero values in finite basis calculations with single-determinant approaches, it may be possible that some higher lying minima could be masked 80 by overlapping variational wells associated with a lower minimum.In other words, some states are inaccessible through the variance because they are hidden by some others, and in this sense, this problem differs from the well-known variational collapse.
The consequence of the non-zero values of the minima may be illustrated by the functional defined in eq 11, which allows one to plot the profile of the variance with respect to ω for any state using its energy and variance; for states inaccessible through σ-SCF, the values of W[Φ] have been evaluated from the corresponding MOM solutions.Here, we focus on the triplet ground state 3 Σ u (blue) and the 3 Σ g (orange) and 3 Π g (yellow) excited states of the hydrogen molecule at zero field, and their profile has been plotted with respect to ω; the results are shown in the left and right panels of Figure 9 in the u-augcc-pVDZ and u-aug-cc-pVQZ basis sets, respectively.As shown in Section 3, the triplet ground state cannot be computed with σ-SCF at zero field in the u-aug-cc-pVQZ basis, in contrast to these two excited states.In the u-aug-cc-pVDZ basis (left panel of Figure 9), the orange and yellow minima are clearly defined, while the blue minimum is situated within the orange well.Despite this, all these states are accessible through a variance minimization.The overlap between these states can be calculated using the natural orbitals associated with the σ-SCF density matrix, revealing that the 3 Σ u and 3 Σ g states are orthogonal, while the 3 Σ u and 3 Π g states are non-orthogonal (with an overlap of 0.31).In the u-aug-cc-pVQZ basis (right panel), the situation is different, since the minimum of the blue well lies just within the yellow well.While the 3 Σ u state may be approached through the energy-targeting function in eq 11 for ω < −1.6 au, the subsequent variance minimization yields the 3 Π g state, which has a lower variance.Interestingly, the solution obtained by the energy-targeting optimization has an energy similar to that of the GS, to which it is expected to be a close approximation; however, its overlap with the 3 Π g state yielded by the subsequent variance minimization is 0.23.This indicates that a state is masked and inaccessible through a variance optimization if its minimum is included in the well of another state to which it is not orthogonal.
The use of W[Φ](ω) also provides an explanation of the fact that all states of interest of the hydrogen molecule are accessible in u-aug-cc-pVQZ in strong magnetic fields.Figure 10 presents the W[Φ](ω) functional for the lowest 3 Σ u , 3 Σ g , 3 Π g , and 3 Π u states at zero field as a solid line and in a magnetic field of 1B 0 parallel to the molecular axis with dashed lines.At zero field, all variance minima are concentrated in a small range of ω values, hiding the minima for 3 Σ u (blue) and 3 Π u (green).Applying a magnetic field results in a separation of the wells in ω, making all states accessible.Furthermore, even though the minimum of 3 Π u is included in the well of 3 Σ u in the presence of the field, both states are orthogonal, and so the 3 Π u state is not masked by the 3 Σ u state.
It is worth noting that while the 3 Σ u state cannot be optimized through the variance in u-aug-cc-pVQZ, it may be targeted through the minimization of the W[Φ](ω) function.Indeed, as shown in the right panel of Figure 9, the well of the 3 Σ u state corresponds to the minimum of W[Φ](ω) for ω < −1.6, and the energy-targeting step may be used to obtain an initial guess for this state.However, some states may be entirely inaccessible such as the 3 Π u , at zero field shown in Figure 10, since the entire green well is masked by the others.Furthermore, unlike in the context of energy optimization, improving the basis set does not necessarily improve the description of states in terms of variance. 94Indeed, the improvement of the basis set is not uniform in all states and may favor the description of some states compared to others.To show this, eq 10 may be expressed as In this equation, ⟨H ̂⟩ is evaluated with an approximate wave function ansatz, a single-determinant in σ-SCF, whereas E i is an FCI eigenvalue.As a result, when the basis set becomes more complete, the variance is not guaranteed to decrease and in some cases may actually increase if the state is poorly described by a determinant.A clear example is shown in Figure 9 with the 3 Σ u state where the variance, given by the minimum of the blue wells, is lower in u-aug-cc-pVDZ with a value of 0.010 than in u-aug-cc-pVQZ, where it has a value of 0.022.

CONCLUSIONS
This work has presented the first implementation and application of self-consistent field single-determinant methods to compute excited states in strong magnetic fields.In this context, two different approaches have been investigated: an energy-based scheme through the MOM/IMOM and a variance-based one with the σ-SCF method.Furthermore, we have presented the first development of the RI approximation for the σ-SCF method.
In this work, the RI σ-SCF implementation allowed for the use of larger basis sets than have been presented previously. 69hese calculations highlighted severe limitations of variancebased optimization methods, with many states becoming inaccessible.This surprising limitation was traced to the masking of states, in which the stationary point of W[Φ] occurs within the variational well of another state.In particular, if the states associated with these two overlapping wells are not orthogonal, then, the variance optimization will always converge to the state with lower variance.This behavior becomes more common in larger basis sets, as more states are accessible and many occur closer together in energy.
This effect was demonstrated for the simple case of the hydrogen molecule in its triplet state, where all low-lying states could be accessed in calculations using the u-aug-cc-pVDZ basis but many states became masked in the u-aug-cc-pVQZ basis.Unfortunately, these issues make the variance-based calculations currently unreliable for determining excitation energies.This was in addition to the previously reported appearance of spurious non-physical solutions in variancebased optimization, observed here in both basis sets. 70owever, while the variance optimization step suffers from the masking issue, the initial energy-targeting step still provides an interesting method to generate initial-guess wave functions; the adoption of this approach into other optimization methods will be the subject of future work.The observations from this work may also be of interest to the QMC community, where variance-based optimization is more frequently used.
In contrast, the MOM/IMOM has been widely used, and the present work shows further promising results for computing excited states in the presence of external magnetic fields.The MOM/IMOM was shown to be relatively reliable for computing excited states in the different examples treated in this work, and the results are in good agreement with those of EOM-CCSD and FCI reference calculations for the largest system, the methylidinium ion.Given the simplicity of the MOM/IMOM scheme and its more favorable computational cost, it provides an attractive approach to study excited states in strong magnetic fields.
The present work has only used the MOM with the HF energy functional, to facilitate a direct comparison with the σ-SCF method.Accounting for correlation, either through post-HF approaches or using current-DFT (CDFT), would be straightforward and should provide a better quantitative agreement with FCI/EOM-CC calculations.The inclusion of correlation may be particularly important for the computation of molecular properties, and further work will be carried out in this direction.For example, work using a recent implementation of geometrical gradients for LAO integrals 95 in combination with the MOM is currently underway to allow for the determination of equilibrium molecular geometries of excited states under strong magnetic fields.

■ APPENDIX σ-SCF Equations and RI Approximations
For the sake of brevity, the working equations are not presented in Section 2 and are given in the present Appendix.Note that in the following, an orthogonal basis set of AOs is assumed.The σ-SCF method is based on two successive optimizations: the energy-targeting step eq 11 and the optimization of the variance eq 8. Since both functions are related by only the form of S[Φ] is given in this Appendix.A key quantity in the σ-SCF method for optimizing S[Φ] is obtained by differentiating eq 8 with respect to the spin density matrix to yield   (18)   where is referred to as a modified Fock matrix in ref 69, the construction of which is simplified here by use of the intermediate RI σ-SCF Equations The RI approximation has been widely used in electronic structure theory and has notably been extended to the strongmagnetic field regime to compute the energy 96,97 and geometrical derivatives. 95As discussed in the present work, the computation of the square of the Hamiltonian and the resulting 8-index summations make the variance-based scheme an expensive approach.The RI approximation, in which twocenter charge distributions are expanded in an atom-centered auxiliary basis, has been applied in this work in order to speed up the σ-SCF method.Four-center integrals are approximated by a contraction of three-center integrals with the inverse of the Coulomb metric of the auxiliary functions 98 where ϕ p represents the orbital basis functions (GAOs or LAOs) and χ U represents the atom-centered auxiliary basis functions (always GAOs), and the Coulomb metric is defined as From these definitions, the two-electron integrals may now be approximated as where The variance may then be expressed as Here, the intermediates Following the same derivation, the modified Fock matrix of eq 18 may be expressed as   (35)   and

Numerical Validation
Table 3 presents numerical applications of the RI equations developed in this work on the singlet ground state of the hydrogen molecule for three casesin the absence of a magnetic field and with a magnetic field of 1B 0 parallel and perpendicular to the bond axis.The study is carried out using the u-aug-cc-pVDZ, u-aug-cc-pVTZ, and u-aug-cc-pVQZ orbital basis, while the size of the auxiliary basis for the RI expansion is increased, from aug-cc-pVDZ-RI to aug-cc-pV6Z-RI, and compared to the standard (non-RI) calculation for each basis.A comparison with the AutoAux auxiliary basis 88 used in this work is also presented for informative purposes.It is clear that the RI calculations approach the non-RI values as the size of the auxiliary basis increases, with the difference decreasing to about 10 −8 E h with u-aug-cc-pVDZ/aug-cc-pV6Z-RI, confirming the accuracy of the present implementation.The AutoAux auxiliary basis used throughout this work offers accuracy, for a given aug-cc-pVNZ orbital basis, comparable to that using an aug-cc-pV(N + 2)Z-RI basis.

Figure 2 .
Figure 2. Total energy (E h ) computed with σ-SCF and the MOM as a function of magnetic field strength in B 0 for the m s = ±1 configurations of He in the 6-311G basis set.

Figure 3 .
Figure 3.Total energies (left panel) computed with σ-SCF and the MOM and energy differences between both (right panel) for the m s = 0 configurations of He in the 6-311G basis.

a
Note that the FCI ⟨S 2⟩ is given here to indicate the states considered, and MOM and σ-SCF yield values close to 1 for the openshell singlet statesreflecting significant spin contamination.Blank entries correspond to doubly excited states that are not accessible from CIS calculations.b The CIS calculation is based on the 1s β 2s β SCF solution.c The CIS calculation is based on the 1s α 1s β SCF solution.d Energies of the open-shell singlets are obtained using the projection in eq 14.

Figure 4 .
Figure 4. Energy spectrum (E h ) in the u-aug-cc-pVDZ/AutoAux basis (left panel) and in the u-aug-cc-pVQZ/AutoAux basis (right panel) of the m s = −1 components computed through CIS, the MOM/IMOM, and σ-SCF in the absence and in the presence of a magnetic field of 1B 0 parallel and perpendicular to the bond axis of the hydrogen molecule.Gray lines correspond to unphysical σ-SCF solutions.

Figure 5 .
Figure 5.Total energies (E h ) computed with EOM-CCSD (top left) and the MOM (top right) for the hydrogen molecule in the u-aug-cc-pVQZ/ AutoAux basis in a magnetic field parallel to the molecular axis and their energy difference with the 1 Σ g + / 1 Σ g state (bottom left and right).Data for the EOM-CCSD calculations are taken from ref 39.

Figure 6 .
Figure 6.Total energies computed with EOM-CCSD (top left) and the MOM (top right) for the hydrogen molecule in the u-aug-cc-pVQZ/ AutoAux basis in a magnetic field perpendicular to the molecular axis and their energy difference with the 1 Σ g + / 1 A g state (bottom).Data for the EOM-CCSD calculations are taken from ref 39.

Figure 7 .
Figure 7.Total energies computed with EOM-CCSD (top left) and the MOM (top right) for CH + in the u-aug-cc-pCVQZ/AutoAux basis in a magnetic field parallel to the molecular axis and their energy difference with the 3 Π (bottom).Data for the EOM-CCSD calculations are taken from ref 39.

Figure 8 .
Figure 8.Total energies computed with FCI (top left) and the MOM (top right) for CH + in the cc-pVDZ/AutoAux basis in a magnetic field perpendicular to the molecular axis and their energy difference with the 1 Σ + /1 1 A′ (bottom).Data for the FCI calculations are taken from ref40.

Figure 9 .
Figure 9. W[Φ] with respect to ω for the 3 Σ u , 3 Σ g , and 3 Π g states of the hydrogen molecule (r H−H = 1.4 bohr) in the absence of a magnetic field in the u-aug-cc-pVDZ/AutoAux basis (left) and the u-aug-cc-pVQZ/AutoAux basis (right).

Figure 10 .
Figure 10.W[Φ](ω) with respect to ω for the 3 Σ u , 3 Σ g ,3 Π u , and 3 Π g states of the hydrogen molecule in the absence of a magnetic field (with r H−H = 1.40 bohr) in solid lines and in the presence of a parallel magnetic field (with r H−H = 1.25 bohr) in dashed lines in the u-aug-cc-pVQZ/AutoAux basis.
jk lq D Q D rs iq lk pj D Q D rs iq pk lj D Q Q rs pi lk jq D Q D rs iq lk pj ( )

1 ( 20 )
Journal of Chemical Theory and Computationwith U and V being the atom-centered auxiliary basis and the three-center integrals, respectively, been introduced, which are expressed in terms of the occupied orbitals m and n through the following decom- Corresponding Author Andrew M. Teale − School of Chemistry, University of Nottingham, Nottingham NG7 2RD, U.K.; Hylleraas Centre for Quantum Molecular Sciences, Department of Chemistry, University of Oslo, Oslo N-0315, Norway; orcid.org/0000-0001-9617-1143; Email: andrew.teale@nottingham.ac.uk

Table 1 .
Mean Absolute Error (mE h ) for the MOM, σ-SCF, and CIS Relative to FCI Calculations for Each State of the Helium Atom in the 6-311G Basis over a Range of Magnetic Fields from 0 to 1B 0 with a Step of 0.05B 0

Table 2 .
Energy Difference (E h ) between the MOM and σ-SCF and the Square Root of the Variance for all Determinants of the Helium Atom in the 6-311G Basis in the Absence of a Magnetic Field )

Table 3 .
Total Energy of the Hydrogen Molecule (r H−H = 1.4 bohr) Computed with σ-SCF