Accurate Computation of Quantum Excited States with Neural Networks
We present a variational Monte Carlo algorithm for estimating the lowest
excited states of a quantum system which is a natural generalization of the
estimation of ground states. The method has no free parameters and requires no
explicit orthogonalization of the different states, instead transforming the
problem of finding excited states of a given system into that of finding the
ground state of an expanded system. Expected values of arbitrary observables
can be calculated, including off-diagonal expectations between different states
such as the transition dipole moment. Although the method is entirely general,
it works particularly well in conjunction with recent work on using neural
networks as variational Ans\"atze for many-electron systems, and we show that
by combining this method with the FermiNet and Psiformer Ans\"atze we can
accurately recover vertical excitation energies and oscillator strengths on a
range of molecules. Our method is the first deep learning approach to achieve
accurate vertical excitation energies, including challenging double
excitations, on benzene-scale molecules. Beyond the chemistry examples here, we
expect this technique will be of great interest for applications to atomic,
nuclear and condensed matter physics.
David Pfau, Simon Axelrod, Halvard Sutterud, Ingrid von Glehn, James S. Spencer
# Natural Quantum Monte Carlo Computation of Excited States

###### Abstract

We present a variational Monte Carlo algorithm for estimating the lowest excited states of a quantum system which is a natural generalization of the estimation of ground states. The method has no free parameters and requires no explicit orthogonalization of the different states, instead transforming the problem of finding excited states of a given system into that of finding the ground state of an expanded system. Expected values of arbitrary observables can be calculated, including off-diagonal expectations between different states such as the transition dipole moment. Although the method is entirely general, it works particularly well in conjunction with recent work on using neural networks as variational Ansatze for many-electron systems, and we show that by combining this method with the FermiNet and Posiformer Ansatze we can accurately recover vertical excitation energies and oscillator strengths on molecules as large as benzene. Beyond the examples on molecules presented here, we expect this technique will be of great interest for applications of variational quantum Monte Carlo to atomic, nuclear and condensed matter physics.

## I Introduction

The computation of excited states properties of quantum systems is a fundamental challenge in chemistry and many branches of physics. Understanding electronic excitations is critical for predicting photochemical phenomena such as fluorescence and conformational changes in the presence of light [1; 2]. In condensed matter physics, excitations determine the optical band gap of semiconductors, which is critical for predicting the behavior of solar cells, photosensors, LEDs and lasers [3]. Excited states are also relevant to understanding nuclear phenomena like metastable isomers and electron capture [4]. Ultimately, the dynamics of quantum systems when stimulated cannot be understood without taking excited states into account. Despite the importance of excited states for quantum phenomena, a full computational account of excited states remains challenging.

Quantum Monte Carlo (QMC) methods [5; 6] are an appealing class of algorithms for computing the behavior of quantum systems due do the favorable scaling with the number of particles, typically \(\mathcal{O}(N^{3})-\mathcal{O}(N^{4})\), and wide applicability. Variational quantum Monte Carlo (VMC) in particular is quite conceptually simple, and consists of finding an explicit functional form for a wavefunction which minimizes a variational bound, but historically was not considered accurate enough on its own for many demanding applications. Recent work using neural networks as a wavefunction Ansatz has reinvigorated interest in VMC [7; 8], and has demonstrated that VMC can be competitive with state-of-the-art methods for ground state calculations.

In this paper, we focus on computing excited states of quantum systems by VMC. When used to optimize ground states, there are only two variational principles for QMC - energy minimization and variance minimization. Innovations in ground state VMC primarily focus on the choice of trial wavefunction [9; 10], or optimization method used to achieve the variational bound [11; 12], but the choice of objective to optimize is well-established. The same cannot be said for variational optimization of excited states.

Approaches for computing excited states by VMC can be broken down into several categories. Most methods are either state-_targeting_, in that they aim to find a single excited state, or state-_averaging_, in that they aim to find the lowest-lying exciting states by minimizing the total weighted energy of many states simultaneously. Among state-targeting methods, there are methods which target specific energy ranges [13; 14], specific symmetries of the system [15], or a specific ordering of the roots (i.e. the \(k\)-th lowest state) [16]. For state-averaging approaches, the different states must be kept orthogonal, which can be achieved by including a penalty term in the variational bound which pushes the states apart [17; 15; 18], or by explicitly constructing orthogonal Ansatze, sometimes repeatedly re-orthogonalizing during optimization [19; 20; 21; 22].

All of these approaches have drawbacks and limitations. Targeting specific symmetries or energy ranges requires prior knowledge about the states of interest which may not be available, and state-targeting by variance minimization can lose track of the desired state [21]. Root-targeting methods are prone to root-flipping, whether they are used for QMC or other computational paradigms [23; 24]. Some methods require solving a generalized eigenvalue problem from stochastic estimates of the Hamiltonian and overlap matrices, which introduces biases into the gradients [16; 25]. Penalty methods often have problems with multiple Ansatze collapsing onto the same state, or have biased gradients [18], and the
strength of the penalty term is a free parameter which must be chosen. Constructing orthogonal Ansatze is usually only possible when the Ansatz is a linear combination of basis set functions [26; 27], which rules out many recently-developed Ansatze based on deep neural networks [28; 29; 30; 7]. Heuristics such as variance matching may be required to achieve good numerical results for all approaches. Despite almost four decades of work on QMC methods for excited states [26; 31], no single variational principle has emerged which has no free parameters, has convergence guarantees when optimizing with noisy Monte Carlo estimates, and is applicable to all possible Ansatze and all excited states, regardless of symmetry.

Here we present a new variational principle for computing the lowest excited states of a quantum system by Monte Carlo which does not suffer from any of these limitations. Our method can be seen as a state-averaging approach with a particular choice of sampling distribution which does not require the states to be orthogonal. This choice of sampling distribution is equivalent to reformulating the problem of finding \(K\) excited states of an \(N\) particle system into the problem of finding the ground state of a \(K\)-fermion system where each fermion is equivalent to \(N\) particles in the original system. Instead of orthogonalizing the states, the local energy is promoted from a scalar to a matrix, which gives unbiased estimates of a matrix whose eigenvalues are the energies of orthogonal states. Because wavefunction optimization can be done by stochastic gradient descent from unbiased noisy estimates of the total energy, the procedure is guaranteed to converge to a local minimum of the total energy over states. Due to the many desirable mathematical properties which follow from the choice of sampling distribution, we refer to our proposed approach as _natural excited states_ for VMC (NES-VMC).

## II Method

### Variational Monte Carlo

First we briefly review ground-state VMC and establish some notation. We will stick to the notation of first quantization and consider a system of \(N\) particles with states \(\mathbf{x}=\mathbf{x}_{1},\ldots,\mathbf{x}_{N}\), although everything we discuss could be applied to variational Ansatze represented in second quantization as well. We aim to find the lowest eigenfunction of a Hamiltonian operator \(\hat{H}\). This can be done by reformulating the eigenfunction problem in variational form, as one of finding the minimum of the Rayleigh quotient:

\[\psi^{*}=\arg\min_{\psi}\frac{\langle\psi\hat{H}\psi\rangle}{\langle\psi^{2} \rangle} \tag{1}\]

where the Ansatz \(\psi\) is not necessarily normalized. Computing this quotient involves taking high-dimensional integrals over all possible particle states \(\mathbf{x}\), and can be approximated by Monte Carlo integration. Many choices of Monte Carlo sampling distribution \(p(\mathbf{x})\) are possible, but if \(p(\mathbf{x})\propto\psi^{2}(\mathbf{x})\), then the Rayleigh quotient take a simple form that allows for unbiased empirical estimation of the energy and gradients of the energy:

\[\frac{\langle\psi\hat{H}\psi\rangle}{\langle\psi^{2}\rangle}=\mathbb{E}_{ \mathbf{x}\sim\psi^{2}}\left[\psi^{-1}(\mathbf{x})\hat{H}\psi(\mathbf{x})\right] \tag{2}\]

For this reason, \(\psi^{2}\) is the natural choice of sampling distribution for ground state estimation. The scalar \(E_{L}(\mathbf{x})\mathbf{\triangleq}\psi^{-1}(\mathbf{x})\hat{H}\psi(\mathbf{ x})\) that appears inside the expectation is the _local energy_, and at any eigenfunction of \(\hat{H}\) it will be constant if \(\hat{H}\) is a local operator.

### Natural Excited States

Going from ground states to excited states, we aim to find the lowest \(K\) eigenfunctions of \(\hat{H}\). We refer to a single set of \(N\) particle states as a _particle set_, and denote different particle sets with an upper index, so that \(\mathbf{x}^{i}\) denotes a set of \(N\) particles \(\mathbf{x}^{i}_{1},\ldots,\mathbf{x}^{i}_{N}\). For the remainder of the article, we will use \(\mathbf{x}\) to denote the complete state of all particle sets \(\mathbf{x}^{1},\ldots,\mathbf{x}^{K}\). Let \(\psi_{i}\) denote a (possibly unnormalized) N-particle wavefunction, then we are trying to find wavefunctions \(\psi_{1},\ldots,\psi_{K}\) which approximate the lowest excited states. Let \(\mathbf{\Psi}(\mathbf{x})\in\mathbb{R}^{K\times K}\) denote the matrix combining all electron sets with all wavefunctions:

\[\mathbf{\Psi}(\mathbf{x})\overset{\triangle}{\equiv}\begin{pmatrix}\psi_{1}( \mathbf{x}^{1})&\ldots&\psi_{K}(\mathbf{x}^{1})\\ \vdots&&\vdots\\ \psi_{1}(\mathbf{x}^{K})&\ldots&\psi_{K}(\mathbf{x}^{K})\end{pmatrix} \tag{3}\]

The determinant of this matrix \(\Psi(\mathbf{x})=\det(\mathbf{\Psi}(\mathbf{x}))\) can be thought of as an unnormalized Slater determinant, except that instead of single-particle orbitals, it is made up of N-particle wavefunctions. We call \(\Psi(\mathbf{x})=\det(\mathbf{\Psi}(\mathbf{x}))\) the _total Ansatz_, while the individual \(\psi_{i}\) are the _single-state Ansatze_.

Rather than optimizing the single-state Ansatze in order from lowest to highest energy, we will only optimize the total Ansatz to minimize the total energy of all states. This is conceptually quite similar to state-averaging approaches in VMC, except that we will not explicitly enforce the orthogonality of the different single-state Ansatze. Note that taking any linear combination of single-state Ansatze \(\psi^{\prime}{}_{i}=\sum_{j}a_{ij}\psi_{j}\) only changes the total Ansatz by a constant factor. Also note that if two single-state Ansatze are the same, the total Ansatz becomes zero. Thus, by representing the total Ansatz as a determinant of single-state Ansatze, we can prevent the collapse of different Ansatze onto the same state, without requiring them to be orthogonal.

For an arbitrary operator \(\hat{\mathcal{O}}\) that acts on \(N\)-particle wavefunctions, let \(\hat{\mathcal{O}}\mathbf{\Psi}(\mathbf{x})\) denote the matrix of all values
of this operator applied to all single-state Ansatze and particle sets:

\[\hat{\mathcal{O}}\boldsymbol{\Psi}(\mathbf{x})\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!
Not only is diagonalizing \(\mathbb{E}_{\Psi^{2}}[\mathbf{E}_{L}(\mathbf{x})]\) sufficient to recover the energies - it also provides us with the necessary change of basis to evaluate other observables \(\hat{\mathcal{O}}\), even off-diagonal observables \(\langle\psi_{i}\hat{\mathcal{O}}\psi_{j}\rangle\) between states. This can be seen due to the identity \(\mathbb{E}_{\Psi^{2}}[\mathbf{\Psi}^{-1}\hat{\mathcal{O}}\mathbf{\Psi}]= \mathbf{S}^{-1}\hat{\mathbf{O}}\), and for single-state Ansatze which are a linear combination of eigenfunctions, \(\mathbf{S}^{-1}\hat{\mathbf{O}}=\mathbf{A}^{-1}\hat{\mathbf{O}}^{*}\mathbf{A}\). So if we accumulate and diagonalize \(\mathbb{E}_{\Psi^{2}}[\mathbf{E}_{L}(\mathbf{x})]\) and use the resulting eigenvectors to compute \(\mathbf{U}^{-1}\mathbb{E}_{\Psi^{2}}[\mathbf{\Psi}^{-1}\hat{\mathcal{O}} \mathbf{\Psi}]\mathbf{U}\), then in the vicinity of the true ground state of the total Ansatz the result will be approximately \(\mathbf{\Sigma}^{-1}\hat{\mathbf{O}}^{*}\mathbf{\Sigma}\). Along the diagonal, this gives exactly the expectations \(\langle\psi_{i}^{*}\hat{\mathcal{O}}\psi_{i}^{*}\rangle\). Off the diagonal, this yields \(\frac{\sigma_{i}}{\sigma_{j}}\langle\psi_{i}^{*}\hat{\mathcal{O}}\psi_{j}^{*}\rangle\). If we multiply the matrix elementwise by its transpose, the \(\sigma_{i}\) terms cancel out, and we recover \(\langle\psi_{i}^{*}\hat{\mathcal{O}}\psi_{j}^{*}\rangle^{2}\), which gives the expectation up to a sign factor. This sign factor is not physically observable however, and in practice for computing quantities like the oscillator strength, only the expectation squared is needed.

### Neural Network Ansatze

The use of variational Monte Carlo for ground state calculations was typically used to find a trial wavefunction for more accurate projector QMC methods like diffusion Monte Carlo [5] or auxiliary field Monte Carlo [33]. However, in recent years, advances in deep neural networks have led to their use as accurate Ansatze for studying spin systems [7], electronic structure [8] and nuclear systems [34], often reaching levels of accuracy rivaling projector QMC methods. This has led to a renewed interest in VMC as a standalone method. While a variety of different neural network architectures can be used depending on the problem, such as restricted Boltzmann machines [7], convolutional neural networks [35], and autoregressive models [36], a number of custom architectures have been developed specifically for many-body electronic structure problems in first quantization [28; 29; 30; 37; 38; 39; 40; 41]. Most of these Ansatze start from a linear combination of Slater determinants:

\[\psi(\mathbf{x})=\sum_{k}\omega_{k}\text{det}\begin{pmatrix}\phi_{1}^{k}( \mathbf{x}_{1})&\ldots&\phi_{N}^{k}(\mathbf{x}_{1})\\ \vdots&&\vdots\\ \phi_{1}^{k}(\mathbf{x}_{N})&\ldots&\phi_{N}^{k}(\mathbf{x}_{N})\end{pmatrix} \tag{11}\]

It has long been recognized [42] that the single-particle orbitals in a Slater determinant can be generalized to depend on _all_ particles, so long as they depend on all but one in a permutation-independent manner:

\[\psi(\mathbf{x})=\sum_{k}\omega_{k}\text{det}\begin{pmatrix}\phi_{1}^{k}( \mathbf{x}_{1};\{\mathbf{x}_{/1}\})&\ldots&\phi_{N}^{k}(\mathbf{x}_{1};\{ \mathbf{x}_{/1}\})\\ \vdots&&\vdots\\ \phi_{1}^{k}(\mathbf{x}_{N};\{\mathbf{x}_{/N}\})&\ldots&\phi_{N}^{k}(\mathbf{x }_{N};\{\mathbf{x}_{/N}\})\end{pmatrix} \tag{12}\]

where \(\{\mathbf{x}_{/i}\}\) denotes the set of all particles _except_\(\mathbf{x}_{i}\). In the event that the particles are spin-assigned, the orbitals can also be expressed as \(\phi_{i}^{k}(\mathbf{x}_{j}^{\dagger};\{\mathbf{x}_{/j}^{\dagger}\},\{\mathbf{ x}^{\dagger}\})\) where the function is only invariant to changing the order of particles of the same spin. Most neural network Ansatze for electrons in real space implement this idea by using permutation-equivariant deep neural networks to represent the orbitals, sometimes with a multiplicative Jastrow factor to account for pairwise interactions [29; 30; 38].

Extending these Ansatze to represent multiple states is quite straightforward. Each state is still expressed as a sum of determinants of generalized neural network orbitals, there are simply more orbitals:

\[\psi_{i}(\mathbf{x})=\sum_{ik}\omega_{ik}\text{det}\begin{pmatrix}\phi_{1}^{ ik}(\mathbf{x}_{1};\{\mathbf{x}_{/1}\})&\ldots&\phi_{N}^{ik}(\mathbf{x}_{1};\{ \mathbf{x}_{/1}\})\\ \vdots&&\vdots\\ \phi_{1}^{ik}(\mathbf{x}_{N};\{\mathbf{x}_{/N}\})&\ldots&\phi_{N}^{ik}(\mathbf{x }_{N};\{\mathbf{x}_{/N}\})\end{pmatrix} \tag{13}\]

Nothing is changed about the neural network architecture itself, just the number of orbitals is increased proportionally to the number of states.

Neural network Ansatze differ from classic Ansatze like the Slater-Jastrow-backflow Ansatz [9] in important ways which make it difficult to apply existing excited state methods. Many methods assume that the Ansatz is a linear combination of orthogonal basis functions like Slater determinants, a necessary assumption for maintaining the orthogonality of states, either through explicit construction or a diagonalization step [19]. Classic Ansatze are usually optimized through a small number of gradient steps, where each gradient step is accumulated over a large number of MCMC steps, so that the gradients are nearly deterministic. Most modern deep neural networks, by constrast, are optimized by stochastic gradient descent using a large number of small, noisy steps [43]. This means bias in the gradients becomes a more significant concern.

Existing work on excited state calculations with neural networks has focused on penalty methods [18; 15], but these still require choosing a free parameter trading off total energy and penalty strength, and may not exactly satisfy orthogonality in the states. Some of these methods also have biased gradients in the penalty term [18] due to nonlinearities meant to push states apart more strongly. By contrast, the NES-VMC method has no free parameters to tune, can be optimized by unbiased gradients that have the same form as for ground state calculations, does not require the states to be orthogonal, and makes no assumption on the functional form of the Ansatz. Thus, while NES-VMC is generally applicable to _all_ excited state VMC calculations, it is particularly well-tailored for use with recently developed neural network Ansatze.

## III Results

While the natural excited states method is fully general and can be applied to any quantum Hamiltonian, our experimental validation is focused on electronic structure in atoms and molecules, due to the abundant experimental and computational literature to compare against. For
all experiments, we are solving the Schrodinger equation in the Born-Oppenheimer approximation [44]:

\[\hat{H}= -\frac{1}{2}\sum_{i}\nabla_{i}^{2}+\sum_{i>j}\frac{1}{|\mathbf{r}_ {i}-\mathbf{r}_{j}|}\] \[-\sum_{iI}\frac{Z_{I}}{|\mathbf{r}_{i}-\mathbf{R}_{I}|}+\sum_{I>J }\frac{Z_{I}Z_{J}}{|\mathbf{R}_{I}-\mathbf{R}_{J}|} \tag{14}\]

where the indices \(i\) and \(j\) are over electrons and \(I\) and \(J\) are over atomic nuclei with fixed locations.

To try to disentangle the effect that the choice of Ansatz has on performance, we investigated two different neural network architectures: the Fermionic Neural Network (FermiNet) [29] and the Wavefunction Transformer (Psiformer) [38]. While the Psiformer has generally been found to be more accurate on large systems, it is also slower, and for ground state calculations up to approximately 15 electrons, no appreciable difference in accuracy between the two has been found.

### Atomic Spectra

As an initial check of the correctness of our method, we investigate the excited states of first-row atoms, from lithium to neon. Atomic spectral lines have been the subject of some of the highest-precision measurements in all of science, and while we do not aim to reach spectroscopic accuracy, we can have high confidence in accuracy of the measurements, and do not need to worry about effects such as adiabatic relaxation and zero-point vibrational energy which affect molecular measurements. All experimental data was taken from the energy level tables in the NIST Handbook of Basic Atomic Spectroscopic Data [32]. Because we are working with the nonrelativistic Schrodinger equation without spin-orbit corrections, we are not able to compute fine or hyperfine structure. To remove the fine structure, experimental energy levels with different total angular momenta are averaged together weighted by the degeneracy \(m_{J}=2J+1\) and treated as a single level. The hyperfine structure is too small to be of concern here. To investigate the effect of the choice of Ansatz as well as the choice of number of states \(k\) to compute, we ran calculations with the FermiNet with both 5 and 10 states, as well as the Psiformer with 10 states. Results are given in Fig. 1, with numerical results (including error bars) in the Appendix in Table 2.

For all atoms, NES-VMC gives results closely matching experiment. From lithium up to oxygen, the error relative to experiment is far less than 1 mHa (27.2 meV) for all but the highest excited state, and is often less than 0.1 mHa, an exceedingly high level of accuracy for a deep neural network Ansatz. On lithium, all Ansatze correctly converge to the \({}^{2}S\) and \({}^{2}P^{\circ}\) states, which are missed by the PauliNet penalty method [18]. The method struggles in some cases to get the highest energy state correct, but this seems to be improved by simply computing more states - for instance, the error in the \({}^{4}P\) states of fluorine is cut in half by increasing the number of states from 5 to 10. In rare cases, the highest state

Figure 1: Excited state energies for first row atoms from lithium to neon. Results from natural excited state VMC applied to the FermiNet (10 states, blue, 5 states, red) are shown on top of experimental results [32]. Spectral lines which match computed states are labeled with electron configurations and atomic term symbols (except for the highest levels of F and Ne, where term symbols are omitted for clarity). For all but the largest systems and highest excited states, there is excellent agreement with experiment. The discrepancy between 5 and 10 excited states is minimal except for the highest excited states of F and Ne, where computing more states increases the accuracy of a given state. Complete numerical results are given in Table 2.
seems to converge to the incorrect state, such as boron with the Psiformer, which seems to converge to the \({}^{2}P^{\circ}\) state rather than the last \({}^{2}D\) state. Fluorine and neon both have relatively large errors on the order of 1-2 mHa for low-lying states, but going from the FermiNet to the more accurate Psiformer Ansatz seems to reduce this error in all cases. The largest errors are in the highest states of fluorine and neon, where the error is significant. In this case we suspect the difficulty is due to the large number of different states with similar electron configurations and energies, and hope that by computing even more states or by using even more expressive Ansatze, the effects of individual states can be disentangled. The excellent performance on low-lying states gives us confidence that NES-VMC is mathematically sound.

### Oscillator Strengths

Going beyond results on single atoms and vertical excitation energies, we are interested in the performance of NES-VMC on more complicated molecular systems, as well as observable quantities other than the energy. The QUEST database [45; 46; 48; 49; 50; 51; 52; 53; 54] is an excellent source of well-controlled benchmark vertical excited states calculations using coupled cluster methods on molecules of various sizes, with consistent geometries and basis set extrapolations. Of particular interest is the subset of QUEST for which oscillator strengths have been computed [45], as oscillator strengths provide a strong test of how well an excited state method can perform on experimentally-observable quantities, and especially as oscillator strength and transition probability calculations are known to be highly sensitive to choices of basis set [55].

Oscillator strengths are a measure of the probability of transition between different states occurring as a result of photon emission or absorption. Under the assumption that the wavelength of the photon is much longer than the system under consideration, so the interaction can be approximated by a constant electric field, the transition dipole moment between two states gives a measure of how that transition will interact with light:

\[\mathbf{d}_{ij}=\left\langle\psi_{i}^{\dagger}\sum_{k}q_{k}\mathbf{r}_{k}\psi_ {j}\right\rangle \tag{15}\]

where the sum over \(k\) is taken over all particles in the system with charge \(q_{k}\) and position \(\mathbf{r}_{k}\). For electrons, \(q_{k}=-e\). The transition dipole moments are vector-valued quantities which include a complex phase factor, and are not directly observable. The oscillator strength

Figure 2: Vertical excitation energies and oscillator strengths for small molecules from Chrayteh _et al._[45]. Singlet states are in blue and triplet states are in gray. NES-VMC results are indicated by markers while theoretical best estimates from Chrayteh _et al._[45] or directly from QUEST [46] are given by the lines. When no data from QUEST is available, no TBE is given. Experimental results from Chrayteh _et al._[45] and references thererin are given by the dashed lines in green. Where available, energies and oscillator strengths from Entwistle _et al._[18] are provided by the black triangles for comparison, with (pointing left) and without (pointing right) variance matching. In most cases, our results on both energies and oscillator strengths agree closely with theoretical best estimates. Complete numerical results are given in Table 3.
of a particular transition can be computed from the transition dipole moment:

\[f_{ij}=\frac{2}{3}\frac{m}{\hbar^{2}}\left(E_{i}-E_{j}\right)|\mathbf{d}_{ij}|^{2} \tag{16}\]

which reduces the transition dipole moment to a dimensionless positive scalar. In natural excited states, we can compute expectations of operators between different states up to an arbitrary sign factor, and that sign factor goes away in the oscillator strength. Computational details are discussed in more detail in Sec. C.3.

We applied NES-VMC to all of the small molecules investigated in Chrayteh _et al._[45], computing the 5 lowest energy states with both the FermiNet and Psiformer. Results are presented in Fig. 2 and Table 3. Wherever possible, we take results from QUEST [45, 46] to be theoretical best estimates (TBEs) for comparison, though for many of the states we converged to, especially triplets, no results exist in QUEST. For molecules with heavier atoms (HCl, H\({}_{2}\)S, H\({}_{2}\)CSi), we found that using pseudopotentials for the heaviest atoms significantly improved the accuracy of the results, likely because the total energy scale was reduced by ignoring core electrons. Where applicable, we also include a comparison against the VMC penalty method of Entwistle _et al._[18]. We omit N\({}_{2}\) because the lowest-lying excited states are all triplets. For all diatomic systems, the \({}^{1}\Pi\) state is doubly-degenerate, and so the baseline oscillator strengths are divided by two to match the computed results.

In almost all cases, both the vertical excitation energies and the oscillator strengths are in excellent agreement with the TBE. The vertical excitation energies are almost all within chemical accuracy (1.6 mHa or 0.04 eV) of the TBE while the oscillators strengths usually diverge from the TBE by at most an amount on the order of 0.001, comparable to the uncertainty in the calculations. The results of Entwistle _et al._, in contrast, often differ noticeably from other theoretical results, even when correction using variance matching are applied. This is particularly noticeable for the oscillator strengths. We note that we do not use variance matching for any of the NES-VMC calculations.

There are a few cases where NES-VMC behaves oddly. While the FermiNet and Psiformer find nearly identical vertical excitation energies for the \({}^{1}\Pi\) state of HCl, and the FermiNet accurately predicts the oscillator strength, the Psiformer mistakenly finds this to be a dark state. On formaldehyde (CH\({}_{2}\)O), both the FermiNet and Psiformer fail to find the \({}^{3}A_{1}\) state at all, and the oscillator strength for the \({}^{1}B_{2}\) state diverges from the TBE by a significant margin, although the Psiformer halves that margin relative to the FermiNet. Vertical excitation energies for systems with heavier atoms, such as H\({}_{2}\)S, and the highest state of thioformaldehydehyde (CH\({}_{2}\)S), are not quite

Figure 3: Excited states of the carbon dimer (C\({}_{2}\)). (a) The symmetries of the different states can be identified by evaluating each single state Ansatz at location \(\mathbf{r}\) and \(-\mathbf{r}\) for parity symmetry (u/g, blue) or by flipping \(\mathbf{r}\) across the x-axis for reflection symmetry (+/–, orange). (b) The vertical and adiabatic energies of excited states of C\({}_{2}\). The green line indicates experimental energies [47] and the red line indicates the energy of the \(B^{1}\Delta_{g}\) state from QUEST [48]. Bright transitions are labelled with their oscillator strength and, when available, their names. (c) Visualization of the 8 lowest natural orbitals of C\({}_{2}\). (d) The occupancy of the different natural orbitals for the different excited states of C\({}_{2}\), identified from the density matrix of each state. The \(a^{3}\Pi_{u}\) through \(A^{1}\Pi_{u}\) states are single excitations while the last two states are double excitations. Complete numerical results are given in Table 4.
as accurate as other results, though in the case of thioformaldehyde we are hopeful that, consistent with the atomic results in the previous section, computing more states will reduce the error in the \({}^{3}B_{2}\) state. For nitroxyl (HNO), the FermiNet fails to converge to the \({}^{1}A^{\prime}\) state, but the Psiformer finds it correctly, albeit with a relatively large error in the vertical excitation energy. This suggests that there are occasional difficulties in getting NES-VMC to converge to all low-lying states, but we are hopeful that improvements in optimization methods can improve this in the future. What is clear is that NES-VMC works well in the large majority of cases, and is far more accurate than alternative methods which have been proposed for neural network Ansatze.

Other QMC methods have also been applied to some of these systems. In particular, the QMC-CIPSI method has been successfully applied to computing the vertical excitation energies of the \({}^{1}A_{2}\) state in formaldehyde and thioformaldehyde to within chemical accuracy, using a conventional Slater-Jastrow Ansatz [56]. While the QMC-CIPSI method cannot be applied to neural network Ansatze, this suggests that good results can still be achieved with VMC with a simple Ansatz, and that the benefit of using NES-VMC relative to the penalty method in Entwistle _et al._ is due to the method rather than the choice of Ansatz.

### Carbon Dimer

In addition to computing observable quantities, it is also desirable to be able to say something about the _nature_ of different excited states - whether a state is a valence or Rydberg or charge transfer excitation, what its symmetries are, whether it is a single or double excitation, etc. As a benchmark system for demonstrating the ability of NES-VMC to characterize different states, we study the carbon dimer (C\({}_{2}\)). Despite its small size, the carbon dimer has a complicated electronic structure with a large number of low-lying excited states [60; 47; 59]. Due to the existence of very strong bands in the visible spectrum, the carbon dimer is frequently detected in astrophysical measurements, and can be observed in comets rich in organic materials [61]. The exact bond order of C\({}_{2}\) is still a subject of some controversy - while molecular orbital theory would classify it as a double bond, valence bond calculations suggest it may be better described as a quadruple bond [62]. And the carbon dimer is one of the smallest molecules to have low-lying double excitations, a class of excited state which other methods often struggle with [48]. Correctly reconstructing the potential energy curves for different low-lying states requires correctly disentangling and characterizing these different states at different geometries.

We compute the 8 lowest-lying states of the carbon dimer at several different bond lengths using NES-VMC and the Psiformer Ansatz, and present the results in Figs. 3 and 4. At equilibrium (1.244 A), we classify the different states by computing their spin magnitude and their symmetries - both parity symmetry (u/g) where the electron positions \(\mathbf{r}\) are replaced by \(-\mathbf{r}\) and reflection symmetry (+/-) where the positions are flipped on the x-axis. We do not compute the orbital angular momentum operator, but confirm that we see the expected degeneracy, for instance \(\Pi\) states are doubly degenerate (one of each reflection symmetry). The oscillator strengths show several bright transitions, which we show in Fig. 3b. Due to the degeneracy of the \(\Pi\) states, we add the oscillator strengths together to give the total strength. We correctly identify the Phillips and Ballik-Ramay systems [63; 64], as well as the unnamed \(B^{1}\Delta_{g}\to A^{1}\Pi_{u}\) transition. We also find that the energy of the \(B^{1}\Delta_{g}\) energy closes matches the TBE in QUEST [48]. The \(A^{1}\Pi_{u}\), \(c^{3}\Sigma_{u}^{+}\) and \(b^{3}\Sigma_{g}^{-}\) states all have nearly the same energy, so correctly identifying the oscillator strengths for these transitions is very challenging.

To better understand the nature of each excited state, we compute the occupancy of the different natural orbitals. We first compute the one-electron reduced density matrix (1-RDM) for each single-state Ansatz in a large basis set and then diagonalize these matrices to find the natural orbitals, as described in more detail in Sec. C.2. In this case, using the Hartree-Fock orbitals as the basis set, we find that all 1-RDMs are nearly diagonal, that is the natural orbitals closely match the Hartree-Fock molecular orbitals. We see in Fig. 3d that all states above the ground state involve excitation of electrons into the \(2p_{z}\sigma_{g}\) orbital. The \(\Pi\) states are well-described by single excitations from one of the \(2p\pi_{u}\) orbitals while the \(c^{3}\Sigma_{u}^{+}\) state promotes an electron from the \(2s\sigma_{u}^{*}\) orbital. Finally, both the \(b^{3}\Sigma_{g}^{-}\) and \(B^{1}\Delta_{g}\) states are double excitations of the \(2p\pi_{u}\) electrons into the \(2s\sigma_{u}^{*}\) orbital, as expected. Not only is NES-VMC able to predict double excitation energies correctly, but by having an explicit functional form

Figure 4: Potential energy curves of the low-lying excited states of C\({}_{2}\) which can be uniquely identified from their symmetries. Complete numerical results are given in
for the wavefunction Ansatz, we can compute quantities such as the reduced density matrices which allow us to derive insight about the nature of electronic excitations.

Predicting experimental excitation energies requires computing the energy difference between different states in their respective lowest energy configurations, the so-called _adiabatic_ excitation energy. To compute this for C\({}_{2}\), we repeated the equilibrium calculations at a number of different bond lengths. Wherever possible, we matched the energy levels at different geometries to the appropriate states based on the same symmetries as in Fig (a)a, and for five states we were able to reconstruct enough of the potential energy curve to identify the minimum energy for each. The results are shown in Fig. 4, smoothed by cubic interpolation. Taking the difference between the minimum energies of each state gives an estimate of the adiabatic excitation energy, which we show in purple in Fig. (b)b, and in 3 out of 4 cases we matched the experimental energy[47] to within roughly 0.01 eV. We did not estimate the zero-point vibrational energies, but believe this may explain the discrepancy in the \(c^{3}\Sigma_{u}^{+}\) state. This shows that not only can NES-VMC match other theoretical calculations of vertical excitation energies, but can predict experimental results to high accuracy.

### Twisted Ethylene

The excited states of ethylene (C\({}_{2}\)H\({}_{4}\)) across its potential energy surface present a challenging benchmark problem for many methods. As the torsion of the carbon double bond is increased, an avoided crossing occurs when the torsion angle is \(90^{\circ}\). Even for ground state calculations, DFT and single-reference coupled cluster calculations predict an unphysical cusp at this location[67]. Starting from the \(90^{\circ}\) torsion and bending the hydrogen atoms on one side inward (so-called "pyramidalization"), ethylene undergoes a conical intersection where the ground state transitions from a \(\pi\) to \(\pi^{*}\) highest occupied orbital (the \(N\) and \(V\) states, with term symbols \({}^{1}A_{g}\) and \({}^{1}B_{1u}\)). Modeling this conical intersection requires fundamentally multireference methods, and while time-dependent density functional theory (TD-DFT) struggles with this system[68], multireference configuration interaction (MR-CI) methods describe it well[58].

We compute the excited states of ethylene as the torsion angle is varied from \(0^{\circ}\) to \(90^{\circ}\), followed by variation of the pyramidalization angle from \(0^{\circ}\) to \(120^{\circ}\), enough to include the conical intersection of the \(N\) and \(V\) states. We try to match the geometry from previous MR-CI studies[58] as closely as possible. Results are shown in Fig. 5. There are also several low-lying triplet states of ethene, the \({}^{3}B_{1u}\) and \({}^{3}B_{3u}\) states, and so we calculated \(K=3\) excited states for all geometries, which we found was enough to find two singlet states for all geometries except at equilibrium, where we used \(K=5\) and took the highest state, as the \({}^{1}B_{3u}\) state has lower energy exclusively at equilibrium. We tried both the FermiNet and Psiformer, and did not find any significant difference in the results, so we show the Psiformer results here (though FermiNet results are included in Table 6). For comparison, in addition to TD-DFT[57] and MR-CI, we also compare against the PauliNet penalty method[18]. For consistency, we show the PauliNet penalty method without variance matching, though the difference is not large. All results are normalized so that the ground state energy at

Figure 5: Excited states and conical intersection of ethylene (C\({}_{2}\)H\({}_{4}\)). Our results (blue) are compared against TD-DFT[57] (purple), MR-CI[58] (green) and a penalty method used with the PauliNet[18] (red). The best estimate of the location of the conical intersection of the V and N states for each method is given by the vertical line in Fig. (b)b. Our method is in close agreement with MR-CI up to a constant shift, and agrees with the location of the conical intersection better than the PauliNet penalty method. Note that the \(\phi=0\) geometry in Fig. (b)b differs slightly from the \(\tau=90\) geometry in Fig. (a)a, as in Barbati _et al.[58]_. Complete numerical results are given in Table 6.
the equilibrium geometry is 0.

Qualitatively, the results from NES-VMC closely match MR-CI. The spurious cusp when the torsion angle is 90\({}^{\circ}\) is avoided, and the error in the ground state relative to MR-CI is smaller than for the PauliNet penalty method across torsion angles. The non-parallelity error in the \(V\) state relative to MR-CI is lower for our method than the PauliNet penalty method, and our predicted location for the conical intersection (\(\sim\)97.5 degrees) is closer to the MR-CI value (\(\sim\)96 degrees) than the predicted PauliNet penalty method value (\(\sim\)100 degrees). There is a nearly constant shift in the energy of the \(V\) state on the order of several tenths of an eV relative to MR-CI, and a shift in the energy of the \(N\) state which grows as the pyramidal angle grows. Increasing the number of excited states and using a different Ansatz did not seem to make a difference. We note that when using the equilibrium geometry for ethylene from QUEST in Sec III.2 as opposed to the geometry from MR-CI, our results agreed with the theoretical best estimates to within chemical accuracy. The overall agreement with experimentally relevant quantities like the location of the conical intersection is in excellent agreement with other highly accurate theoretical studies, and so we are confident that NES-VMC is able to capture the important behavior of this system across the potential energy surface.

### Benzene

Finally, as a demonstration of the ability of our method to scale to larger systems, we applied NES-VMC with both the FermiNet and Psiformer to benzene. Benzene is a common benchmark for excited state methods for medium-sized molecules, so there is abundant data for us to compare against. For VMC, in addition to the penalty method of Entwistle _et al._[18], there is also the penalty method of Pathak _et al._[17], which is used in combination with a traditional Slater-Jastrow Ansatz, and uses a different form of the penalty function which allows for unbiased gradients. On top of VMC results and coupled-cluster-based TBEs from QUEST, we also compare against CASPT2[66] and TD-DFT with the PBE0 functional[65]. Results are shown in Fig. 6, with complex numerical results in Table 7. For our calculations, we used the same geometry as in QUEST[50].

As can be seen in Fig 5(a), NES-VMC with the Psiformer comes very close to reaching the TBE for all computed states. The FermiNet is not quite as accurate, and struggles with the highest energy \({}^{3}B_{2u}\) state. Inspection of the spin magnitude reveals that the highest excited state of the FermiNet converges to a mixture of a triplet and singlet state, which suggests that contamination from the \({}^{1}B_{1u}\) state is affecting the performance. The Psiformer is known to be much more accurate for ground state calculations on systems as large as benzene[38], so it is not surprising that the Psiformer is also better suited for computing the relative energy between states at this

Figure 6: Excited states of benzene. The NES-VMC results (green and blue) are compared against theoretical best estimates from QUEST[50; 54] alongside TD-DFT-PBE0[65], CASPT2[66], DMC with a Slater-Jastrow Ansatz and penalty method[17], and the PauliNet with a penalty method[18]. NES-VMC with the Psiformer Ansatz is competitive with state-of-the-art methods. All excitations are \(\pi\rightarrow\pi^{*}\) excitations, and the orbitals involved are visualized in Fig 5(b). Complete numerical results are given in Table 7.
scale. CASPT2 and TD-DFT methods are less accurate across the board, though this is not surprising as density functional methods are generally less accurate than wavefunction methods, and CASPT2 is generally intermediate in accuracy between DFT and coupled cluster. The penalty method of Entwistle _et al._ in combination with the PauliNet was only applied to the lowest excited state, and even on that, it only reaches CASPT2-level accuracy, even with variance matching (which we do not use in NES-VMC). The penalty method of Pathak _et al._, however, is much more accurate, generally reaching comparable levels of accuracy to NES-VMC with the Psiformer. We suspect this is due to the unbiased gradients in the method of Pathak _et al._. Additionally, the results reported in Pathak _et al._ include a diffusion Monte Carlo correction, which we omit, though this only reduces the root mean squared error by \(\sim\)0.1 eV. We note that NES-VMC with a sufficiently expressive Ansatz not only reaches levels of accuracy near coupled cluster, but does so without any free parameters to tune, unlike penalty methods.

To better understand the nature of the excitations computed, we inspected the density matrices of the respective states, similarly to the analysis of the carbon dimer in Sec III.3 and Fig. 3c. The natural orbitals are well described by the Hartree-Fock orbitals, and so the density matrices in the Hartree-Fock basis are all nearly diagonal. All five excited states for benzene we computed are single excitations from a \(\pi\) to \(\pi^{*}\) orbital, but interestingly, in the natural orbital picture, they are best described by exciting half an electron from two distinct \(\pi_{g}\) orbitals into two distinct \(\pi_{u}^{*}\) orbitals. These orbitals are visualized in Fig 6b. The ability to easily evaluate and analyze properties of wavefunctions other than just the energy is one of the advantages of explicitly computing the functional form of the wavefunction in VMC. Overall, our results on benzene show that NES-VMC can be usefully applied to larger systems and still produce accurate results, so long as the correct Ansatz is used.

## IV Discussion

We have presented a novel method for calculating excited state properties of quantum systems by variational quantum Monte Carlo (VMC), the Natural Excited States method (NES-VMC). NES-VMC has no free parameters to tune, and allows for unbiased estimation of energies and gradients, by reformulating a state-averaging approach as the problem of finding the ground state of an extended system. In much the same way that sampling from \(\psi^{2}\) is the natural way to compute ground state properties by VMC, we believe that NES-VMC is the natural variational principle for computing excited state properties. Additionally, it dovetails well with recent work on neural network Ansatze for many-body systems.

We have demonstrated the effectiveness of NES-VMC on a number of benchmark problems ranging from small atoms and molecules up to the benzene molecule. In all cases, NES-VMC is competitive with theoretical best estimates using coupled cluster methods for estimating energies and oscillator strengths, and can capture the behavior of double excitations and conical intersections. The optimized Ansatz can be used in downstream analyses to characterize the nature of the electronic structure of different excited states. We are confident that NES-VMC is as effective as any other method for computing excited states with QMC, with the added benefit of simplicity and generality.

Neural network Ansatze can be quite computationally expensive, which puts an upper limit on the scale of systems we considered. We believe that recent work on scaling and accelerating neural network Ansatze for many-electron problems [41] can be usefully applied to NES-VMC as well, which could allow these methods to be applied to problems for which no satisfactory solution exists today. While we focused on applications using neural network Ansatze, classic Ansatze like the Slater-Jastrow Ansatz can be scaled to much larger systems [25; 69]. Although our results suggest that more accurate Ansatze are quite important for achieving good performance, we look forward to finding out how well NES-VMC works in conjuction with these classic Ansatze on large problems.

Finally, while our experiments in this paper focused on molecular systems, that is, many-electron problems with open boundary conditions, NES-VMC is fully general and can be applied to _any_ quantum Hamiltonian. Excited state calculations with QMC are an important tool for studying nuclear physics [6], optical band gaps in condensed matter physics [13; 70], many properties of spin systems, as well as time dynamics and finite temperature phenomena. We are excited to see how NES-VMC can be applied to many of the most challenging open problems in many-body quantum mechanics in the future.

###### Acknowledgements.

 The authors would like to thank Matthew Foulkes, Denis Jacquemin, Michael Bearpark, Aron Cohen and Alex Gaunt for helpful discussions, and James Kirkpatrick, Annette Obika, Ali Eslami and Pushmeet Kohli for support.
