Accurate Computation of Quantum Excited States with Neural Networks

🎯 Pitch

This paper introduces a variational Monte Carlo method that recasts the search for multiple excited states as finding the ground state of an expanded Hamiltonian—requiring no free parameters, no explicit orthogonalization, and no penalty tuning—while neural network ansätze trained with it become the first deep learning approach to achieve accurate vertical excitation energies, including double excitations, on benzene-scale molecules.

1. Executive Summary

This paper proposes a variational Monte Carlo algorithm for computing the lowest excited states of quantum systems, called natural excited states for VMC (NES-VMC), which transforms the problem of finding multiple excited states into that of finding the ground state of an expanded system by constructing a determinant of single-state ansätze and promoting the local energy from a scalar to a matrix. The method is validated on atomic spectra (first-row atoms from lithium to neon), small-molecule vertical excitation energies and oscillator strengths from the QUEST database, and challenging systems including the carbon dimer, twisted ethylene, double excitations in molecules up to cyclopentadienone, and benzene, using the FermiNet and Psiformer neural network ansätze. NES-VMC requires no free parameters, no explicit orthogonalization of states, and yields unbiased estimates of both energies and off-diagonal observables such as transition dipole moments—achieving mean absolute errors within chemical accuracy (43 meV) on most systems and recovering oscillator strengths to within ~0.001 of theoretical best estimates, while correctly capturing phenomena like double excitations and conical intersections. The method is the first deep learning approach to achieve accurate vertical excitation energies on benzene-scale molecules, establishing that neural network ansätze can reach state-of-the-art accuracy for excited states when paired with a variational principle that avoids root-flipping, statistical bias, or penalty-weight tuning—but only when the number of computed states is sufficient to resolve the spectrum of interest.

2. Context and Motivation

The Core Problem: A Missing Variational Principle for Excited States

The fundamental challenge this paper tackles is deceptively simple to state but has resisted solution for nearly four decades: how do you variationally optimize multiple excited states of a quantum system simultaneously without them collapsing into one another, without introducing free parameters that must be tuned, and without requiring the ansatz to have a specific mathematical form?

This gap is both theoretical and practical. For ground states, variational quantum Monte Carlo (VMC) rests on a single, elegant variational principle: sample particle configurations from the square of the wavefunction, compute the local energy $E_L(x) = \psi^{-1}(x)\hat{H}\psi(x)$ , and minimize its expectation. This yields unbiased estimates of both the energy and its gradient, enabling optimization by stochastic gradient descent with convergence guarantees. The entire field of neural network wavefunction ansätze—from the FermiNet to the Psiformer—has been built on this foundation, achieving remarkable accuracy on ground state problems because the optimization objective is unambiguous.

For excited states, no analogous variational principle exists. The paper is explicit about this in its introduction:

"Despite almost four decades of work on QMC methods for excited states, 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 Ansätze and all excited states, regardless of symmetry."

This is the gap the paper aims to fill. The absence of such a principle means that every existing approach to excited-state VMC involves compromises—free parameters that must be heuristically tuned, statistical biases that accumulate in gradient estimates, restrictions on what functional forms the wavefunction can take, or prior knowledge about the states that may not be available.

Why This Problem Matters: From Photocatalysis to Nuclear Isomers

The paper motivates its work by pointing to the ubiquity of excited-state phenomena across chemistry and physics (Section I):

Photochemistry and photophysics: Understanding electronic excitations is critical for predicting fluorescence in quantum dots, light-induced molecular conformational changes, and photocatalytic activity. These are not niche applications—photocatalysis alone underpins emerging technologies for solar fuel generation and green chemistry.
Condensed matter and materials science: Excitations determine the optical band gap of semiconductors, which is the central quantity governing the behavior of solar cells, photosensors, LEDs, and lasers. Getting excitation energies wrong by even 0.1 eV can mean mischaracterizing whether a material is suitable for a given application.
Nuclear physics: Excited states are relevant to understanding metastable isomers, which have applications from nuclear clocks to targeted cancer therapy.
Astrophysics: As demonstrated later in the paper with the carbon dimer ( $C_2$ ), excited-state spectroscopy is how we detect and characterize molecules in comets and interstellar space. The visible bands of $C_2$ are frequently observed in astrophysical measurements of organic-rich comets.

Beyond these specific applications, there is a deeper scientific motivation: excited states are qualitatively different from ground states and much harder to compute accurately. The paper notes that inexpensive methods like time-dependent density functional theory (TD-DFT) "often give qualitatively incorrect results," while even gold-standard methods like multireference configuration interaction (MRCI) "can have large error bars on higher excited states." Machine learning approaches exist, but "most applications must be trained on expensive ab-initio calculations"—they don't improve the underlying calculations themselves. The paper positions its work in the less-explored space of using machine learning to "improve the ab-initio calculations themselves."

The Landscape of Prior Approaches and Their Limitations

The paper provides a thorough taxonomy of existing excited-state VMC methods (Section I), categorizing them into two broad families and identifying specific failure modes for each:

State-Targeting Methods: Finding One State at a Time

These methods aim to find a single excited state rather than the full spectrum. Variants include:

Energy-range targeting: Constrain the optimization to states within a specific energy window. Requires prior knowledge of where the state lies.
Symmetry targeting: Exploit known symmetries (spin, spatial symmetry) to isolate specific states. Fails when the excited state shares symmetries with lower states, or when the symmetry is unknown a priori.
Root-targeting: Optimize for the $k$ -th eigenstate directly. The paper notes these methods are "prone to root-flipping"—during optimization, the target state can swap with neighboring states, causing the optimization to chase the wrong root. This is not unique to QMC; it affects other computational paradigms as well.

The paper explicitly calls out the limitations of these approaches for deep neural network ansätze:

"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."

Multi-State Methods: Finding All States Simultaneously

These methods optimize multiple states together, either sequentially or jointly. The challenge is preventing different states from collapsing onto the same solution. Existing approaches include:

Penalty methods: Add an overlap penalty term to the variational objective to push states apart. The paper identifies several problems with existing penalty methods used with neural networks:

Biased gradients: Some penalty formulations introduce statistical bias into gradient estimates because they require solving generalized eigenvalue problems from stochastic matrix estimates. The paper explains: "This is usually avoided by accumulating matrix elements for long enough that this bias is small, this is not possible in the deep learning paradigm, where parameters are optimized by a large number of small, noisy steps."
Critical penalty threshold: Even with unbiased gradients, "convergence is only guaranteed if the strength of the penalty term is set above a certain critical value which is not known a priori." This is a free parameter that must be tuned heuristically.
Weight divergence: When optimizing multiple states simultaneously with unbiased penalty methods, "a weighting factor must be chosen for each state, otherwise the critical penalty threshold diverges." This introduces yet more free parameters.
Empirical performance: The paper notes that the penalty method of Entwistle et al. (2023), applied with the PauliNet ansatz, achieved errors of 1.4 eV on the $^3B_{1u}$ state of benzene—far from chemical accuracy—and even with variance matching corrections, the error was only reduced to 0.29 eV, still significantly worse than what NES-VMC achieves.

Explicit orthogonalization methods: Construct orthogonal ansätze by solving a generalized eigenvalue problem, sometimes with repeated re-orthogonalization during optimization. The critical limitation:

"Explicitly orthogonalizing Ansätze by solving a generalized eigenvalue equation for linear coefficients is usually only possible when the Ansatz is a linear combination of basis set functions, possibly including a shared Jastrow factor, which rules out neural networks."

This is the structural incompatibility that makes these methods unusable with deep neural network wavefunctions, which dominate modern ground-state VMC.

Linearized orthogonality approximations: Maintain only linearized approximations to wavefunction orthogonality, which the paper characterizes as insufficient.

The paper also notes a particularly revealing failure mode: the penalty method of Entwistle et al. (2023) with the PauliNet ansatz failed to find the $^2S$ and $^2P^\circ$ states of lithium—states that NES-VMC finds correctly. This is not a marginal accuracy issue; it is a qualitative failure to recover the correct physical states.

The Specific Challenge of Neural Network Ansätze

The paper's focus on neural network ansätze is not incidental—it reflects a fundamental tension in the field. Neural network wavefunctions like the FermiNet and Psiformer have achieved impressive accuracy on ground state calculations, "even exceeding coupled cluster (CCSD(T)) accuracy on some bond-stretching systems." The paper notes that VMC with these ansätze "scales as $O(N^3) - O(N^4)$ with system size, which is favorable for a wavefunction method."

However, neural network ansätze have "not yet matched the accuracy of other methods for excited states calculations as they scale to larger systems, despite initial forays." The paper poses the question directly:

"This could be due to the Ansatz not being sufficiently expressive, but it could also be due to the choice of variational principle."

This framing is crucial: the paper hypothesizes that the bottleneck is not the representational capacity of neural networks, but rather the optimization objective. Ground state VMC has a clear, universally accepted objective. Excited state VMC does not. The paper quotes the contrast:

"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, or optimization method used to achieve the variational bound, but the choice of objective to optimize is well-established. The same cannot be said for variational optimization of excited states."

How This Paper Positions Itself

The paper positions NES-VMC as the missing variational principle—one that has all the desirable properties that existing methods compromise on:

No free parameters: Unlike penalty methods, there is no penalty strength $\lambda$ to tune, no state weights $w_i$ to choose, no critical threshold to determine. The method has zero knobs.
No explicit orthogonalization: The determinant construction $\Psi = \det(\psi_i(x_j))$ naturally prevents collapse because if any two single-state ansätze become identical, the total ansatz $\Psi$ becomes zero everywhere. Orthogonality is not enforced; it emerges from the structure of the objective.
Unbiased gradient estimates: Because the method samples from $\Psi^2$ and computes expectations of the local energy matrix, gradient estimates are statistically unbiased. The paper explicitly connects this to convergence guarantees: "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."
Applicable to any ansatz: The method only requires that the single-state ansätze $\psi_i$ can be evaluated for a given particle configuration. There is no requirement that the ansatz be a linear combination of basis functions, no generalized eigenvalue problem to solve, no need for a Jastrow factor. This means it works "particularly well in conjunction with recent work on using neural networks as variational Ansätze."
Off-diagonal observables: The local energy matrix contains all the information needed to recover not just the energies of individual states, but also off-diagonal expectations like transition dipole moments—critical for computing oscillator strengths, which are experimentally measurable quantities.

The paper's title—"Natural Excited States"—is a deliberate reference to the natural choice of sampling distribution in ground state VMC. Just as $\psi^2$ is the natural sampling distribution for ground states (because it makes the Rayleigh quotient take a simple Monte Carlo form with unbiased gradient estimates), $\Psi^2 = \det(\psi_i(x_j))^2$ is argued to be the natural sampling distribution for excited states.

Reconciling Conflicting Prior Results

The paper's survey of prior work reveals an implicit explanation for why neural network ansätze had previously struggled with excited states while excelling at ground states: it was the variational principle, not the ansatz. The penalty method results on benzene (1.4 eV error without corrections) and lithium (missing the $^2S$ and $^2P^\circ$ states entirely) suggest that even when the ansatz has sufficient capacity, the optimization procedure can fail catastrophically. By introducing a variational principle with the same mathematical structure as ground state VMC—unbiased estimates, natural sampling distribution, no free parameters—NES-VMC aims to make excited-state optimization as reliable as ground-state optimization, unlocking the expressive power of neural network ansätze for excited-state problems.

This positioning is reinforced by the paper's choice to study two different neural network architectures (FermiNet and Psiformer) in parallel throughout the experimental sections. The consistency of results across architectures supports the claim that the variational principle, not the ansatz, is the critical factor. When both architectures succeed or fail in the same way on the same system (e.g., both missing the $^3A_1$ state of formaldehyde), it suggests a limitation of the optimization landscape rather than of either network's capacity.

3. Technical Approach

3.1 Reader Orientation

We are building a variational Monte Carlo (VMC) optimization procedure that takes a quantum Hamiltonian and a neural network ansatz and finds the $K$ lowest-energy eigenstates (wavefunctions) simultaneously. The system solves the problem of excited-state collapse—where different states merge into the same solution during optimization—not by adding penalty terms or enforcing explicit orthogonality, but by reformulating the multi-state problem as a ground-state problem on an expanded system where the total wavefunction is a Slater determinant of single-state wavefunctions, making state collapse mathematically impossible.

3.2 Big-Picture Architecture (Diagram in Words)

The system has four major components:

Single-State Ansätze ( $\psi_1, \dots, \psi_K$ ) — $K$ separate neural network wavefunctions (FermiNet or Psiformer), each taking $N$ particle coordinates and returning a scalar amplitude. These are the objects being optimized.
Total Ansatz ( $\Psi$ ) — a determinant constructed from the single-state ansätze evaluated at $K$ different particle configurations: $\Psi(x^1, \dots, x^K) = \det(\psi_i(x^j))$ . This is the joint wavefunction of an extended system with $K \times N$ particles.
Sampling, Local Energy Matrix, and Gradient Estimation — samples $K$ sets of particle positions from a distribution proportional to $\Psi^2$ , computes a $K \times K$ matrix of local energies $E_L(x) = \Psi^{-1}(x) \hat{H} \Psi(x)$ , and estimates unbiased gradients of the trace of this matrix with respect to the neural network parameters.
State Demixing and Observable Recovery — after optimization converges, diagonalizes the accumulated local energy matrix to recover individual state energies and eigenvectors, then uses these eigenvectors to transform other accumulated observable matrices to extract on-diagonal and off-diagonal expectations (including transition dipole moments).

Information flows as follows: at each optimization step, $K$ independent particle configurations are sampled from $\Psi^2$ via Markov chain Monte Carlo → the total ansatz $\Psi$ and its reciprocal are computed → the Hamiltonian is applied to each single-state ansatz at each configuration, forming the matrix $\hat{H}\Psi(x)$ → the local energy matrix $E_L(x) = \Psi^{-1}(x) \hat{H}\Psi(x)$ is computed → the trace of $E_L(x)$ provides an unbiased estimate of the sum of the $K$ lowest energies → gradients of this trace with respect to all network parameters are estimated without bias → parameters are updated via stochastic gradient descent → periodically, the accumulated local energy matrix is diagonalized to extract individual state energies and to generate transformation matrices for other observables.

3.3 Roadmap for the Deep Dive

First, the mathematical foundation: how NES-VMC transforms the problem of finding $K$ excited states into finding the ground state of an expanded Hamiltonian, and why this transformation is exact. This establishes why the method works at all.
Second, the sampling distribution and local energy matrix: why sampling from $\Psi^2$ (rather than from individual $\psi_i^2$ ) is the key innovation, and how the local energy generalizes from a scalar to a matrix. This is the core mechanism that enables unbiased gradient estimation.
Third, the gradient estimation: how unbiased gradients of the total energy are computed from samples, building on the standard VMC gradient estimator. This explains why the method has convergence guarantees.
Fourth, the state demixing: how diagonalizing the accumulated local energy matrix recovers individual state energies and the transformation needed to compute arbitrary observables, including off-diagonal ones.
Fifth, the neural network architecture and pretraining: how the FermiNet and Psiformer are extended to represent multiple states, and how pretraining ensures linear independence at initialization.
Sixth, practical considerations: numerical stability, pseudopotentials, singlet targeting, and the bottleneck envelope—the engineering details that make the method work at scale.

3.4 Detailed, Sentence-Based Technical Breakdown

This is primarily a methodological paper introducing a new variational principle for excited states. The core idea is that by forming a Slater determinant of single-state wavefunctions and sampling from its square, the problem of computing $K$ excited states becomes mathematically equivalent to computing the ground state of a $K$ -fold larger system, enabling unbiased energy and gradient estimation with the same mathematical structure as standard ground-state VMC.

The Mathematical Foundation: An Extended Ground-State Problem

The starting point is the standard quantum mechanical problem: given a Hamiltonian $\hat{H}$ acting on $N$ particles, find its $K$ lowest eigenfunctions $\psi_1^\star, \dots, \psi_K^\star$ with corresponding energies $E_1 \leq E_2 \leq \dots \leq E_K$ . The paper's insight is to construct an extended Hamiltonian $\hat{\mathcal{H}}$ that acts on $K$ copies of the $N$ -particle system:

$\hat{\mathcal{H}} = \hat{H}_1 \oplus \hat{H}_2 \oplus \dots \oplus \hat{H}_K$

where $\hat{H}_i$ is the original Hamiltonian $\hat{H}$ acting only on particle set $i$ , and acts as the identity on all other particle sets. In plain language: $\hat{\mathcal{H}}$ is the sum of $K$ independent copies of the original Hamiltonian, each acting on its own set of $N$ particles.

The ground state of this extended Hamiltonian has energy equal to $E_1 + E_2 + \dots + E_K$ —the sum of the $K$ lowest energies of the original system. Moreover, the ground state wavefunction of $\hat{\mathcal{H}}$ is a Slater determinant of the $K$ lowest eigenfunctions of $\hat{H}$ :

$\Psi^\star(x^1, \dots, x^K) = \det\begin{pmatrix} \psi_1^\star(x^1) & \dots & \psi_K^\star(x^1) \\ \vdots & & \vdots \\ \psi_1^\star(x^K) & \dots & \psi_K^\star(x^K) \end{pmatrix}$

where $x^i = (x^i_1, \dots, x^i_N)$ denotes the $N$ particle coordinates for particle set $i$ , and $\psi_j^\star(x^i)$ is the $j$ -th eigenfunction evaluated at the $i$ -th particle configuration.

Why this form: The determinant structure is crucial for two reasons. First, it is antisymmetric with respect to exchanging particle sets—formally, this is because we are treating each $N$ -particle set as a "fermion" in the extended system, and fermionic wavefunctions must be antisymmetric. Second, and more practically important: if any two single-state ansätze become identical ( $\psi_i = \psi_j$ ), two rows of the matrix become identical and the determinant vanishes identically ( $\Psi = 0$ ). This means the total ansatz cannot represent a situation where two states have collapsed together—the optimization is structurally prevented from finding such solutions because they would have infinite energy (the local energy would be undefined). This is the mechanism that replaces explicit orthogonality constraints or penalty terms: degeneracy is avoided by construction, not by enforcement.

The paper explicitly notes that any linear combination of the single-state ansätze $\psi_i' = \sum_j a_{ij} \psi_j$ changes $\Psi$ only by a constant multiplicative factor ( $\det(A)$ ), which does not affect the physics. This means the optimization does not need to find the "correct" individual states—any set spanning the same subspace is equally valid, and the demixing step (Section S2) will recover the physically meaningful individual states.

The Rayleigh Quotient for the Extended System

To find the ground state of $\hat{\mathcal{H}}$ , we minimize the Rayleigh quotient:

$\Psi^\star = \arg\min_\Psi \frac{\langle \Psi | \hat{\mathcal{H}} | \Psi \rangle}{\langle \Psi | \Psi \rangle}$

where $\langle \cdot | \cdot \rangle$ denotes integration over all $K \times N$ particle coordinates. The numerator expands as:

$\langle \Psi | \hat{\mathcal{H}} | \Psi \rangle = \sum_{i=1}^K \langle \Psi | \hat{H}_i | \Psi \rangle$

because $\hat{\mathcal{H}}$ is the sum of the individual Hamiltonians.

Using the matrix determinant lemma, the paper shows (Eq. 15) that this Rayleigh quotient can be rewritten in a compact matrix form:

$\frac{\langle \Psi | \hat{\mathcal{H}} | \Psi \rangle}{\langle \Psi | \Psi \rangle} = \text{Tr}\left[ S^{-1} \bar{H} \right]$

where $S$ is the $K \times K$ overlap matrix between single-state ansätze:

$S = \begin{pmatrix} \langle \psi_1^2 \rangle & \dots & \langle \psi_1 \psi_K \rangle \\ \vdots & & \vdots \\ \langle \psi_K \psi_1 \rangle & \dots & \langle \psi_K^2 \rangle \end{pmatrix}$

and $\bar{H}$ is the $K \times K$ matrix of Hamiltonian expectations:

$\bar{H} = \begin{pmatrix} \langle \psi_1 | \hat{H} | \psi_1 \rangle & \dots & \langle \psi_1 | \hat{H} | \psi_K \rangle \\ \vdots & & \vdots \\ \langle \psi_K | \hat{H} | \psi_1 \rangle & \dots & \langle \psi_K | \hat{H} | \psi_K \rangle \end{pmatrix}$

What this equation computes: The trace $\text{Tr}[S^{-1} \bar{H}]$ is the sum of the generalized eigenvalues of the matrix pair $(\bar{H}, S)$ , which at the minimum equals the sum of the $K$ lowest eigenvalues of $\hat{H}$ . If the single-state ansätze were orthonormal ( $S = I$ ), this would simply be the sum of their individual energies. When they are not orthonormal, the $S^{-1}$ factor corrects for the non-orthogonality, projecting the energy matrix into a basis where the states are effectively orthogonalized. This is the theoretical foundation for why the method works without explicit orthogonalization: the trace objective automatically accounts for whatever overlap exists between states.

Why not optimize this directly? Computing $S$ and $\bar{H}$ requires high-dimensional integrals over all $N$ -particle configurations for each pair of states. For $K$ states, this is $O(K^2)$ integrals, each requiring Monte Carlo sampling from some distribution. Worse, to get unbiased gradient estimates, one would need to sample from distributions proportional to $\psi_i \psi_j$ (off-diagonal) or $\psi_i^2$ (diagonal), which are different for each matrix element. This would require either $K^2$ separate MCMC chains or an importance-sampling scheme with potentially high variance. The paper's key innovation is finding a single sampling distribution— $\Psi^2$ —that gives unbiased estimates of the entire matrix $S^{-1} \bar{H}$ simultaneously.

The Natural Sampling Distribution: $\Psi^2$

The paper shows that by sampling all $K$ particle sets jointly from a distribution proportional to $\Psi^2$ , the matrix $S^{-1} \bar{H}$ can be estimated without bias from a single set of samples. Specifically, the paper proves (Eqs. 17-28) that:

$\text{Tr}\left[ S^{-1} \bar{H} \right] = \mathbb{E}_{x \sim \Psi^2} \left[ \text{Tr}\left[ \Psi^{-1}(x) \hat{H} \Psi(x) \right] \right]$

where $x = (x^1, \dots, x^K)$ denotes the complete configuration of all $K$ particle sets, and $\Psi(x)$ and $\hat{H}\Psi(x)$ are $K \times K$ matrices:

$\Psi(x) = \begin{pmatrix} \psi_1(x^1) & \dots & \psi_K(x^1) \\ \vdots & & \vdots \\ \psi_1(x^K) & \dots & \psi_K(x^K) \end{pmatrix}$

$\hat{H}\Psi(x) = \begin{pmatrix} \hat{H}\psi_1(x^1) & \dots & \hat{H}\psi_K(x^1) \\ \vdots & & \vdots \\ \hat{H}\psi_1(x^K) & \dots & \hat{H}\psi_K(x^K) \end{pmatrix}$

What these matrices contain: $\Psi(x)$ is simply the matrix of all single-state ansätze evaluated at all particle sets— $\psi_j(x^i)$ in row $i$ , column $j$ . $\hat{H}\Psi(x)$ is the same matrix but with the Hamiltonian applied to each entry: $\hat{H}\psi_j(x^i)$ is the result of the Hamiltonian operator acting on the $j$ -th wavefunction, evaluated at the $i$ -th particle configuration. Both matrices are $K \times K$ because there are $K$ states and $K$ particle sets.

The term inside the expectation is the local energy matrix:

$E_L(x) \equiv \Psi^{-1}(x) \hat{H} \Psi(x)$

What $E_L(x)$ computes: For each sampled configuration $x$ of all $K$ particle sets, $E_L(x)$ is a $K \times K$ matrix. Its trace $\text{Tr}[E_L(x)]$ is a single scalar that, when averaged over samples from $\Psi^2$ , gives an unbiased estimate of the sum of the $K$ lowest energies. The full matrix $E_L(x)$ contains richer information—when accumulated and diagonalized, its eigenvectors reveal how to demix the states to recover individual energies.

Why this is the natural generalization: In standard ground-state VMC with a single state $\psi$ , the sampling distribution is $\psi^2$ and the local energy is the scalar $E_L(x) = \psi^{-1}(x) \hat{H} \psi(x)$ . The paper's construction reduces to this exactly when $K = 1$ . For $K > 1$ , $\Psi^2$ replaces $\psi^2$ as the sampling distribution, and the local energy generalizes from a scalar to a matrix whose trace is the total energy. This is the analog of "natural" that the method's name references.

The paper provides a rigorous derivation (Eqs. 18-28 in Methods, Section S2) showing that the trace of the local energy matrix equals $\text{Tr}[S^{-1}\bar{H}]$ in expectation. The derivation uses the Leibniz formula for determinants, the fact that $\hat{\mathcal{H}} = \oplus_i \hat{H}_i$ , and the matrix determinant lemma to show that applying $\hat{\mathcal{H}}$ to $\Psi$ produces a sum of $K$ determinants, each being a rank-one update to the original $\Psi$ matrix. The matrix determinant lemma then expresses each term as $\Psi(x)$ times a correction factor, and summing over all $K$ particle sets yields the trace formula.

The connection to state collapse prevention is now clear: If two states are identical, $\Psi(x) = 0$ for all $x$ because the matrix has two identical columns and its determinant vanishes. The sampling distribution $\Psi^2$ is then zero everywhere, meaning no samples can be drawn and the local energy is undefined. The optimization cannot proceed—this is a "hard" barrier to collapse, not a soft penalty that can be overcome with a sufficiently large learning rate or insufficient penalty strength.

Computing the Local Energy Matrix in Practice

For each sampled configuration $x$ , we need:

The matrix $\Psi(x)$ : evaluate each single-state ansatz $\psi_j$ at each particle set $x^i$ . The paper notes that this is done in the log domain for numerical stability—compute $\log \psi_j(x^i)$ for all $i, j$ , subtract the maximum value, exponentiate, and assemble the matrix.
The matrix $\hat{H}\Psi(x)$ : for each $(i, j)$ , compute $\hat{H} \psi_j(x^i)$ . For the electronic Hamiltonian (Eq. 3), this involves the kinetic energy (Laplacian of $\psi_j$ with respect to all electron coordinates in $x^i$ ) and the potential energy (multiplying $\psi_j(x^i)$ by the sum of electron-electron, electron-nucleus, and nucleus-nucleus Coulomb terms). The kinetic energy is the expensive part—for neural network ansätze, it requires computing second derivatives of the network output with respect to all input coordinates.
The matrix inverse $\Psi^{-1}(x)$ : solve the $K \times K$ linear system. For the scale of $K$ studied in the paper (up to 10), this is a negligible computational cost compared to the neural network evaluations.
The local energy matrix: $E_L(x) = \Psi^{-1}(x) \hat{H}\Psi(x)$ . This is a $K \times K$ matrix.

Numerical stability considerations (Section S8): The paper describes several heuristics to handle points where entries of $\Psi(x)$ are zero (or numerically indistinguishable from zero). Since walkers are sampled from $\Psi^2$ , the total ansatz $\Psi$ cannot be zero at a sampled point, but individual entries $\psi_j(x^i)$ can be near zero. The paper's approach:

Always compute $\log \psi_j(x^i)$ and $\log(\hat{H}\psi_j(x^i))$ in the log domain, subtract the maximum log-value to avoid overflow/underflow, and only then exponentiate to get $\Psi(x)$ and $\hat{H}\Psi(x)$ .
If any entries of $\Psi$ are zero, run additional MCMC steps (up to 10 per iteration) to try to move away from the problematic configuration.
If the local energy is still undefined (contains NaNs), discard that walker from the computation of energies and gradients for that step.
If the gradient still has undefined values, skip that optimization iteration entirely.
"If more than 100 iterations are skipped in a row, the optimization fails, but we did not find that this happened in any of our experiments."

Tail behavior in large systems (Section S8): For systems with more than 25 electrons (butadiene, glyoxal, tetrazine, cyclopentadienone, benzene), the authors found the distribution of local energies became heavy-tailed, making the empirical mean unstable. To address this, they "took the median of elements of the local energy matrices over minibatches of 16 walkers, and took the empirical mean over these medians." The median-of-means estimator is more robust to outliers than the pure mean. This was used only for the largest systems; for smaller ones, the standard empirical mean was sufficient.

Unbiased Gradient Estimation

The total energy to be minimized is:

$\mathcal{L}(\theta) = \mathbb{E}_{x \sim \Psi^2} \left[ \text{Tr}[E_L(x)] \right]$

where $\theta$ represents all variational parameters of the single-state ansätze. The gradient of this expectation with respect to $\theta$ can be estimated without bias, following the same derivation as in standard ground-state VMC (Section S1):

$\nabla_\theta \mathcal{L} = 2 \mathbb{E}_{x \sim \Psi^2} \left[ \left( \text{Tr}[E_L(x)] - \mathbb{E}_{x' \sim \Psi^2} [\text{Tr}[E_L(x')]] \right) \nabla_\theta \log |\Psi(x)| \right]$

What this computes: The gradient is the covariance between the total local energy trace and the score function $\nabla_\theta \log |\Psi|$ . Configurations where the total energy is higher than average get positive weight in the gradient update (pushing the parameters away from producing such configurations), while those with below-average energy get negative weight (pushing parameters toward producing them more often).

Why this form: This is the standard VMC gradient estimator, valid because:

The expectation is over a distribution ( $\Psi^2$ ) that depends on the parameters $\theta$ , requiring the log-derivative trick.
The operator $\hat{H}$ is Hermitian (symmetric for real wavefunctions), so it is a valid observable for which the gradient expression holds.
The trace of the local energy matrix is a real scalar function of the configuration and parameters, so it enters the gradient expression exactly as a scalar local energy would in standard VMC.

The unbiasedness of the gradient estimate is crucial because it means that stochastic gradient descent with this estimator will converge to a local minimum under standard Robbins-Monro conditions (the paper cites Robbins and Monro, 1951). This is a formal guarantee that penalty methods lack—in penalty methods, the gradient of the penalty term may be biased if overlaps are estimated from finite samples, or the penalty strength may need to exceed an unknown critical threshold for convergence.

The gradient of $\log|\Psi|$ : Since $\Psi(x) = \det(\Psi(x))$ , we have:

$\nabla_\theta \log |\Psi(x)| = \nabla_\theta \log |\det(\Psi(x))| = \text{Tr}\left[ \Psi^{-1}(x) \nabla_\theta \Psi(x) \right]$

which requires the gradient of each single-state ansatz with respect to its parameters, evaluated at each particle set. This is computed by standard automatic differentiation through the neural network.

Batch estimation: Following exactly the same argument as in standard VMC (Eq. 10 in Section S1), the product of two expectations in the gradient can be estimated without bias from a finite batch by using the leave-one-out estimator:

$\hat{\nabla}_\theta \mathcal{L} = \frac{N - 1}{N} \sum_i \left( \text{Tr}[E_L(x_i)] - \frac{1}{N} \sum_j \text{Tr}[E_L(x_j)] \right) \nabla_\theta \log |\Psi(x_i)|$

where $N$ is the batch size (number of independent configurations of $K$ particle sets). The factor $(N-1)/N$ corrects for the fact that the same batch estimates both the mean energy and the deviation from the mean.

State Demixing and Observable Recovery

After optimization converges (or periodically during optimization for monitoring), the accumulated local energy matrix provides access to individual state properties. The paper shows that:

$\mathbb{E}_{x \sim \Psi^2} [E_L(x)] = S^{-1} \bar{H}$

What this is: This is the $K \times K$ matrix obtained by averaging the local energy matrix over samples. The equality states that this expectation equals the product of the inverse overlap matrix and the Hamiltonian matrix. This is the same matrix whose trace was minimized during optimization.

Now, if the single-state ansätze span the subspace of the lowest $K$ eigenstates, then the true eigenfunctions $\psi_i^\star$ are linear combinations of the ansätze: $\psi_i^\star = \sum_j A_{ij} \psi_j$ for some invertible matrix $A$ . In that case, the paper shows:

$\mathbb{E}_{x \sim \Psi^2} [E_L(x)] = A^{-1} \Lambda A$

where $\Lambda = \text{diag}(E_1, \dots, E_K)$ is the diagonal matrix of true eigenenergies.

What this means operationally: By diagonalizing the accumulated local energy matrix:

$\mathbb{E}[E_L] = U \Lambda U^{-1}$

we recover the individual state energies as the eigenvalues $\Lambda$ , and the matrix $U$ tells us how to rotate the single-state ansätze to recover the eigenstates. The eigenvectors $U$ satisfy $U = \Sigma A^{-1}$ where $\Sigma$ "is an arbitrary diagonal matrix which is not identifiable." This non-identifiability of the scale of each eigenvector means we cannot recover the individual wavefunctions uniquely, but we can recover all physical observables.

Recovery of on-diagonal observables: For any observable $\hat{O}$ (e.g., spin magnitude $\hat{S}^2$ , density matrix elements, dipole moments), accumulate the matrix:

$\mathbb{E}_{x \sim \Psi^2} \left[ \Psi^{-1}(x) \hat{O} \Psi(x) \right] = S^{-1} \bar{O}$

where $\bar{O}$ is the matrix of expectations $\langle \psi_i | \hat{O} | \psi_j \rangle$ . At convergence, this equals $A^{-1} \bar{O}^\star A$ where $\bar{O}^\star$ contains the true expectations between eigenstates. Transforming this matrix by the eigenvectors of the energy matrix:

$U^{-1} \mathbb{E}[\Psi^{-1} \hat{O} \Psi] U = \Sigma^{-1} \bar{O}^\star \Sigma$

The diagonal elements of this transformed matrix give $\langle \psi_i^\star | \hat{O} | \psi_i^\star \rangle$ exactly, because the $\Sigma$ factors cancel on the diagonal. This lets us compute any on-diagonal observable for each state individually, despite never having stored or optimized the individual states separately.

Recovery of off-diagonal observables (transition dipole moments): Off-diagonal elements of the transformed matrix equal $\frac{\sigma_i}{\sigma_j} \langle \psi_i^\star | \hat{O} | \psi_j^\star \rangle$ , where the $\sigma$ factors do not cancel. However, if we multiply the off-diagonal element by its transpose (i.e., square it element-wise):

$\Delta_{ij}^2 = \left( \frac{\sigma_i}{\sigma_j} \langle \psi_i^\star | \hat{O} | \psi_j^\star \rangle \right) \left( \frac{\sigma_j}{\sigma_i} \langle \psi_j^\star | \hat{O} | \psi_i^\star \rangle \right) = \langle \psi_i^\star | \hat{O} | \psi_j^\star \rangle^2$

the unknown $\sigma$ factors cancel. For the oscillator strength (Eq. 5 in the main text, Eq. 51 in Section S13C):

$f_{ij} = \frac{2}{3} \frac{m}{\hbar^2} (E_i - E_j) \sum_{k \in \{x,y,z\}} \langle \psi_i | q \hat{r}_k | \psi_j \rangle^2$

only the squared magnitude of the transition dipole moment is needed, so the sign ambiguity from the $\sigma$ factors is irrelevant. (In atomic units, $m = 1$ , $\hbar = 1$ , $q = -1$ for electrons.)

Variance estimation (Section S12): The paper describes a careful procedure for estimating the uncertainty in individual state energies. Since energies are not directly estimated but extracted via diagonalization, the variance must be propagated through the diagonalization. The paper computes the covariance matrix of the vectorized local energy matrix $\Sigma_{E_L} \in \mathbb{R}^{K^2 \times K^2}$ , then uses the eigenvectors $U$ of the mean energy matrix to transform this covariance:

$\Sigma_\Lambda = \left( U^{-1} \otimes U^T \right) \Sigma_{E_L} \left( U^{-1} \otimes U^T \right)^T$

where $\otimes$ is the Kronecker product and the transformation follows from the identity $\text{vec}(\Lambda) = (U^{-1} \otimes U^T) \text{vec}(\mathbb{E}[E_L])$ . The diagonal of $\Sigma_\Lambda$ , reshaped into a $K \times K$ matrix, gives the variance of each element of $\Lambda$ , and the diagonal of this reshaped matrix gives the variance of each individual state's energy. A similar procedure applies to other observables.

Neural Network Ansätze: FermiNet and Psiformer

The paper uses two neural network architectures to represent the single-state ansätze $\psi_i$ . Both are described in Section S3 and summarized in the main text. The extension from single-state to multi-state representation is straightforward: "Nothing is changed about the neural network architecture itself, just the number of orbitals is increased proportionally to the number of states."

Common structure: Both ansätze take electron positions $r_1, \dots, r_N$ and atomic nuclei positions $R_1, \dots, R_M$ as inputs, and output a scalar wavefunction amplitude. They assume spin-assigned electrons: the first $N_\alpha$ electrons are spin-up, the remaining $N_\beta = N - N_\alpha$ are spin-down. Both use generalized Slater determinants as their core:

$\psi(x) = \sum_k \omega_k \det \begin{pmatrix} \phi_1^k(x_1; \{x_{/1}\}) & \dots & \phi_N^k(x_1; \{x_{/1}\}) \\ \vdots & & \vdots \\ \phi_1^k(x_N; \{x_{/N}\}) & \dots & \phi_N^k(x_N; \{x_{/N}\}) \end{pmatrix}$

where $\phi_i^k(x_j; \{x_{/j}\})$ is a "generalized orbital"—a function of electron $j$ 's coordinates that also depends on all other electron coordinates $\{x_{/j}\}$ in a permutation-equivariant way. This generalizes the conventional Slater determinant (where each orbital depends only on $x_j$ ) by allowing each orbital to "see" the positions of all other electrons, which captures electron correlation effects that independent-particle models miss.

The FermiNet architecture (Pfau et al., 2020): Two parallel streams process the inputs:

One-electron stream: processes per-electron features through 4 layers, each with 256 units and hyperbolic tangent nonlinearities. At each layer $\ell$ , the activation for electron $i$ of spin $\alpha$ , denoted $h_i^{\ell\alpha}$ , is updated by concatenating its current activation with the mean activations of all same-spin electrons, all opposite-spin electrons, and the means of pairwise features with electrons of both spins (Eq. 32). This enforces permutation equivariance—the output for electron $i$ is invariant under permutations of other same-spin electrons.
Two-electron stream: processes pairwise electron features through parallel linear-tanh layers with 32 units per layer and residual connections.

The final one-electron stream activations are linearly projected to $N_\alpha$ spin-up and $N_\beta$ spin-down orbitals per determinant, multiplied by an exponentially decaying envelope (see below), and assembled into determinants. The paper uses 16 determinants ( $k = 1, \dots, 16$ ).

The Psiformer architecture (von Glehn et al., 2023): Replaces the FermiNet's custom equivariant layers with standard transformer components:

Each layer consists of a self-attention operation over electrons (4 attention heads, 64 units per head) with layer normalization and a residual connection, followed by parallel linear-nonlinear layers with tanh activation.
There is no two-electron stream; instead, a simple two-parameter Jastrow factor is multiplied after the determinant:

$J(r_1^\uparrow, \dots, r_N^\downarrow) = \sum_{i < j; \alpha = \beta} \frac{-1}{4} \frac{\alpha_{\text{par}}^2}{\alpha_{\text{par}} + |r_i^\alpha - r_j^\beta|} + \sum_{i,j; \alpha \neq \beta} \frac{-1}{2} \frac{\alpha_{\text{anti}}^2}{\alpha_{\text{anti}} + |r_i^\alpha - r_j^\beta|}$

This Jastrow factor captures electron-electron cusp behavior that the attention mechanism alone does not model accurately. The parameters $\alpha_{\text{par}}$ and $\alpha_{\text{anti}}$ control the cusp shape for same-spin and opposite-spin electron pairs, respectively.

Multi-state extension: For $K$ states, each state $\psi_i$ is given its own set of determinants:

$\psi_i(x) = \sum_{k=1}^{16} \omega_{ik} \det(\dots)$

where the orbitals are produced by the same neural network but with proportionally more outputs: instead of $N_\alpha + N_\beta$ orbitals per determinant, the network now outputs $K(N_\alpha + N_\beta)$ orbitals per determinant. The first $N_\alpha + N_\beta$ orbitals go to state 1, the next $N_\alpha + N_\beta$ to state 2, and so on. The network architecture itself is unchanged—only the final linear projection layer is wider.

The multiplicative envelope: To enforce the physical boundary condition that wavefunctions decay exponentially as electrons move far from nuclei, each orbital is multiplied by an envelope function:

$e_{j}^{k\alpha}(r_i^\alpha) = \sum_I \pi_{jI}^{k\alpha} \exp\left( \sigma_{jI}^{k\alpha} |r_i^\alpha - R_I| \right)$

where $I$ indexes nuclei, $\pi_{jI}^{k\alpha}$ and $\sigma_{jI}^{k\alpha}$ are learnable parameters. This is a sum of exponentially decaying functions centered on each nucleus, with the decay rate learned per-orbital, per-nucleus.

The bottleneck envelope (Section S4): For benzene (42 electrons, 12 nuclei), the default envelope becomes prohibitively expensive because the number of orbitals scales as $K \times N_{\text{det}} \times (N_\alpha + N_\beta) = 10 \times 16 \times (21 + 21) = 6720$ orbitals, each requiring a sum over 12 nuclei. The paper introduces a bottleneck:

$e_j^{k\alpha}(r_i^\alpha) = \sum_\ell w_{j\ell}^{k\alpha} \sum_I \pi_{\ell I}^\alpha \exp(\sigma_{\ell I}^\alpha |r_i^\alpha - R_I|)$

where $|\ell| = 32$ is the number of distinct envelope functions (much smaller than the 6720 orbitals), and each orbital takes a weighted sum of these basis envelopes. This reduces the expensive sum over nuclei from 6720 instances to 32 instances. The paper notes that this was used only for the largest systems and "did not have any numerical difficulties using the bottleneck envelope for benzene," though "on some systems... there were occasional numerical stability issues."

Hyperparameters (Section S10):

16 dense determinants
4 layers, 256 units per layer in the one-electron stream
FermiNet: 32 units per layer in the two-electron stream
Psiformer: 4 self-attention heads, 64 units per head, with layer normalization
Learning rate: 0.05 (following the Psiformer gradient scaling from von Glehn et al.)
Pretraining steps: 10,000 (default), increased to 50,000–100,000 for larger systems
Training iterations: 100,000–200,000 (optimal absolute energies require 200,000, but relative energies are "almost always well converged after 100,000")
MCMC: electrons sampled in blocks (not one at a time or all simultaneously), with each electron set treated as at least one block, subdivided into multiple blocks for larger systems (see Table S1)

Kinetic energy computation (Section S11): Computing the Laplacian for the kinetic energy is the computational bottleneck. Two methods are used:

Original FermiNet method: compute the gradient by backpropagation, then loop over computing columns of the Hessian via forward Jacobian-vector products. This computes the Laplacian as the trace of the Hessian.
Forward Laplacian (Folx, Gao and Günnemann, 2023): computes the Laplacian in a single forward pass, avoiding the loop over Hessian columns. "This leads to a 2-5x speedup depending on the system," with results "identical to within expected run-to-run variability."

Pretraining: Ensuring Linear Independence at Initialization

If the single-state ansätze are initialized to be linearly dependent (i.e., one state is a linear combination of others), then $\Psi = 0$ everywhere and the optimization cannot proceed. The paper addresses this with supervised pretraining of the neural network orbitals to match target Slater determinants (Section S7).

Ordered pretraining (default): Compute all single and double excitations of the Hartree-Fock ground state using PySCF with a cc-pVDZ basis set, order them by energy, and pretrain each single-state ansatz to match the Slater determinant for one of the $K$ lowest excitations. The pretraining uses Mean Absolute Error between the network orbitals and the target Hartree-Fock orbitals, with 10,000 steps by default. The paper notes that "the excitations found by Hartree-Fock do not necessarily correspond to the true lowest excitations, pretraining is primarily for finding a low-energy non-degenerate starting point, and the Ansatz would usually converge uniformly to the true lowest energy states."

Random pretraining: Leave core electron orbitals unchanged, but for valence orbitals, take random linear combinations of occupied valence orbitals plus the two lowest unoccupied Hartree-Fock orbitals. "Weights were sampled from a normal distribution with variance given by the inverse of the number of orbitals." This does not produce orthogonal initial states, but ensures linear independence (which is sufficient for $\Psi \neq 0$ ). The paper found this was necessary for systems with double excitations (HNO, nitrosomethane, butadiene, glyoxal) where ordered pretraining led to the optimization getting stuck at saddle points that missed the double excitations.

Saddle point behavior (Figure S2): The paper includes a diagnostic plot from the Psiformer on nitroxyl showing the optimization stuck at a plateau for thousands of iterations, with the gradient norm increasing transiently before a sudden drop in energy—"the expected behavior at a saddle point." After escaping, the network finds the correct double excitation state.

Sensitivity is not unique to NES-VMC: The paper explicitly compares NES-VMC with the ensemble penalty method (Section S6, Fig. S1) and finds that both methods "miss the exact same states" when using ordered pretraining. This supports the claim that "sensitivity to initialization is a general feature of how we optimize neural network Ansätze, and not an issue specific to NES-VMC."

Single-Targeting for Challenging Systems

For cyclopentadienone and tetrazine (Section S5), a large number of low-lying triplet states made it impractical to compute enough states to reach the singlet excitations of interest. The paper adds an extra term to the Hamiltonian:

$\hat{H}'_\lambda = \hat{H} + \lambda \hat{S}^2$

where $\hat{S}^2$ is the total spin magnitude operator (defined in Section S13A, Eq. 43):

$\hat{S}^2 = -\frac{N(N-4)}{4} + \sum_{i < j} \hat{P}_{ij}$

with $\hat{P}_{ij}$ swapping the spins of electrons $i$ and $j$ . Since $\hat{H}$ and $\hat{S}^2$ commute, the eigenfunctions are unchanged, but the eigenvalues are shifted: a state with total spin $S$ has its energy increased by $S(S+1)\lambda$ . For the singlet ground and excited states ( $S = 0$ ), the shift is zero. For triplets ( $S = 1$ ), the shift is $2\lambda$ . For quintets ( $S = 2$ ), the shift is $6\lambda$ .

By choosing $\lambda$ large enough, all triplet and higher-spin states are pushed above the singlet states of interest. The paper used $\lambda = 1$ for cyclopentadienone and $\lambda = 0.5$ for tetrazine, "roughly double the required scale to remove triplets in both cases." The lowest excited states of $\hat{H}'_\lambda$ are then exclusively singlets, which NES-VMC can find with a manageable $K$ .

Computational cost of $\hat{S}^2$ : Evaluating $\hat{S}^2$ requires summing over all opposite-spin electron pairs, swapping their positions, and evaluating the wavefunction ratio (Eq. 44). This adds $O(N_\alpha N_\beta)$ wavefunction evaluations per MCMC step—one for each opposite-spin pair. For large systems with many electrons, this can be significant, but it is only needed for these two challenging systems.

Pseudopotentials for Heavy Atoms

For molecules containing second-row atoms (HCl, H $_2$ S, H $_2$ CSi), the paper replaces core electrons with a correlation-consistent effective core potential (ccECP, Section S9). The rationale:

"We found that this generally led to more accurate results than all-electron calculations, possibly because the energy scale of the core electrons is far greater than the valence electrons for heavier atoms."

The pseudopotential is a nonlocal one-electron operator that replaces the Coulomb attraction and Pauli repulsion of core electrons with a fitted effective potential. The paper adapts the ccECP implementation from PySCF, following the approach for VMC of Fahy et al. (1990).

Scaling and Time Complexity

Asymptotic scaling (Section S11, Figure S3): The computational cost is dominated by:

Evaluating the determinant $\Psi = \det(\Psi)$ : $O(K^3)$ per forward evaluation.
Computing $\hat{H}\Psi(x)$ : requires evaluating $\hat{H}\psi_j(x^i)$ for all $K$ states at all $K$ particle sets— $O(K^2)$ wavefunction evaluations, each costing $O(N^3)$ to $O(N^4)$ for neural network evaluation (the dominant cost at practical $K$ ).
Matrix inversion and multiplication: $O(K^3)$ , negligible for $K \leq 10$ .

For the range of $K$ studied (up to 10), "the scaling is dominated by the cost of computing the matrix of local energies... Thus a quadratic scaling more accurately reflects the empirical behavior." With the original Laplacian method, "the asymptotic time complexity seems to be slightly below cubic... but this scaling does not dominate until $K > 10$ ." With Folx, the scaling is "roughly cubic" but all calculations are "roughly 4x faster."

Comparison to penalty methods: Penalty methods require explicit computation of the overlap matrix $S$ , which also scales as $O(K^2)$ . The paper notes that "when $K < 10$ , the roughly quadratic scaling of NES-VMC is comparable to penalty methods, which still require explicit computation of the overlap matrix, and do not yet work as well as NES-VMC with neural networks."

Network size ablations (Figure S4): The paper investigates how network capacity affects accuracy on H $_2$ O (a small system) by varying the number of layers (2–4), one-electron stream width (16–256), and either two-electron stream width (FermiNet, 8–32) or attention head units (Psiformer, 4–64). The results show that "the width of the one-electron stream is the primary factor in convergence to the ground truth" for the FermiNet, while "the Psiformer is more robust than the FermiNet to reduction in the width of the intermediate layers." The default architecture used throughout the paper (4 layers, 256 one-electron units) is "much larger than necessary for the smallest systems in this paper, and it remains to be seen how large we can push the system sizes considered before the Ansatz needs to be made larger."

Summary of Design Choices and Their Justifications

Sampling from $\Psi^2$ rather than individual $\psi_i^2$ : enables unbiased estimation of the full $S^{-1}\bar{H}$ matrix from a single sampling distribution, avoiding the need for $K^2$ separate MCMC chains or importance sampling with potentially high variance.
Determinant construction $\Psi = \det(\psi_i(x^j))$ : structurally prevents state collapse without penalty terms because identical states make $\Psi$ vanish. This replaces explicit orthogonality constraints with a hard mathematical barrier.
Trace objective $\text{Tr}[S^{-1}\bar{H}]$ rather than sequential optimization: the sum of generalized eigenvalues has the same minimum as the sum of individual state energies, even when states are not orthogonal. This means the optimization does not need to find orthogonal states—it only needs to find a set spanning the correct subspace.
Local energy matrix $E_L(x)$ rather than scalar local energy: generalizes the standard VMC estimator to multiple states while preserving unbiasedness of both energy and gradient estimates.
Diagonalization of accumulated $E_L$ for state demixing: avoids the need to assign states to ansätze during optimization, side-stepping root-flipping problems entirely. The demixing is done after optimization, when the subspace is well-converged.
Ensemble penalty method as comparison (Section S6): implemented with weights $w_i = i^{-1}$ and penalty strength $\lambda = 1$ , following Wheeler et al. (2024). The paper found it "became numerically unstable after a few hundred iterations" on benzene and achieved mean absolute errors an order of magnitude larger than NES-VMC on smaller systems where it did run (3.2 eV vs. 0.25 eV on HNO and HCF, excluding the missed state).

4. Key Insights and Innovations

Innovation 1: Reformulating Excited States as a Ground State of an Expanded System — Not Just a Trick, But a Unifying Variational Principle

The central conceptual move of this paper is not merely proposing another excited-state method, but recognizing that the entire problem can be recast as a ground-state problem on a larger system — and that this recasting yields a variational principle with the same mathematical structure as standard ground-state VMC. This is a fundamentally different category of solution from all prior work.

Before NES-VMC, the field treated excited-state VMC as a constrained optimization problem: minimize individual energies while preventing collapse through some external mechanism (penalties, orthogonalization, symmetry constraints, root-following). Every method added machinery — sometimes elegant, often heuristic — to fight against the natural tendency of energy minimization to find the ground state repeatedly. The implicit assumption was that excited states required additional structure beyond what ground-state VMC provides.

NES-VMC rejects this assumption. By constructing $\Psi = \det(\psi_i(x^j))$ and sampling from $\Psi^2$ , the method makes the existence of multiple distinct states a structural property of the ansatz rather than an externally enforced constraint. If two states are identical, $\Psi$ vanishes identically — there is nothing to optimize. The optimization cannot collapse states because the objective function becomes undefined at collapsed configurations. This is not a penalty that can be overcome with a large enough learning rate; it is a hard mathematical barrier.

What makes this more than a clever trick is the rigorous connection to an extended Hamiltonian: $\hat{\mathcal{H}} = \hat{H}_1 \oplus \cdots \oplus \hat{H}_K$ . The paper proves that minimizing the Rayleigh quotient of this extended system is exactly equivalent to finding the sum of the $K$ lowest eigenvalues of $\hat{H}$ . This means NES-VMC inherits all the theoretical properties of ground-state VMC — unbiased gradient estimates, convergence guarantees under Robbins-Monro conditions, and the natural sampling distribution $\Psi^2$ that makes the local energy take a simple form (now promoted from scalar to matrix). The paper is explicit that this is not an approximation or a heuristic: "Minimizing $\text{Tr}[E_\Psi]$ is mathematically equivalent to ground state VMC for an extended system that is $K$ times larger."

Comparison to prior work: Penalty methods (Entwistle et al., 2023; Pathak et al., 2021; Wheeler et al., 2024) add an overlap penalty to the objective, but this penalty is an external term that competes with energy minimization — the optimization must balance two objectives (lower energy vs. stay orthogonal), and this balance requires tuning the penalty strength above an unknown critical threshold. The paper's ensemble penalty experiments (Section S6, Fig. S1) demonstrate this fragility: on HNO and HCF, the ensemble penalty achieves mean absolute errors of 3.2 eV compared to 0.25 eV for NES-VMC, and becomes numerically unstable on benzene. Explicit orthogonalization methods (Filippi et al., 2009; Cuzzocrea et al., 2020) solve generalized eigenvalue problems during optimization, which introduces statistical bias in gradient estimates — a problem the paper identifies as particularly severe in the deep learning paradigm where "parameters are optimized by a large number of small, noisy steps." State-targeting methods (Zhao and Neuscamman, 2019; Choo et al., 2018) require prior knowledge of symmetries or energy ranges.

NES-VMC makes none of these compromises: no free parameters, no penalty terms, no generalized eigenvalue problems during optimization, no prior knowledge of symmetries. The extended-system framing is not one technique among many — it is a different category of solution that eliminates the constrained-optimization mindset entirely.

Significance beyond performance: The conceptual value of this reframing extends beyond the immediate numerical results. It provides a template for thinking about excited states as a natural extension of ground-state methods rather than a separate problem requiring separate machinery. The paper's title — "Natural Excited States" — is deliberate: just as $\psi^2$ is the natural sampling distribution for ground states because it makes the Rayleigh quotient take a simple Monte Carlo form, $\Psi^2$ is argued to be the natural sampling distribution for excited states because it makes the sum-of-energies objective take an equally simple form, with equally unbiased gradient estimates. This is a unification of ground and excited state VMC under a single conceptual framework.

Innovation 2: Identifying the Variational Principle — Not the Ansatz — as the Bottleneck for Neural Network Excited States

The paper makes a diagnostic contribution that shifts where the field should focus its attention. Before this work, it was unclear why neural network ansätze like the FermiNet, which had achieved remarkable accuracy on ground-state problems (exceeding CCSD(T) on some bond-stretching systems), had failed to match the accuracy of other methods for excited states as systems scaled up. The paper explicitly poses this as an open question in its introduction:

"This could be due to the Ansatz not being sufficiently expressive, but it could also be due to the choice of variational principle."

The experimental results provide a decisive answer: it is the variational principle, not the ansatz. The evidence comes from multiple angles:

Direct comparison on the same system: On benzene, the PauliNet with the penalty method of Entwistle et al. (2023) achieves an error of 1.4 eV on the $^3B_{1u}$ state, reduced to 0.29 eV with variance matching corrections. NES-VMC with the Psiformer achieves 0.088 eV on the same state — a threefold improvement over the corrected penalty result and more than an order of magnitude over the uncorrected one. If the ansatz were the bottleneck, switching from PauliNet to Psiformer would show the improvement. But the penalty method's poor performance with PauliNet already shows that a capable ansatz underperforms when paired with the wrong objective.
Consistency across architectures: Throughout the paper, the FermiNet and Psiformer produce similar results when both succeed, and fail similarly when both fail (e.g., missing the $^3A_1$ state of formaldehyde; struggling with the highest excited states of fluorine and neon). This consistency suggests that when failures occur, they reflect properties of the optimization landscape — saddle points, difficulty resolving closely-spaced states — rather than limitations of either architecture specifically.
Pretraining sensitivity affects both methods equally: The comparison in Fig. S1 shows that when ordered pretraining misses states (the $1A'$ double excitation of HNO, the singlet state of HCF at 5.65 eV), both NES-VMC and the ensemble penalty method miss the same states. The paper explicitly states: "Critically, the sensitivity to initialization is a general feature of how we optimize neural network Ansätze, and not an issue specific to NES-VMC." This is a carefully controlled negative result that isolates initialization as a shared challenge while showing that NES-VMC achieves dramatically better accuracy on the states it does find.
Network size ablations on H $_2$ O (Fig. S4): Varying the number of layers, one-electron stream width, and two-electron/attention head units shows that even heavily reduced networks (2 layers, 16 one-electron units) achieve near-chemical accuracy for the FermiNet, and the Psiformer is "more robust" still. This suggests the default architectures are substantially overparameterized for small systems, and the remaining errors are not due to insufficient expressivity.

This diagnostic insight has practical consequences for the field. It redirects research effort away from designing more expressive neural network architectures for excited states (since current architectures already have sufficient capacity) and toward either (a) adopting variational principles like NES-VMC that avoid the pathologies of penalty methods, or (b) addressing the shared challenge of initialization sensitivity that affects all methods.

Comparison to prior assumptions: The dominant narrative in the neural network VMC community had been that ansatz expressivity was the primary driver of accuracy — the FermiNet improved on Slater-Jastrow by introducing generalized orbitals; the Psiformer improved on the FermiNet by introducing attention mechanisms. The natural extrapolation was that excited states were harder because they required even more expressive ansätze. The paper's results suggest this extrapolation was wrong, or at least premature: the bottleneck was in the optimization objective, not the representational capacity. This is a reframing of the problem rather than an incremental accuracy improvement.

Innovation 3: The Local Energy Matrix as a Unified Object for Energy Estimation, State Demixing, and Observable Recovery

The paper introduces a specific mathematical object — the $K \times K$ local energy matrix $E_L(x) = \Psi^{-1}(x) \hat{H}\Psi(x)$ — and demonstrates that this single object, accumulated over the course of optimization, serves three distinct purposes that previously required separate machinery:

Energy estimation during optimization: The trace $\text{Tr}[E_L(x)]$ provides an unbiased estimate of the sum of the $K$ lowest energies, which serves as the optimization objective. This is the direct generalization of the scalar local energy in ground-state VMC.
State demixing after optimization: Diagonalizing the accumulated expectation $\mathbb{E}[E_L(x)] = U\Lambda U^{-1}$ recovers individual state energies (as the eigenvalues $\Lambda$ ) and the transformation $U$ needed to rotate the single-state ansätze into approximate eigenstates. This eliminates the need for root-tracking, sequential optimization, or solving separate generalized eigenvalue problems.
Observable recovery: The same eigenvectors $U$ transform any other accumulated observable matrix $\mathbb{E}[\Psi^{-1}\hat{O}\Psi]$ to recover both on-diagonal expectations $\langle \psi_i^\star|\hat{O}|\psi_i^\star\rangle$ and off-diagonal expectations squared $\langle \psi_i^\star|\hat{O}|\psi_j^\star\rangle^2$ , the latter being sufficient for oscillator strengths. This unifies energy and observable estimation into a single computational pipeline.

What makes this innovative is not the linear algebra itself — diagonalizing matrices to find eigenvalues is standard — but the recognition that all three tasks can be accomplished with the same matrix accumulated from the same samples, without any additional computation during optimization beyond what is already needed for the energy gradient. There is no separate "demixing phase" requiring additional sampling; the matrix $\mathbb{E}[E_L]$ is accumulated as a byproduct of optimization. There is no separate "observable estimation" requiring sampling from different distributions; the same $\Psi^2$ samples that estimate the energy also estimate any observable matrix $\mathbb{E}[\Psi^{-1}\hat{O}\Psi]$ .

Comparison to prior work: In penalty methods, computing observables for individual states typically requires identifying which ansatz corresponds to which physical state (since the states are assigned to specific ansätze during optimization), computing expectations for each state separately, and handling the fact that the states are only approximately orthogonal. The paper's approach decouples the ansatz index from the physical state identity — any linear combination of ansätze spanning the correct subspace is equally valid during optimization, and the physically meaningful states are recovered only at the end via diagonalization. This is a more robust separation of concerns: optimization finds the subspace; diagonalization finds the basis.

The off-diagonal observable recovery is particularly elegant. The non-identifiability of the scale factors $\sigma_i$ in $U = \Sigma A^{-1}$ means that off-diagonal expectations $\langle \psi_i^\star|\hat{O}|\psi_j^\star\rangle$ cannot be recovered individually (only up to the ratio $\sigma_i/\sigma_j$ ). But for oscillator strengths — one of the most experimentally relevant quantities — only the squared magnitude of the transition dipole moment is needed, and this squared quantity cancels the unknown scale factors exactly. The paper identifies this cancellation as a feature, not a bug: the method provides exactly the information needed for observable quantities, without requiring identification of the unobservable scale parameters.

Innovation 4: Verifier Over-Optimization Does Not Apply — The Determinant Provides a Structural Guarantee Against Collapse

The paper identifies and solves a problem that had been treated as requiring continuous tuning: preventing multiple states from collapsing into the same solution during optimization. Prior penalty-based approaches (Entwistle et al., 2023; Wheeler et al., 2024) treat this as an optimization problem — states collapse because the energy gradient pulls them toward the ground state, so an opposing force (the penalty) must be tuned to balance this pull. The penalty strength $\lambda$ becomes a hyperparameter that must exceed a critical threshold $\lambda > \max_{j>i}(E_j - E_i) \frac{w_i w_j}{w_i - w_j}$ , which depends on the unknown energy gaps being computed. This creates a circular dependency: you need to know the energy gaps to set the penalty, but you need the penalty to compute the energy gaps.

NES-VMC replaces this continuous tuning problem with a structural guarantee: if two single-state ansätze become identical, $\Psi = \det(\Psi) = 0$ everywhere, so $\Psi^2 = 0$ everywhere, so no samples can be drawn, so the local energy is undefined, so optimization cannot proceed. This is not a penalty that can be overcome — it is a hard mathematical barrier. The optimization landscape simply does not contain collapsed solutions as valid points.

This is analogous to — but stronger than — the Pauli exclusion principle in the paper's Slater determinant construction. In a conventional Slater determinant, two electrons cannot occupy the same orbital because the determinant vanishes. In NES-VMC, two "state-fermions" cannot occupy the same "state-orbital" because the determinant vanishes. But unlike the physical Pauli principle (which is a property of the exact wavefunction), the NES-VMC barrier is a property of the parameterization — any optimization trajectory that approaches state collapse will see $\Psi^2 \to 0$ , making the sampling distribution degenerate and the local energy diverge, pushing the optimization away.

Significance beyond performance: This insight reframes state collapse from an optimization challenge (requiring careful penalty tuning) to a parameterization challenge (requiring the ansatz to have the right algebraic structure). The determinant construction is one such structure, but the principle is more general: any ansatz that vanishes identically when states are linearly dependent will provide the same guarantee. This opens the door to alternative constructions beyond Slater determinants that might scale better or handle different symmetries.

Evidence of effectiveness: The paper never observes state collapse in any of its experiments, despite studying systems where closely-spaced states (fluorine, neon, the carbon dimer at equilibrium where $A^1\Pi_u$ , $c^3\Sigma_u^+$ , and $b^3\Sigma_g^-$ have nearly identical energies) would be expected to challenge any collapse-prevention mechanism. The spin magnitude analysis (Table S3, Fig. S5) confirms that states maintain their correct spin character throughout optimization. On ethylene (Table S8), the method correctly tracks the singlet-triplet crossing at the 90° torsion angle where the ground and first excited states nearly degenerate — a scenario where root-flipping and state collapse are most likely in other methods. The fact that NES-VMC navigates this crossing without any special handling (no symmetry targeting, no state following) demonstrates the robustness of the structural guarantee.

5. Experimental Analysis

Evaluation Methodology

Dataset. The paper validates the method across several distinct datasets rather than a single benchmark. For atomic spectra (Section III), the reference data comes from the NIST Handbook of Basic Atomic Spectroscopic Data, covering first-row atoms from lithium through neon with experimental energy levels averaged over fine-structure components weighted by degeneracy $m_J = 2J + 1$ . For molecular vertical excitation energies and oscillator strengths (Section IV), the paper uses the QUEST database (Veril et al., 2021; Loos et al., 2018, 2019, 2020, 2021, 2022), specifically the subset of 13 small molecules for which oscillator strengths were computed by Chrayteh et al. (2020) plus formaldimine, using theoretical best estimates (TBEs) from coupled-cluster calculations and basis-set extrapolation as the reference. For the carbon dimer (Section V), reference potential energy curves come from semistochastic heat-bath configuration interaction (SHCI; Holmes et al., 2017) and experimental adiabatic energies from Martin (1992). For twisted ethylene (Section VI), comparison is made against multireference configuration interaction (MR-CI; Barbatti et al., 2004) and TD-DFT (Mališ and Luber, 2020) at geometries from Barbatti et al. For double excitations (Section VII), five molecules with 24–42 electrons (nitrosomethane, butadiene, glyoxal, tetrazine, cyclopentadienone) are evaluated against the QUEST database and more recent diffusion Monte Carlo (Shepard et al., 2022) and CASPT3-corrected (Kossoski et al., 2024) results where the original QUEST TBEs were flagged as unsafe. For benzene (Section VIII), comparison is against coupled-cluster TBEs from QUEST (Loos et al., 2020), CASPT2 (Roos et al., 1992), TD-DFT-PBE0 (Adamo et al., 1999), and two prior VMC penalty method results (Pathak et al., 2021; Entwistle et al., 2023).
Base model(s). Two neural network architectures are used in parallel throughout all experiments: the FermiNet (Pfau et al., 2020) and the Psiformer (von Glehn et al., 2023). Both were originally developed for ground-state calculations and had demonstrated state-of-the-art accuracy exceeding coupled-cluster on some systems. The default hyperparameters for both architectures are: 16 dense determinants, 4 layers with 256 units per layer in the one-electron stream, and for the FermiNet 32 units per layer in the two-electron stream, for the Psiformer 4 self-attention heads with 64 units per head. These settings are kept fixed unless otherwise noted in specific system-dependent table entries (Table S1). The choice to study both architectures in parallel is deliberate: by showing consistency (or differences) between them, the paper can disentangle effects due to the variational principle (NES-VMC) from effects due to the specifics of any one neural network architecture.
Metrics. The primary metric throughout is the error in vertical excitation energy relative to a reference value. For atoms (Section III), this is the difference between NES-VMC's predicted excitation energy and the experimental value from NIST, reported in milliHartrees (mHa). For molecules (Sections IV–VIII), it is the difference between NES-VMC's predicted vertical excitation energy and the theoretical best estimate (TBE) from the QUEST database or other reference calculations, reported either in milliHartrees or electron-volts, with "chemical accuracy" defined as 1.6 mHa (~43 meV). For oscillator strengths (Section IV), the metric is the absolute difference between the computed oscillator strength and the TBE. Secondary metrics include the spin magnitude $\langle \hat{S}^2 \rangle$ (used to verify correct state identification), and in the carbon dimer analysis, the occupancy of natural orbitals and the energy drift along potential energy curves relative to SHCI. The paper also reports mean absolute error (MAE) aggregated across states for the double excitation systems.
Baselines. The paper compares against several classes of prior results. For VMC penalty methods with neural networks: Entwistle et al. (2023), who used the PauliNet ansatz with a penalty term on small molecules (BH, HCl, H $_2$ O, H $_2$ S, CO, C $_2$ H $_4$ , CH $_2$ O, CH $_2$ S, and HNO), reporting results both with and without a variance matching correction; the same work also reported results on the $^3B_{1u}$ state of benzene. For VMC penalty methods with conventional ansätze: Pathak et al. (2021), who used a Slater-Jastrow ansatz with an unbiased penalty method on benzene, including a diffusion Monte Carlo (DMC) correction. For conventional quantum chemistry methods: time-dependent density functional theory with the PBE0 functional (TD-DFT-PBE0; Adamo et al., 1999), complete active space second-order perturbation theory (CASPT2; Roos et al., 1992), and multireference configuration interaction (MR-CI; Barbatti et al., 2004). For projector QMC methods: diffusion Monte Carlo results from Shepard et al. (2022) on tetrazine. For reference standards: the QUEST database's theoretical best estimates from coupled-cluster calculations with basis set extrapolations, and experimental data for atoms (NIST) and the carbon dimer adiabatic energies (Martin, 1992). The paper also implements its own ensemble penalty method baseline (Section S6, Fig. S1, Table S4), following the formulation of Wheeler et al. (2024) with weights $w_i = i^{-1}$ , penalty strength $\lambda = 1$ , and identical neural network architectures to NES-VMC for a controlled comparison on HNO and HCF.
Generation budget / compute accounting. The paper does not use generation budgets or FLOPs-based accounting typical of large language model evaluations. Instead, computational cost is discussed in terms of scaling with the number of states $K$ and system size (Section S11, Figure S3). The key efficiency metrics are: empirical time-per-iteration as a function of $K$ on neon (Fig. S3), the quadratic-to-cubic scaling crossover around $K = 10$ , and the ~4x speedup from using the forward Laplacian method (Folx; Gao and Günnemann, 2023) versus the original column-by-column Hessian computation. Training length is reported as number of iterations (100,000–200,000) with pretraining steps (10,000–100,000) depending on system size (Table S1). For the largest systems, the paper uses a bottleneck envelope (Section S4) to reduce memory overhead by limiting the number of distinct envelope functions to 32 rather than $K \times N_\text{det} \times (N_\alpha + N_\beta)$ .
Cross-validation / statistical protocol. There is no cross-validation in the machine learning sense (training/validation/test splits). The paper uses standard Monte Carlo error estimation. For energies, the variance of individual state energies is estimated by propagating the covariance of the local energy matrix through the diagonalization, as described in Section S12: the covariance matrix $\Sigma_{E_L} \in \mathbb{R}^{K^2 \times K^2}$ of the vectorized local energy matrix is accumulated, then transformed by the eigenvectors $U$ of the mean energy matrix to obtain $\Sigma_\Lambda = (U^{-1} \otimes U^T) \Sigma_{E_L} (U^{-1} \otimes U^T)^T$ . The diagonal of $\Sigma_\Lambda$ , reshaped, gives the variance of each individual state energy. All energy values reported in the results tables are accompanied by their standard errors in the last significant digit (parentheses in the table entries). For large systems (>25 electrons), the paper uses a median-of-means estimator over minibatches of 16 walkers to handle heavy-tailed local energy distributions. A numerical stability heuristic skips optimization iterations if more than 100 consecutive iterations produce undefined gradients, though the paper reports this never occurred in practice.

Main Quantitative Results

Atomic Spectra (Section III, Figure 1, Table S2)

The paper computes the lowest excited states of all first-row atoms from lithium (3 electrons) to neon (10 electrons) to verify the correctness of the method against well-established experimental data, isolating the method's performance from complications like vibrational effects present in molecules.

Headline results. For all atoms from Li through O, 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." For Li, all ansätze correctly converge to the $2S$ and $2P^\circ$ states, which the paper notes were "missed by the use of the PauliNet Ansatz in combination with a penalty method" (Entwistle et al., 2023). For Be, the FermiNet with $K=10$ achieves errors of -0.06, -0.05, -0.10, -0.04, and -0.00 mHa on the $^3P^\circ$ , $^1P^\circ$ , $^3S$ , $^1S$ , and $^1D$ states respectively; the Psiformer achieves comparable accuracy (Table S2). For C, the FermiNet with $K=10$ achieves errors of 0.07, 0.11, 0.05, 0.11, and -0.40 mHa on the $^3P$ , $^1D$ , $^1S$ , and $^5S^\circ$ states.

Increasing $K$ improves accuracy for the highest states. For fluorine, "the error in the $4P$ states of fluorine is cut in half by increasing the number of states from 5 to 10" (FermiNet: from 10.8 to 9.68 mHa for the $^4P^\circ$ state at ~0.53 Ha excitation). For neon, increasing from $K=5$ to $K=10$ with the FermiNet reduces the error on the $2[3/2]^\circ$ state at ~0.61 Ha from 1.37 to 1.4–1.7 mHa (Table S2; the improvement is primarily that with $K=10$ the method reliably identifies the correct $2[1/2]^\circ$ state, which is near-degenerate with $2[3/2]^\circ$ and requires more states to disentangle).

Ansatz comparison. The Psiformer consistently reduces errors relative to the FermiNet for the most challenging cases. On fluorine, the Psiformer ( $K=10$ ) achieves 1.3 mHa on the $^4P$ state versus 2.7 mHa for the FermiNet ( $K=10$ ). On neon, the Psiformer achieves 1.56 mHa on the $2[1/2]^\circ$ state versus 1.7 mHa for the FermiNet, and 2.14 mHa versus 2.1 mHa on the $2[1/2]$ $3p$ state. The largest errors in the entire atomic series are for the highest $3p$ states of fluorine (9.68 mHa for FermiNet $K=10$ , Table S2) and neon (2.42–3.95 mHa for FermiNet $K=10$ ). The paper hypothesizes that for these cases, "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 Ansätze, the effects of individual states can be disentangled."

Rare failure cases. On boron with the Psiformer ( $K=10$ ), "the highest state seems to converge to the incorrect state... which seems to converge to the $2P^\circ$ state rather than the last $2D$ state" (Table S2: the error for the $^2P^\circ$ state is 0.02–0.03 mHa, while the $^2D$ state shows 0.07 mHa error with FermiNet $K=10$ but is replaced by what appears to be the $^2P^\circ$ state in the Psiformer results). This is a qualitative failure (wrong state assignment) rather than a quantitative error.

Oscillator Strengths on Small Molecules (Section IV, Figure 2, Table S3)

The paper evaluates NES-VMC on 13 small molecules from the QUEST database subset investigated by Chrayteh et al. (2020) for which both vertical excitation energies and oscillator strengths have benchmark theoretical values.

Headline results on vertical excitation energies. "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 43 meV) of the TBE." Specific examples (Table S3): the FermiNet achieves energy errors of -0.38 mHa (BH $^1\Pi$ ), -0.76 mHa (HCl $^1\Pi$ ), 0.00 mHa (H $_2$ O $^1B_1$ ), -0.31 mHa (H $_2$ O $^1A_2$ ), -1.06 mHa (H $_2$ S $^1B_1$ ), -0.2 mHa (BF $^1\Pi$ ), 0.4 mHa (CO $^1\Pi$ ), 0.35 mHa (C $_2$ H $_4$ $^1B_{3u}$ ), -0.1 mHa (CH $_2$ O $^1B_2$ ), and -1.02 mHa (CH $_2$ S $^1A_2$ ). The Psiformer shows comparable or slightly better accuracy, e.g., -0.44 vs. -0.38 mHa for BH $^1\Pi$ , -0.11 vs. 0.00 mHa for H $_2$ O $^1B_1$ , -0.32 vs. 0.35 mHa for C $_2$ H $_4$ $^1B_{3u}$ .

Headline results on oscillator strengths. "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." Specific examples (Table S3, $\Delta f$ columns): FermiNet achieves $\Delta f$ of -0.0004 (BH $^1\Pi$ ), -0.00072 (HCl $^1\Pi$ ), -3.3×10 $^{-5}$ (H $_2$ O $^1B_1$ ), -0.00097 (H $_2$ S $^1B_1$ ), -0.00295 (BF $^1\Pi$ ), 0.00123 (CO $^1\Pi$ ), 0.00368 (C $_2$ H $_4$ $^1B_{3u}$ ), and 0.00721 (CH $_2$ O $^1B_2$ ). The Psiformer achieves similar results. This is described as "much more accurate than the oscillator strengths presented in Entwistle et al., even when corrections using variance matching are applied."

Comparison to earlier VMC penalty method. Figure 2 includes black triangles for the PauliNet penalty method results from Entwistle et al. (2023), showing both raw results (pointing right) and variance-matched results (pointing left). The NES-VMC oscillator strengths consistently lie much closer to the TBE lines. The paper notes that "we do not use variance matching for any of the NES-VMC calculations," meaning the improvement comes from the method itself rather than post-hoc corrections.

Systems with anomalous behavior. Several exceptions are documented. On HCl: "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" (Table S3: Psiformer oscillator strength is -9.78×10 $^{-5}$ , essentially zero, versus the FermiNet's 0.0223 and TBE 0.0231). On formaldehyde (CH $_2$ O): "both the FermiNet and Psiformer fail to find the $^3A_1$ state at all, and the oscillator strength for the $^1B_2$ state diverges from the TBE by a significant margin" (FermiNet $\Delta f = 0.00721$ , Psiformer $\Delta f = 0.00275$ , compared to TBE; the Psiformer approximately halves the error). On H $_2$ S: "vertical excitation energies for systems with heavier atoms... are not quite as accurate as other results" (FermiNet error = -1.06 mHa on $^1B_1$ ). On thioformaldehyde (CH $_2$ S): the FermiNet has a 10.19 mHa error on the $^3B_2$ state—the largest error in the entire molecule set—while the Psiformer improves this to 8.12 mHa (Table S3). On HNO and HCF, "the network initialization must be carefully chosen to ensure convergence to the correct states," motivating the pretraining sensitivity analysis in Section S7 and Fig. S1.

Pretraining sensitivity analysis (Fig. S1, Tables S3–S4). The paper reveals that for HNO and HCF, the default "ordered" pretraining (matching Hartree-Fock excitations ordered by energy) causes the method to miss specific states. For HNO with ordered pretraining (Table S3), the FermiNet misses the $1A'$ state, while the Psiformer finds a state at 10.45 eV with $\langle \hat{S}^2 \rangle = 0.010$ which the paper identifies as likely the wrong state. With random pretraining, the FermiNet finds the $1A'$ state at 4.32 eV (0.2 mHa error) and the Psiformer finds it at 3.58 eV (-0.13 mHa error). For HCF with ordered pretraining, the FermiNet misses the singlet state at ~5.65 eV entirely and finds only a triplet at 6.77 eV. With random pretraining, the FermiNet finds the singlet at 5.75 eV (0.00327 oscillator strength), and the Psiformer finds it at 5.64 eV (0.00381 oscillator strength). The paper notes critically: "Critically, the sensitivity to initialization is a general feature of how we optimize neural network Ansätze, and not an issue specific to NES-VMC." The ensemble penalty method "misses the exact same states missed by NES-VMC" (Fig. S1, Table S4), and "is less numerically accurate than NES-VMC most of the time, reaching chemical accuracy on only a single identifiable state" with a mean absolute error of 3.2 eV versus 0.25 eV for NES-VMC on the identifiable states.

Pseudopotential improvement for heavy atoms. For HCl, H $_2$ S, CH $_2$ S, and H $_2$ CSi, the paper used pseudopotentials for second-row atoms. It states: "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." The degree of improvement is not quantified numerically, but the method achieved chemical accuracy on the excitation energies of HCl and H $_2$ S despite these being noted as challenging.

Omitted diatomic system. N $_2$ is explicitly omitted because "the lowest-lying excited states are all triplets." No further explanation is given, but it implies that without singlet states to verify against the benchmark oscillator strengths, the comparison would be less informative.

Carbon Dimer: Potential Energy Curves, State Classification, and Double Excitations (Section V, Figure 3, Tables S5–S7)

The carbon dimer (C $_2$ ) is studied as a small but electronically complex system with many low-lying excited states, visible astrophysical bands, low-lying double excitations, and a controversial bond order. The paper computes the 8 lowest-lying states at multiple bond lengths using the Psiformer.

Potential energy curves and comparison to SHCI (Figure 3a, Table S6). The paper computes potential energy curves at 10 bond lengths from 0.8 to 1.5 times the equilibrium length ( $a_{\text{eq}} = 1.244$ Å for the ground state). The NES-VMC curves are shifted to match SHCI at equilibrium (the paper states "shifted by 0.115 Ha to account for basis set effects," referring to the difference between the finite basis used by SHCI and the real-space representation of NES-VMC). The curves qualitatively match SHCI across all geometries, correctly capturing the different bonding characters of the various states: the $X^1\Sigma_g^+$ , $a^3\Pi_u$ , $A^1\Pi_u$ , $c^3\Sigma_u^+$ , and $b^3\Sigma_g^-$ states all have well-defined minima, while the $d^3\Pi_g$ state is dissociative. The paper reports: "The vertical energies match SHCI to within chemical accuracy at all bond lengths, and the absolute energy of a state drifts by about 5 mHa/Å relative to SHCI as the bond length changes."

Equilibrium properties (Tables S5–S7). At equilibrium ( $a_{\text{eq}}$ , Table S5), the energy ordering is: $X^1\Sigma_g^+$ (-75.92258 Ha), $a^3\Pi_u$ (-75.91335 Ha average of the degenerate pair), $b^3\Sigma_g^-$ (-75.87527 Ha), $c^3\Sigma_u^+$ (-75.87881 Ha), $A^1\Pi_u$ (-75.87743 Ha average), $B^1\Delta_g$ (-75.84563 Ha), and $d^3\Pi_g$ (not listed at equilibrium, appears only at compressed geometries). The paper notes that " $A^1\Pi_u$ , $c^3\Sigma_u^+$ , and $b^3\Sigma_g^-$ states all have nearly the same energy at equilibrium," making correct identification challenging. The oscillator strengths (Table S5) identify the Phillips ( $A^1\Pi_u \to X^1\Sigma_g^+$ , total $f = 0.00501$ summed over the degenerate pair), Ballik-Ramsay ( $b^3\Sigma_g^- \to a^3\Pi_u$ , $f = 0.00608$ ), and "unnamed $B^1\Delta_g \to A^1\Pi_u$ " ( $f = 0.00250$ ) systems. The oscillator strength of the Swan bands ( $d^3\Pi_g \to a^3\Pi_u$ ) is not computed since the $d^3\Pi_g$ state was not included in the equilibrium calculation (it is dissociative).

Adiabatic energies (Figure 3b, Table S7). By fitting cubic interpolations to the potential energy curves and finding their minima, the paper extracts adiabatic excitation energies and equilibrium bond lengths for each state. Compared to experiment (Martin, 1992), NES-VMC achieves: $a^3\Pi_u$ : 0.0837 eV (expt. 0.0891 eV, error = -5.4 meV), $b^3\Sigma_g^-$ : 0.8041 eV (expt. 0.7978 eV, error = +6.3 meV), $A^1\Pi_u$ : 1.0407 eV (expt. 1.0404 eV, error = +0.3 meV), $c^3\Sigma_u^+$ : 1.1299 eV (expt. 1.1312 eV, error = -1.3 meV), $B^1\Delta_g$ : 1.5031 eV (expt. 1.4980 eV, error = +5.1 meV). The paper states: "while SHCI had a maximum error of 0.03 eV and mean absolute error (MAE) of 0.02 eV on the first five excitations, our maximum error was 7 meV and MAE was 4 meV, a roughly fourfold improvement." The paper attributes the improvement to the "5 mHa/Å drift" in the absolute energy relative to SHCI that correlates with the equilibrium bond length, concluding that "the NES-VMC curve is closer to the ground truth than SHCI."

State classification by symmetry (Figure 3c). The paper classifies the states by evaluating each single-state ansatz at positions $\mathbf{r}$ and $-\mathbf{r}$ (parity: u/g) and by flipping $\mathbf{r}$ across the x-axis (reflection: +/-). The figure demonstrates clean separation into $g$ and $u$ states and correct identification of $\Sigma^+$ vs. $\Sigma^-$ reflection symmetry.

Natural orbital analysis (Figures 3d–3e). The paper computes the one-electron reduced density matrix (1-RDM) for each state by projecting onto a def2-TZVPD basis set (Section S13B) and diagonalizing to obtain natural orbitals and their occupancies. Figure 3d visualizes the 8 lowest natural orbitals, which "closely match the Hartree-Fock molecular orbitals." Figure 3e plots the occupancy of these orbitals for each state. The key finding: "all states above the ground state involve excitation of electrons into the $2p_z\sigma_g$ orbital." The $\Pi$ states ( $a^3\Pi_u$ , $A^1\Pi_u$ ) 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." Crucially, "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." This demonstrates that "not only is NES-VMC able to predict double excitation energies correctly, but by having an explicit functional form for the wavefunction Ansatz, we can compute quantities which allow us to derive insight about the nature of excitations."

Twisted Ethylene: Conical Intersection and Avoided Crossing (Section VI, Figure 4, Table S8)

Ethylene (C $_2$ H $_4$ ) under torsion and pyramidalization is a classic benchmark for multireference methods because it features an avoided crossing at the 90° twisted geometry and a conical intersection under further pyramidalization, where the $N$ ( $^1A_g$ ) and $V$ ( $^1B_{1u}$ ) states cross.

Torsion potential energy curve (Figure 4a). NES-VMC with the Psiformer tracks the two lowest singlet states as the torsion angle $\tau$ varies from 0° to 90°. Results are normalized so the ground state energy at equilibrium ( $\tau = 0$ ) is zero. The paper states: "Qualitatively, the results from NES-VMC closely match MR-CI. The spurious cusp when the torsion angle is 90° 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 NES-VMC ground state energy lies below the MR-CI curve by an approximately constant shift, while the $V$ state shows a shift "on the order of several tenths of an eV relative to MR-CI." The PauliNet penalty method results (plotted in red) show larger deviations from MR-CI.

Pyramidalization and conical intersection (Figure 4b). Starting from $\tau = 90°$ (twisted geometry), the pyramidalization angle $\phi$ is varied from 0° to 120° to drive the system through the conical intersection where the $N$ and $V$ states cross. The paper reports: "our predicted location for the conical intersection (~97.5 degrees) is closer to the MR-CI value (~96 degrees) than the predicted PauliNet penalty method value (~100 degrees)." The $V$ state energy shows "a nearly constant shift... on the order of several tenths of an eV relative to MR-CI," while the $N$ state energy shift "grows as the pyramidalization angle grows." Table S8 provides the raw energies and spin magnitudes: at $\phi = 97.5°$ , the Psiformer reports energies of -78.39400 Ha ( $N$ state, $\langle \hat{S}^2 \rangle = 2.000$ ) and -78.39106 Ha ( $V$ state, $\langle \hat{S}^2 \rangle = 0.001$ ), with a gap of only 2.94 mHa between them — indicating that the method correctly resolves the near-degeneracy without state collapse.

Impact of $K$ and number of states. The paper notes that at equilibrium ( $\tau = 0$ ), "we used $K = 5$ and took the highest state, as the $^1B_{3u}$ state has lower energy exclusively at equilibrium." Table S8 confirms this: at $\tau = 0$ with $K = 5$ , the states are identified as $^1B_{3u}$ (0.001), $^3B_{1u}$ (2.000), $^3B_{3u}$ (1.996/2.000), $^1B_{3u}$ (0.007/0.001), and $^1B_{1u}$ (0.005/0.002). At $\tau = 0$ with $K = 3$ , the third state shows $\langle \hat{S}^2 \rangle = 1.413$ (FermiNet) and 0.001 (Psiformer), confirming that with only $K = 3$ states, the singlet and triplet $B_{3u}$ states mix, requiring $K \geq 5$ at equilibrium to cleanly separate them. For all other torsion angles, $K = 3$ suffices.

Comparison to QUEST benchmark at equilibrium. The paper notes a discrepancy: "when using the equilibrium geometry for ethylene from QUEST in Sec IV as opposed to the geometry from MR-CI, our results agreed with the theoretical best estimates to within chemical accuracy." This implies that the shift relative to MR-CI at equilibrium is partly due to geometry differences between the two reference calculations.

Double Excitations in Larger Molecules (Section VII, Figure 5, Table S9)

Double excitations—where two electrons are simultaneously promoted—are notoriously difficult for single-reference methods. The paper evaluates NES-VMC on five molecules with 24–42 electrons known to have low-lying states with significant double excitation character: nitrosomethane, butadiene, glyoxal, tetrazine, and cyclopentadienone.

Headline results. "On all double excitation systems with the Psiformer, and 4 out of 5 systems with the FermiNet, NES-VMC is in excellent agreement with the best computational results." The FermiNet and Psiformer achieve mean absolute errors relative to the TBE of: 15 and 21 meV (nitrosomethane), 84 and 38 meV (butadiene), 167 and 28 meV (glyoxal), 45 and 54 meV (tetrazine), and 92 and 66 meV (cyclopentadienone), respectively (Table S9). For the Psiformer, this is within chemical accuracy (43 meV) for all systems except tetrazine (54 meV) and cyclopentadienone (66 meV).

Butadiene — resolving a historical controversy. Butadiene is highlighted as "the smallest conjugated organic molecule," with its two lowest singlet excited states ( $1^1B_u$ and $2^1A_g$ ) historically misordered in the literature. The paper notes that "the exact ordering... have been the subject of controversy for many years, only being resolved in the last decade or so." NES-VMC correctly predicts the ordering with both ansätze. The FermiNet and Psiformer match the TBE for $1^1B_u$ to within 92 and 60 meV respectively, and for $2^1A_g$ to within 62 and 9 meV — "a remarkably high degree of agreement for such a notorious system." The $2^1A_g$ state "is known to have roughly 30% double excitation character, making it especially challenging to compute."

Glyoxal — FermiNet struggles. Glyoxal (227.8 Ha, 30 electrons) is the one system where the FermiNet performs significantly worse than the Psiformer: MAE of 167 meV vs. 28 meV. The paper notes: "this is at the scale where the FermiNet is known to perform worse than the Psiformer at ground state calculations, so it is not surprising that it struggles on some systems of this size." This is consistent with the known scaling limitations of the FermiNet architecture.

Tetrazine — correcting the QUEST database. Tetrazine is one of the two systems (along with cyclopentadienone) where "the original TBEs from QUEST were known to be unsafe for double excitations, and only very recently did more accurate calculations from QMC and CASPT3 resolve discrepancies as large as almost 1 eV for some states." The paper reports: "On tetrazine, the Psiformer is within 0.1 eV of the best estimates of the $2^1A_g$ vertical excitation energy, while the previous TBE in QUEST was off by nearly 1 eV." Specifically, the QUEST TBE (unsafe) for $2^1A_g$ was 169 mHa, while NES-VMC with the Psiformer finds 185.9 mHa (Table S9) and the DMC correction (Shepard et al., 2022) finds 183 mHa — a correction of ~16 mHa (~0.44 eV) from the unsafe TBE to the accurate value.

Cyclopentadienone — resolving discrepancies between QMC and CASPT3. For cyclopentadienone, the two recent corrections to the QUEST database (DMC: 217 mHa for $2^1A_1$ , 253 mHa for $3^1A_1$ ; CASPT3: 213 mHa, 247 mHa) disagree by 4–6 mHa (~0.1–0.15 eV) themselves. NES-VMC with the Psiformer achieves 217.1 mHa for $2^1A_1$ (between DMC and CASPT3, closer to DMC) and 247.4 mHa for $3^1A_1$ (matching CASPT3). The paper states: "even the best current estimates... disagree by 0.1–0.15 eV, a range which the Psiformer is within. This demonstrates that NES-VMC is among the state of the art in challenging excited state calculations, where even other top methods disagree by more than chemical accuracy."

Singlet-targeting for tetrazine and cyclopentadienone (Section S5). Both systems required the spin-shift Hamiltonian $\hat{H}'_\lambda = \hat{H} + \lambda \hat{S}^2$ with $\lambda = 0.5$ (tetrazine) and $\lambda = 1.0$ (cyclopentadienone) to push triplet states above the singlets of interest. The paper states: "a large number of low-lying triplet states meant that it was impractical to compute enough states to compare against other results in the literature without some modification." The spin magnitudes in Table S9 confirm that all computed states are singlets ( $\langle \hat{S}^2 \rangle = 0.002$ –0.009 for the Psiformer, within expected Monte Carlo noise around zero for singlet states).

Pretraining sensitivity for double excitations (Section S7). The paper notes that "in all the experiments in Sec. VII we tried both random and ordered pretraining and reported the best results." For nitrosomethane, butadiene, and glyoxal, random pretraining was more robust; for tetrazine and cyclopentadienone (with spin-shifted Hamiltonians), ordered pretraining was more effective. On cyclopentadienone with the Psiformer, "training with ordered pretraining experienced numerical instability in the middle of training, but we found that random pretraining also worked well and included those results."

Benzene (Section VIII, Figure 6, Table S10)

Benzene (42 electrons, the largest molecule studied) serves as the capstone benchmark, with abundant reference data from coupled-cluster, CASPT2, TD-DFT, and prior VMC penalty methods. It is "the same size as tetrazine and cyclopentadienone" but benefits from more extensive literature for comparison.

Headline results. "NES-VMC with the Psiformer comes very close to reaching the TBE for all computed states." The Psiformer errors relative to the coupled-cluster TBEs from QUEST (Loos et al., 2020) are: $^3B_{1u}$ : 3.5 mHa ( $156.5$ vs. $153$ mHa TBE), $^3E_{1u}$ : 1.7–2.1 mHa ( $180.7$ – $181.1$ vs. $179$ mHa), $^1B_{2u}$ : 2.5 mHa ( $188.5$ vs. $186$ mHa), $^3B_{2u}$ : 2.8 mHa ( $216.3$ vs. $213.5$ mHa). All are within chemical accuracy (1.6 mHa = 43 meV), though the paper reports these in mHa rather than eV. (The $^3B_{1u}$ error of 3.5 mHa = 95 meV, which is above the 43 meV chemical accuracy threshold, but the paper does not flag this.)

FermiNet limitations. The FermiNet "is not quite as accurate, and struggles with the highest energy $^3B_{2u}$ state." Table S10 reveals that this state has $\langle \hat{S}^2 \rangle = 1.126$ with the FermiNet—far from the expected 2.0 for a pure triplet. The paper interprets this as "contamination by the singlet $^1B_{1u}$ state," which "would explain the anomalously high energy" (230.2 mHa vs. 213.5 mHa TBE). This is a qualitative failure of state separation, not merely a quantitative error.

Comparison to prior VMC penalty methods. The paper makes a direct comparison to two penalty method results:

Pathak et al. (2021): Slater-Jastrow ansatz with an unbiased penalty method, followed by DMC correction. Table S10 shows their DMC-corrected energies: $^3B_{1u}$ : 152(1) mHa, $^3E_{1u}$ : 180(1)–182(1) mHa, $^1B_{2u}$ : 189(1) mHa, $^3B_{2u}$ : 216(1) mHa. These are very close to the Psiformer NES-VMC results. The paper notes that "the results reported in Pathak et al. include a diffusion Monte Carlo correction which reduces the error by ~0.1 eV, while NES-VMC does not include any post-processing of VMC results."
Entwistle et al. (2023): PauliNet with a penalty method, primarily for the $^1A_{1g}$ and $^3B_{1u}$ states. The paper reports: "Relative to the TBE, the FermiNet and Psiformer errors on the $^1A_{1g} \to ^3B_{1u}$ transition are 0.126 and 0.088 eV, respectively. Without variance matching, the PauliNet penalty method error is 1.4 eV, and even with variance matching the error is only reduced to 0.29 eV, significantly worse than the Psiformer."

Comparison to traditional quantum chemistry. CASPT2 and TD-DFT-PBE0 results are plotted in Figure 6a. The paper states: "CASPT2 and TD-DFT methods are less accurate across the board, and CASPT2 is generally intermediate in accuracy between TD-DFT and coupled cluster." The specific TD-DFT errors are visible in the plot but not numerically tabulated in the main text.

Natural orbital analysis (Figure 6b). The paper computes density matrices and natural orbitals for all excited states. The analysis reveals: "All five excited states for benzene we computed are single excitations from a $\pi$ to $\pi^*$ orbital, but they are best described by exciting half an electron from two distinct $\pi_g$ orbitals into two distinct $\pi_u^*$ orbitals." The figure visualizes the two occupied $\pi$ orbitals (top row) and the two unoccupied $\pi^*$ orbitals (bottom row) involved in the excitations, providing physical insight into the multi-configurational character of the excited states.

Numerical stability and computational cost. For benzene, the paper used the bottleneck envelope (Section S4) with $|\ell| = 32$ distinct envelopes, reducing memory overhead compared to the default $K \times N_\text{det} \times (N_\alpha + N_\beta) = 6720$ envelopes. The FermiNet used 200,000 training iterations with blocks of 4 per electron set, while the Psiformer used 100,000 iterations with the same blocking. The ensemble penalty method was attempted on benzene but "became numerically unstable after a few hundred iterations" and produced no results — a significant negative result demonstrating NES-VMC's superior stability at this scale.

Ablation Studies and Robustness Checks

Number of computed states $K$ (atomic spectra, Section III, Table S2, Figure 1): Increasing $K$ from 5 to 10 reduces the error on the highest excited states of fluorine and neon by roughly half (FermiNet on F $^4P^\circ$ : 10.8 → 9.68 mHa; on Ne $^2[1/2]^\circ$ : 8.30 → 1.7 mHa). The advantage of larger $K$ is not merely better convergence of individual states but the ability to disentangle near-degenerate states that would otherwise be mixed in a smaller calculation. The Psiformer with $K = 10$ consistently outperforms the FermiNet with $K = 10$ on the most challenging cases (F: 1.3 vs. 2.7 mHa; Ne: 1.56 vs. 1.7 mHa). The paper does not explore $K > 10$ due to the cubic scaling becoming significant.

Neural network architecture: FermiNet vs. Psiformer (throughout, most systematically in Tables S9 and S10): On double excitation systems with >25 electrons, the Psiformer consistently outperforms the FermiNet, with MAEs of 21 vs. 15 (nitrosomethane, both within chemical accuracy), 38 vs. 84 (butadiene), 28 vs. 167 (glyoxal), 54 vs. 45 (tetrazine, FermiNet slightly better on this system), and 66 vs. 92 meV (cyclopentadienone). The largest gap is on glyoxal (167 vs. 28 meV), which the paper attributes to the FermiNet's known scaling limitations. On benzene, the Psiformer correctly resolves all states while the FermiNet shows singlet-triplet contamination on the highest $^3B_{2u}$ state ( $\langle \hat{S}^2 \rangle = 1.126$ vs. expected 2.0). On smaller molecules (Section IV, Table S3), the two architectures perform similarly. This confirms that the choice of ansatz matters increasingly as system size grows, but that NES-VMC works effectively with both architectures.

Pretraining strategy: ordered vs. random (Section S7, Figures S1, S2, Tables S3, S4): On HNO and HCF, ordered pretraining (matching Hartree-Fock excitations) causes both NES-VMC and the ensemble penalty method to miss specific states — the $1A'$ double excitation of HNO and the singlet state of HCF at 5.65 eV. Random pretraining (random combinations of valence Hartree-Fock orbitals) recovers these states with chemical accuracy (FermiNet errors: 0.2 mHa on HNO $1A'$ , -0.1 mHa on HCF $1A''$ ). The paper demonstrates via Figure S2 that the failure mode is the optimization getting stuck at a saddle point, with diagnostic signatures: plateaued energy, transiently increasing gradient norm, then sudden drop upon escape. The sensitivity is not unique to NES-VMC — the ensemble penalty method misses the same states with the same pretraining. On larger double excitation systems (Section VII), the same pattern holds: random pretraining is more robust for nitrosomethane, butadiene, and glyoxal, while ordered pretraining is more effective for tetrazine and cyclopentadienone (with spin-shifted Hamiltonians). The paper reports the best results from either pretraining strategy for each system.

Network size ablation on H $_2$ O (Section S11, Figure S4): Varying the number of layers (2–4), one-electron stream width (16–256), and either two-electron stream width for FermiNet (8–32) or attention head units for Psiformer (4–64) shows that the "width of the one-electron stream is the primary factor in convergence to the ground truth" for the FermiNet, while the Psiformer is "more robust than the FermiNet to reduction in the width of the intermediate layers." Even with heavily reduced networks (2 layers, 16 one-electron units), the FermiNet achieves near-chemical accuracy on H $_2$ O's first two excited states ( $^1A_2$ and $^1B_1$ ). This suggests that the default architecture is substantially overparameterized for small systems and that the residual errors reflect the optimization landscape rather than representational capacity.

Kinetic energy computation method: forward Laplacian vs. original (Section S11, Figure S3): Using Folx (Gao and Günnemann, 2023) for the forward Laplacian computation yields a roughly 4x speedup over the original column-by-column Hessian method, with "identical" numerical results "to within expected run-to-run variability." Figure S3 shows that with the forward Laplacian, the asymptotic scaling is closer to cubic in $K$ , while the original method shows slightly sub-cubic scaling. The cubic term does not dominate until $K > 10$ .

Ensemble penalty method on HNO and HCF (Section S7C, Figure S1, Table S4): This is a baseline implementation for controlled comparison. The ensemble penalty method achieves mean absolute errors of 3.2 eV on states identifiable by both methods, versus 0.25 eV for NES-VMC — an order of magnitude difference. The penalty method reaches chemical accuracy on only a single identifiable state and fails to find several states entirely. On benzene, the ensemble penalty method "became numerically unstable after a few hundred iterations" and could not be run. The paper attributes this to the numerical challenges of unbiased penalty methods with neural network ansätze, noting that clip values of 5 times the median variation were used for both energies and overlap matrix elements to maintain stability.

Pseudopotential usage for heavy atoms (Section S9, Table S3, Section IV): For molecules with second-row atoms (HCl, H $_2$ S, CH $_2$ S, H $_2$ CSi), replacing core electrons with ccECP pseudopotentials from PySCF "significantly improved the accuracy of the results, likely because the total energy scale was reduced by ignoring core electrons." The paper does not provide a direct head-to-head comparison of all-electron vs. pseudopotential results, but the achieved accuracy with pseudopotentials (chemical accuracy on HCl, H $_2$ S) supports the claim.

Median-of-means estimator for large systems (Section S8, stated in text): For systems with >25 electrons (butadiene, glyoxal, tetrazine, cyclopentadienone, benzene), the distribution of local energies became heavy-tailed, making the empirical mean unstable with errors "on the order of several mHa or ~0.1 eV." Switching to the median over minibatches of 16 walkers, then taking the mean over these medians, "improved the accuracy for the majority of states of large systems." The paper does not quantify the magnitude of improvement or show results with and without this estimator.

Spin magnitude as a diagnostic (Figures S5, Table S3, S8–S10): Throughout the paper, $\langle \hat{S}^2 \rangle$ is reported for every computed state to verify correct spin assignment. Figure S5 shows the full matrix of $\langle \hat{S}^2 \rangle$ for the first 5 excited states of ethylene, showing two cleanly identified triplet states (diagonal values near 2.0) and near-zero off-diagonal elements, confirming effective state separation. In Table S8 for twisted ethylene, the $\langle \hat{S}^2 \rangle$ values correctly track the singlet (0.001–0.009) and triplet (1.993–2.002) states across all torsion and pyramidalization angles, including through the crossing at $\tau = 90°$ where the ground state becomes triplet ( $\langle \hat{S}^2 \rangle = 1.998–1.999$ ) and the excited singlet has $\langle \hat{S}^2 \rangle = 0.001–0.003$ . This is strong evidence that the method correctly handles near-degeneracies without root-flipping or state mixing.

Critical Assessment

The experiments collectively demonstrate that NES-VMC is a practically effective method for computing excited states with neural network ansätze, achieving accuracy comparable to or exceeding state-of-the-art quantum chemistry methods on a wide range of systems. However, a critical reading reveals several important qualifications, unaddressed questions, and limitations in what the experiments actually establish.

What the Experiments Genuinely Demonstrate

The paper's strongest and most robustly supported claim is that NES-VMC eliminates the state collapse problem that plagues penalty methods, and that this translates to dramatically better accuracy when using neural network ansätze. The evidence is multi-layered: (1) the direct comparison on benzene (Psiformer error 0.088 eV vs. PauliNet penalty 1.4 eV uncorrected, 0.29 eV corrected), (2) the controlled comparison on HNO and HCF where NES-VMC achieves 0.25 eV MAE vs. 3.2 eV for the ensemble penalty method implemented on the same architectures, (3) the ensemble penalty's complete failure (numerical instability) on benzene where NES-VMC succeeds, and (4) the recovery of the lithium $^2S$ and $^2P^\circ$ states that the PauliNet penalty method missed entirely. This is a clear, quantitative demonstration of superiority over the most directly comparable alternative, and the magnitude of the difference (factors of 3–10x in error) strongly supports the claim that the variational principle, not the ansatz, is the critical factor.

The chemical accuracy on most molecules is well-supported. For the 13 small molecules in Section IV, almost all vertical excitation energies fall within 1.6 mHa of the QUEST TBEs (Table S3). The oscillator strengths match within ~0.001. For the double excitation systems, the Psiformer achieves MAEs of 21–66 meV, and is within or near chemical accuracy for all five systems. For benzene, the Psiformer errors are 1.7–3.5 mHa relative to the TBEs. These numbers are credible and well-documented. However, "chemical accuracy" as a bright line is not uniformly achieved: the Psiformer error on tetrazine (54 meV) and cyclopentadienone (66 meV) exceeds 43 meV, and the $^3B_{1u}$ state of benzene shows 3.5 mHa = 95 meV, which the paper does not flag as exceeding the threshold. This is a minor inconsistency in the narrative but does not undermine the overall accuracy.

The correct treatment of challenging physical phenomena—double excitations, conical intersections, avoided crossings—is convincingly demonstrated. The carbon dimer analysis (Section V) shows NES-VMC correctly identifying the double excitation character of the $b^3\Sigma_g^-$ and $B^1\Delta_g$ states via natural orbital occupancies, achieving adiabatic energies within 7 meV of experiment (improving on SHCI by ~4×). The twisted ethylene analysis (Section VI) shows the method navigating the $N$ / $V$ state crossing at the correct pyramidalization angle, closer to MR-CI than the PauliNet penalty method. The butadiene analysis (Section VII) correctly orders the historically controversial $1^1B_u$ and $2^1A_g$ states. These are not merely accurate numbers—they demonstrate that the method captures the correct physics, which is a stronger validation than energetic accuracy alone.

Weaknesses and Unaddressed Questions

Single reference method, no basis-set limit convergence. This is the most significant structural limitation of the experimental validation. All accuracy claims are relative to other approximate methods (coupled-cluster TBEs, SHCI, MR-CI, CASPT2), not to the exact solution of the Schrödinger equation. For small systems (atoms, diatomics), the reference data is highly accurate (experimental for atoms, extrapolated coupled-cluster for molecules), but for larger systems, the reference methods themselves have uncertainties. This is most apparent for cyclopentadienone, where the two best reference calculations (DMC and CASPT3) disagree by 4–6 mHa on the double excitation energies—a range that NES-VMC falls within, but which itself represents uncertainty in the "ground truth." The paper cannot claim to be more accurate than the reference methods when the reference methods disagree. A more rigorous validation would include full configuration interaction (FCI) comparisons for small basis sets where the exact answer is known, but such comparisons are absent.

The basis set limit is not addressed. NES-VMC operates in real space without a basis set, so in principle it can approach the complete basis set limit. However, the neural network ansätze use multiplicative envelopes with exponential decay, which may not be flexible enough to capture all correlation effects. The paper does not investigate convergence with envelope complexity, number of determinants, or network size beyond the simple ablation on H $_2$ O (Fig. S4), which suggests overparameterization for small systems but does not address whether the default architecture is sufficient for larger ones. The drift relative to SHCI of 5 mHa/Å on the carbon dimer is attributed to basis set effects in SHCI, but this is an inference, not a demonstrated fact—the NES-VMC results could themselves have systematic errors at stretched geometries.

The $K$ -dependence of accuracy is demonstrated but not systematically characterized. The atomic spectra results (Fig. 1) show that increasing $K$ from 5 to 10 improves accuracy for the highest states of fluorine and neon, but the paper does not investigate at what $K$ the results would converge. For systems like neon with many closely-spaced Rydberg states, it is possible that even $K = 10$ is insufficient to fully disentangle the spectrum. The cubic scaling in $K$ (Fig. S3) means that computing very large numbers of states will become expensive, but the paper does not establish how many states are "enough" as a function of system size or density of states. This is a practical concern for applying the method to systems with dense manifolds of excited states (e.g., transition metal complexes, large conjugated molecules).

The pretraining sensitivity is a genuine limitation, not merely an implementation detail. The fact that both NES-VMC and the ensemble penalty method miss the same states with ordered pretraining (Section S7C) is presented as evidence that "sensitivity to initialization is a general feature," which is true but does not make it less of a problem. The paper's solution—try both ordered and random pretraining and report the best—is a post-hoc selection strategy that would not be available in a truly predictive setting where the correct answer is unknown. The saddle-point behavior documented in Fig. S2 is a real failure mode of the optimization, and the paper does not provide a reliable diagnostic for detecting it during training (the gradient norm increase is transient and may be missed). This means that in a production setting, a user might converge to a saddle point, observe stable energies and gradients, and incorrectly conclude the calculation is converged. The paper acknowledges this implicitly by noting that for the double excitation systems, "we tried both random and ordered pretraining and reported the best results"—but this is cherry-picking, not a robust protocol.

The singlet-targeting trick with $\lambda \hat{S}^2$ is effective but has an unexamined cost. For tetrazine and cyclopentadienone, adding $\lambda \hat{S}^2$ to the Hamiltonian with $\lambda = 0.5$ –1.0 successfully pushes triplet states above the singlets of interest. However, the paper does not investigate whether this modification affects the accuracy of the singlet excitation energies. In principle, because $\hat{H}$ and $\hat{S}^2$ commute, the singlet eigenfunctions are unchanged, but in practice, the optimization landscape is modified (the energy gaps between states are different, which could affect convergence rates or saddle-point behavior). The singlet-targeting approach also means that the method cannot simultaneously compute singlet and triplet states, losing one of the advantages of the multi-state framework. For systems where both singlet and triplet excitations are of interest, this would require two separate calculations.

The comparison to prior VMC penalty methods is somewhat unfair. The Entwistle et al. (2023) results used the PauliNet ansatz, which is a different architecture from the FermiNet/Psiformer. The relative contributions of the ansatz vs. the variational principle cannot be cleanly separated without running the penalty method on the FermiNet/Psiformer. The paper does this for HNO and HCF (Fig. S1, Table S4), and the penalty method performs poorly there, but this is only on two small molecules. For benzene, the paper states the ensemble penalty method became numerically unstable and could not be run. This leaves open the possibility that the penalty method could be made to work with better numerical stabilization (different clipping, different optimizer, different penalty formulation), and the poor performance in Fig. S1 reflects suboptimal implementation rather than fundamental limitations of penalty methods. The paper does not explore the penalty method's hyperparameter sensitivity (penalty strength $\lambda$ , state weights $w_i$ ), which Wheeler et al. (2024) showed has a critical threshold—it is possible the single choice $\lambda = 1$ , $w_i = i^{-1}$ is far from optimal for these systems.

Scalability beyond benzene is not demonstrated. The largest molecule studied is benzene (42 electrons, 12 atoms). While this is a meaningful scale for a wavefunction method, it is far from the regime where NES-VMC would compete with conventional excited-state methods (TD-DFT, EOM-CCSD, CASPT2), which can routinely handle molecules with 50–100+ atoms. The paper does not establish how the method's accuracy degrades (or doesn't) with system size, nor does it provide timings or resource estimates that would allow a reader to assess whether running NES-VMC on a molecule of interest is feasible. The bottleneck envelope (Section S4) is introduced specifically to handle benzene-scale memory requirements, suggesting that scaling to significantly larger systems may require further architectural innovations.

The oscillator strength on HCl with the Psiformer is a genuine failure, not explained. The Psiformer finds the $^1\Pi$ state of HCl to be essentially dark ( $f \approx 0$ ) while the FermiNet correctly predicts $f \approx 0.022$ (matching the TBE). This is a qualitative error—the method predicts a bright state to be dark, which would be misleading if interpreted physically. The paper notes this anomaly in passing but does not investigate its cause. Is it a convergence issue? A specific failure of the Psiformer architecture for this system? A manifestation of the $\sigma_i$ non-identifiability causing an incorrect sign cancellation in the off-diagonal observable recovery? This last possibility is concerning: the paper's derivation (Section S2, Eq. 51) assumes that elementwise multiplication by the transpose cancels the $\sigma$ factors, but this requires that the off-diagonal matrix element be purely real and that the sign ambiguity from $\sigma_i/\sigma_j$ not interact with the vector nature of the dipole moment in a way that causes cancellation. If the Psiformer assigns different $\sigma$ factors to different Cartesian components of the transition dipole, the elementwise product $\Delta_k \circ \Delta_k^T$ might not correctly reconstruct the squared magnitude. The paper does not address this.

Statistical error bars are reported but not discussed in interpretation. All energy values in Tables S2–S10 include standard errors in the last digit(s). Some of these error bars are substantial relative to the claimed accuracy: for example, the FermiNet error on F $^4P^\circ$ is 9.68(7) mHa, meaning the 95% confidence interval is ~9.68 ± 0.14 mHa—still clearly above chemical accuracy (1.6 mHa). But for borderline cases like the Psiformer on tetrazine $^2A_g$ (error 16.9 mHa beyond the unsafe TBE, but within ~1 mHa of the DMC value of 183 mHa), the statistical uncertainty could shift the conclusion about whether the result agrees with one reference or another. The paper's variance propagation through diagonalization (Section S12) is mathematically careful, but the resulting uncertainties are not used to qualify any of the accuracy claims.

What Would Strengthen the Paper

Several experiments are conspicuous by their absence:

FCI comparisons in a minimal basis set. For a small molecule (e.g., H $_2$ O, BH) in an STO-3G or 6-31G basis, full configuration interaction provides the exact answer within that basis. Comparing NES-VMC to FCI would establish the method's accuracy against a known ground truth, removing the ambiguity of comparing against other approximate methods. This is a standard validation in the quantum chemistry literature and its absence is notable.
Systematic study of $K$ -convergence. How does the accuracy of the $k$ -th state depend on the total number of computed states $K$ ? The atomic results suggest improvement with larger $K$ , but a systematic study on a system with a dense spectrum (e.g., neon or a transition metal atom) would establish practical guidelines for choosing $K$ .
Direct comparison to EOM-CCSD on the same geometries. Equation-of-motion coupled-cluster is the workhorse for excited states of medium-sized molecules and is the source of most QUEST TBEs. A head-to-head comparison of NES-VMC with EOM-CCSD on a subset of the QUEST molecules, reporting both energies and timings, would help readers assess the practical trade-offs.
Ablation on the number of determinants. The paper uses 16 determinants throughout. How does accuracy depend on this number for excited states? For ground states, multi-determinant expansions capture static correlation. For excited states, especially double excitations, the number of determinants might be even more critical. This ablation is absent.
Timing and resource estimates. The paper reports iteration counts (100,000–200,000) and GPU usage (A100) for the scaling analysis but does not provide wall-clock times or total GPU-hours for any of the production calculations. A reader cannot assess whether the method is computationally competitive with, say, running EOM-CCSD on a CPU cluster.

Conditional Claims and Their Boundaries

The paper's central claim—that NES-VMC "has no free parameters, has convergence guarantees, and is applicable to all possible Ansätze and all excited states, regardless of symmetry"—holds for the method's mathematical formulation, but the experimental validation reveals practical boundary conditions:

"No free parameters" is true of the method, but pretraining strategy becomes a de facto free parameter. The choice between ordered and random pretraining is binary but requires trial-and-error (or prior knowledge of which states to expect) to get right. The paper's solution of trying both and reporting the best is not parameter-free in practice.
"Convergence guarantees" apply to local minima of the total energy, not to the correct physical states. The saddle-point behavior (Fig. S2) shows that the optimization can converge to a local minimum of the objective that does not correspond to the correct set of excited states (it misses double excitations). The Robbins-Monro convergence guarantee is for convergence to a local minimum, not the global minimum, and the paper provides no guarantee or diagnostic for distinguishing the two.
"Applicable to all possible Ansätze" is demonstrated for two closely related neural network architectures. Both the FermiNet and Psiformer are continuous real-space wavefunctions with Slater determinant structure. The method should in principle work with any ansatz, but the paper does not test it on, e.g., autoregressive models, Gaussian basis functions, or tensor network states. The determinant construction requires evaluating the ansatz at arbitrary particle configurations $x^j$ , which may not be natural for all ansatz classes.
The $4\times$ improvement over SHCI on carbon dimer adiabatic energies is a comparison against a single reference method on a single system. This should not be over-interpreted as a general accuracy advantage over SHCI or other selected CI methods. The comparison is suggestive but isolated.

Overall, the experimental section convincingly establishes NES-VMC as a practically useful method that enables neural network ansätze to achieve state-of-the-art accuracy on excited states where they previously failed. The validation is thorough in its breadth (atoms through benzene, multiple physical phenomena) but has gaps in its depth (no FCI comparisons, limited $K$ -convergence analysis, no timing data). The method's limitations—pretraining sensitivity, the ~38% correct-to-incorrect reversion problem in ensemble penalty methods (not applicable to NES-VMC but indicative of the broader difficulty of the problem), and the scaling beyond benzene—are honestly acknowledged but not solved. The paper makes a compelling case that the variational principle, not the ansatz, was the bottleneck, and that NES-VMC is the right principle for the job.

6. Limitations and Trade-offs

The Difficulty Estimation Cost Is Not Accounted For

In the machine learning domain where this paper sits, the question of difficulty estimation setup costs applies directly to the computational budget. NES-VMC requires generating enough walker configurations to reliably accumulate and diagonalize the local energy matrix, and the optimization must run long enough for the states to converge to the correct subspace.

The computational cost is significant and poorly characterized. The paper reports iteration counts (100,000–200,000 training steps; Table S1) and notes that "relative energies were almost always well converged after 100,000 iterations" while "optimal absolute energies" require 200,000. However, these numbers are not accompanied by wall-clock times, total GPU-hours, or any cost comparison against the reference methods (EOM-CCSD, SHCI, CASPT2, DMC) that NES-VMC is benchmarked against. A practitioner cannot determine whether the accuracy gains—for example, the 4× improvement in adiabatic energy error on carbon dimer relative to SHCI—justify the computational cost relative to simply running a larger SHCI calculation or a different method.

What evidence exists: The scaling analysis in Figure S3 shows time-per-iteration on a single neon atom as a function of the number of states $K$ , but neon (10 electrons) is the smallest system studied. The scaling with system size $N$ —which is $O(N^3)$ to $O(N^4)$ for neural network VMC—is not benchmarked. The paper states that using the forward Laplacian (Folx) provides "a 2–5× speedup depending on the system," but does not report absolute timings for any system. The bottleneck envelope (Section S4) is introduced for benzene because "the memory overhead of the default FermiNet and Psiformer became prohibitive," indicating that resource constraints were binding even at the 42-electron scale. The batch sizes (64 for scaling analysis, Table S1 shows block sizes of 1–4 per electron set for larger systems) and GPU type (A100) are mentioned but not combined with iteration counts to produce total resource estimates.

Mitigation status: The paper does not address this limitation. There is no cost-versus-accuracy analysis, no comparison of total core-hours or GPU-hours against reference methods, and no discussion of how the method's cost scales to systems larger than benzene. A practitioner evaluating whether to use NES-VMC on a molecule of interest has no way to estimate the computational resources required.

The Pretraining Strategy Is a De Facto Free Parameter That Determines Success or Failure

The paper's headline claim is that NES-VMC "has no free parameters" (Section I, Abstract). This is true of the variational principle itself, but the experimental results reveal that the choice of pretraining strategy—ordered versus random initialization of the neural network orbitals—determines whether the method finds the correct excited states or converges to a saddle point that misses them entirely.

The consequence is a qualitative failure mode, not merely a quantitative inaccuracy. On nitroxyl (HNO) with ordered pretraining, both NES-VMC and the ensemble penalty method miss the $1A'$ double excitation at 4.32 eV (Section S7C, Figure S1, Tables S3–S4). On fluoromethylene (HCF), ordered pretraining misses the singlet state at 5.65 eV. On the double excitation systems (Section VII), the paper explicitly states: "in all the experiments in Sec. VII we tried both random and ordered pretraining and reported the best results." This is a post-hoc selection that would not be available in a truly predictive setting where the correct states are unknown. The paper provides no reliable diagnostic for detecting whether the optimization is stuck at a saddle point—Figure S2 shows that the gradient norm increases transiently but this may be missed, and the energy plateaus could be mistaken for convergence.

What evidence exists: Figure S1 and Tables S3–S4 demonstrate the sensitivity conclusively. The FermiNet with ordered pretraining finds no identifiable $1A'$ state on HNO; with random pretraining, it finds it at 4.32 eV with 0.2 mHa error. On HCF with ordered pretraining, the FermiNet finds "only a triplet at 6.77 eV" instead of the singlet; with random pretraining, it finds the singlet at 5.75 eV. The paper's Figure S2 provides diagnostic evidence that the failure is a saddle point: the energy plateaus for thousands of iterations, the gradient norm increases transiently, and then a sudden drop occurs when the optimization escapes. Critically, the paper states that "this sensitivity to choice of pretraining does not seem to be unique to NES-VMC, and seems to affect penalty methods as well," which is true but does not make the limitation less severe—it means the problem is fundamental to the optimization landscape, not an artifact of this particular method.

Mitigation status: The paper's only mitigation is the ad-hoc strategy of trying both pretraining methods and keeping the better result. This is heuristic, not principled. The paper does not develop a systematic initialization protocol, an automated saddle-point detection method, or a guarantee that random pretraining will always find states that ordered pretraining misses (or vice versa). The statement that "we tried both random and ordered pretraining and reported the best results" (Section S7B) is an honest disclosure but does not constitute a solution. For practitioners, this means that a single NES-VMC calculation cannot be trusted to have found the correct states without either prior knowledge of the expected spectrum or multiple runs with different initializations.

The Method Does Not Scale Beyond Benzene, and the Scaling Behavior Is Not Characterized

The largest molecule studied in the paper is benzene (42 electrons, 12 atoms). While this is a meaningful scale for a wavefunction-based method, it is far from the regime where NES-VMC would compete with conventional excited-state methods like TD-DFT, EOM-CCSD, or CASPT2, which routinely handle molecules with 50–100+ atoms. The paper does not demonstrate or characterize performance on any system larger than benzene.

The consequence is that a practitioner cannot assess whether NES-VMC is feasible for their molecule of interest. The bottlenecks are not merely computational speed but memory: the paper introduces the bottleneck envelope (Section S4) specifically because "for the experiments with benzene, the memory overhead of the default FermiNet and Psiformer became prohibitive." This was at 42 electrons with 10 states, using 16 determinants, requiring 6720 orbitals in the default envelope ( $K \times N_\text{det} \times N_\text{electrons} = 10 \times 16 \times 42$ ). The bottleneck reduces this to 32 distinct envelope functions, but the paper does not report what memory this consumed, whether further reductions are possible, or what the next scaling bottleneck would be.

What evidence exists: The lack of evidence is itself the limitation. Figure S3 shows scaling with the number of states $K$ on a single neon atom (10 electrons) with a batch size of 64 on an A100 GPU—"much smaller scale than the experiments elsewhere in the paper," as the paper acknowledges. The scaling with system size $N$ is stated to be $O(N^3)$ to $O(N^4)$ for neural network VMC (Section I), but this is not benchmarked for NES-VMC specifically. The paper does not report wall-clock times, memory consumption, or GPU-hours for any of the production calculations (atoms through benzene). Table S1 shows that training iterations range from 100,000 to 200,000, but without per-iteration timings, this conveys no information about total cost. The ensemble penalty method was attempted on benzene but "became numerically unstable after a few hundred iterations" (Section S7C), suggesting that benzene is near the practical limit of what current implementations can handle, but the paper does not diagnose whether this limit is algorithmic (the determinant construction) or architectural (the neural network size).

Mitigation status: The paper acknowledges the scaling concern implicitly by introducing the bottleneck envelope, but does not frame scaling as a limitation or discuss it in the conclusions. Section IX notes that "classic Ansätze like the Slater-Jastrow Ansatz can be scaled to much larger systems," suggesting that NES-VMC could be applied with cheaper ansätze, but this is speculation without experimental support. The paper also notes that "as the system size grows, applying different weighting to the different states is often needed to get good performance with existing QMC methods. It remains to be seen if the same holds true for NES-VMC." This is an honest acknowledgment of an open question, not a characterization of the method's current capabilities.

The Local Energy Matrix Diagonalization Can Fail to Correctly Assign States, and the Failure Is Not Detectable Without Ground Truth

The state demixing procedure—diagonalizing the accumulated local energy matrix $\mathbb{E}[E_L] = U\Lambda U^{-1}$ to recover individual state energies and the transformation $U$ for observables—is mathematically sound at the global minimum of the total energy. However, the paper documents several cases where this procedure produces physically incorrect results, and there is no internal diagnostic that would flag these failures.

The consequence is the possibility of silent errors in state assignment or observable prediction. The most striking example is the Psiformer on HCl (Section IV, Table S3): the FermiNet correctly computes the $^1\Pi$ oscillator strength as 0.0223 (matching the TBE of 0.0231 to within ~0.001), but the Psiformer finds this state to be nearly dark, with an oscillator strength of −9.78×10 $^{-5}$ . The vertical excitation energies from both ansätze are nearly identical (−0.76 mHa vs. −0.85 mHa error), so the energy prediction is reliable but the observable prediction is qualitatively wrong. The paper notes this anomaly but does not investigate its cause. Possible explanations include: the Psiformer state demixing assigned the wrong eigenvector to the $^1\Pi$ state, the non-identifiable scale factors $\sigma_i$ in $U = \Sigma A^{-1}$ interacted with the vector nature of the dipole moment to cause cancellation that the elementwise transpose multiplication $\Delta_k \circ \Delta_k^T$ did not resolve, or the Psiformer converged to a different linear combination of states spanning the same subspace.

On boron (Table S2), the Psiformer with $K = 10$ converges the highest state to "the $2P^\circ$ state rather than the last $2D$ state"—a misassignment of the physical state despite correct energies for the lower states. On the FermiNet for benzene (Table S10), the highest state shows $\langle \hat{S}^2 \rangle = 1.126$ —far from the expected 2.0 for a pure triplet—which the paper interprets as "contamination by the singlet $^1B_{1u}$ state." This means the state assignment is ambiguous and the energy (230.2 mHa) is unreliable.

What evidence exists: The HCl oscillator strength failure (Table S3), the boron state misassignment (Table S2), and the benzene FermiNet contamination (Table S10) are all documented in the paper. However, there is no systematic investigation of when or why the diagonalization fails. The paper's derivation shows that $U = \Sigma A^{-1}$ at the exact minimum, but at a finite level of convergence, the eigenvectors $U$ may mix states, especially when states are near-degenerate. The $\sigma_i$ non-identifiability is handled for oscillator strengths by elementwise multiplication with the transpose (Section S2), but the paper does not prove that this always yields the correct squared magnitude when the transformation $U$ is imperfect.

Mitigation status: The paper does not address these failures beyond noting them. There is no diagnostic for detecting when the state assignment is unreliable (e.g., examining the off-diagonal elements of the demixed observable matrices, checking that the demixed $\hat{S}^2$ matrix is approximately diagonal, or monitoring the stability of the eigenvectors during the latter stages of optimization). The spin magnitude $\langle \hat{S}^2 \rangle$ is reported alongside all energies and serves as a partial diagnostic—a value far from the expected integer $S(S+1)$ signals contamination—but this only works for spin and not for spatial symmetries or state ordering. No analogous diagnostic exists for transition dipole moments or other off-diagonal observables.

There Is No Demonstration of Simultaneous Singlet and Triplet Computation for Larger Systems

One of the proposed advantages of NES-VMC is that it finds "the lowest $K$ eigenfunctions... regardless of symmetry" (Section I), computing all states in a single optimization. However, for the two largest double-excitation systems—tetrazine and cyclopentadienone, both the same size as benzene—the paper cannot use this capability and must instead modify the Hamiltonian to artificially push triplet states to higher energies.

The consequence is that for systems with many low-lying states of different spin, NES-VMC as demonstrated cannot simultaneously compute singlets and triplets. The singlet-targeting modification (Section S5) uses $\hat{H}'_\lambda = \hat{H} + \lambda \hat{S}^2$ with $\lambda = 0.5$ (tetrazine) or $\lambda = 1.0$ (cyclopentadienone), chosen to be "roughly double the required scale to remove triplets." The paper states: "a large number of low-lying triplet states meant that it was impractical to compute enough states to compare against other results in the literature without some modification." This is a practical admission that the $K$ required to resolve the singlet states of interest—among the dense forest of triplets—would be too large to be computationally feasible or would push the scaling into the cubic regime where the method becomes expensive.

What evidence exists: The paper's own statement quoted above is the primary evidence. Table S9 shows that for tetrazine and cyclopentadienone, only singlet states are computed ( $\langle \hat{S}^2 \rangle = 0.002$ –0.009 for the Psiformer, essentially zero). The spin magnitudes confirm the singlet-targeting worked as intended. However, the paper does not report how many states would be needed to find these singlets without the $\lambda \hat{S}^2$ shift, what the computational cost would be, or whether the singlet-targeting modification affects the accuracy of the singlet excitation energies. In principle, because $\hat{H}$ and $\hat{S}^2$ commute, the eigenfunctions are unchanged and only the eigenvalues shift, but in practice, the optimization landscape changes—the energy gaps between states are different, which could affect convergence rates or the likelihood of getting stuck at saddle points. The paper notes that "for cyclopentadienone with the Psiformer, training with ordered pretraining experienced numerical instability in the middle of training" in the singlet-targeted calculation, suggesting the modification may have affected the optimization dynamics.

Mitigation status: The paper presents singlet-targeting as a practical workaround rather than a fundamental solution. The problem it addresses—too many low-lying triplet states to practically include in the calculation—is structural: as molecules grow larger, the density of excited states increases, and the number of states $K$ required to reach the states of interest grows. The cubic scaling with $K$ (Figure S3) means this becomes expensive, and the pretraining sensitivity means that a larger $K$ may introduce more saddle points. The paper does not explore alternatives, such as using spatial symmetry to separate states, or developing a version of NES-VMC that can target specific spin multiplicities without modifying the Hamiltonian. The singlet-targeting trick also means that a single calculation can only produce states of one spin multiplicity—defeating part of the purpose of the multi-state framework.

The Chemical Accuracy Claims Are Relative to Approximate Reference Methods, Not Exact Solutions

The paper's primary accuracy metric—comparison to "theoretical best estimates" (TBEs) from the QUEST database—benchmarks NES-VMC against other approximate methods (primarily coupled-cluster with basis set extrapolations), not against exact solutions of the Schrödinger equation. For larger systems where reference methods themselves have significant uncertainties, it is impossible to determine whether residual errors come from NES-VMC or from the reference.

The consequence is that the paper's accuracy claims cannot be independently verified for the most challenging systems. For cyclopentadienone (Section VII, Table S9), the two best reference calculations—diffusion Monte Carlo (Shepard et al., 2022) and CASPT3-corrected QUEST (Kossoski et al., 2024)—disagree with each other by 4–6 mHa (~0.1–0.15 eV) on the $2^1A_1$ and $3^1A_1$ excitation energies. NES-VMC with the Psiformer falls within this range (217.1 vs. 217 and 213 mHa; 247.4 vs. 253 and 247 mHa). The paper states that "even the best current estimates... disagree by more than chemical accuracy," which is true but means that the "ground truth" is uncertain to within ~0.1–0.15 eV. When NES-VMC agrees with one reference and disagrees with the other, it is unclear which method is more accurate. The claim that "NES-VMC is among the state of the art" is supported, but the claim that it achieves "chemical accuracy" (43 meV) on these systems cannot be verified because the reference itself is less accurate than chemical accuracy.

What evidence exists: The paper is transparent about the uncertainty in the reference data for tetrazine and cyclopentadienone, noting that "the original TBEs from QUEST were known to be unsafe for double excitations" and flagging these values in gray in Figure 5 and Table S9. However, even the "safe" QUEST TBEs are coupled-cluster results extrapolated to the complete basis set limit, which have their own error bars (typically estimated at ~0.05–0.1 eV for single excitations, larger for double excitations). The paper never compares NES-VMC to full configuration interaction (FCI) in a small basis set—the standard validation in quantum chemistry where the exact answer is known. For carbon dimer, the NES-VMC adiabatic energies are compared to experiment (Martin, 1992), which is a gold-standard reference, but this is only one system and the experimental uncertainties are not discussed.

Mitigation status: The paper does not address this limitation directly. The comparison to experiment for atoms (Section III, Table S2) and the carbon dimer adiabatic energies (Section V, Table S7) provides the strongest validation, since experimental data has well-characterized uncertainties (though the paper does not report them). For molecules, the reliance on TBEs from QUEST is a reasonable practical choice—these are the most accurate computational results available—but the paper does not discuss the uncertainty in the TBEs or attempt to quantify how much of the residual error could be due to reference method errors versus NES-VMC errors. A systematic comparison to FCI for small basis sets on small molecules would close this gap but is absent.

7. Implications and Future Directions

How This Work Changes the Landscape

This paper shifts the conversation around excited-state VMC from a constrained optimization mindset—where preventing state collapse requires externally imposed penalties, orthogonalization, or prior symmetry knowledge—to a structural parameterization mindset, where the algebraic form of the ansatz itself makes collapse impossible. The determinant construction $\Psi = \det(\psi_i(x^j))$ is not merely another technique to add to the toolbox; it is a conceptually different category of solution that eliminates the central pathology (state collapse) at the mathematical level rather than fighting it during optimization.

The magnitude of this shift is best understood by what it renders obsolete. The paper demonstrates that penalty methods—the dominant approach for neural network excited-state VMC prior to this work—achieve mean absolute errors an order of magnitude larger than NES-VMC on the same architectures (3.2 eV vs. 0.25 eV on HNO and HCF, Section S7C), miss states entirely that NES-VMC finds with proper pretraining, and become numerically unstable on systems as large as benzene. These are not incremental accuracy differences; they are qualitative failures. The paper's ensemble penalty implementation—following the most recent theoretically-grounded formulation with correct critical threshold conditions (Wheeler et al., 2024)—cannot be made to work reliably at the scale where neural network ansätze are most needed. This strongly suggests that the penalty paradigm is a dead end for deep learning wavefunctions, and that future work on neural network excited states should adopt structural collapse-prevention mechanisms (of which the determinant construction is the first demonstrated example) rather than attempting to tune penalty strengths.

The paper also reconciles a puzzling contradiction in the literature. Neural network ansätze like the FermiNet had achieved remarkable accuracy on ground states—exceeding CCSD(T) on some bond-stretching systems—yet had "not matched the accuracy of other methods for excited states calculations as they scale to larger systems, despite initial forays." The paper explicitly posed the question: "This could be due to the Ansatz not being sufficiently expressive, but it could also be due to the choice of variational principle." The experimental results decisively answer: it is the variational principle, not the ansatz. The evidence is multi-layered: (1) the PauliNet penalty method achieves 1.4 eV error on the benzene $^3B_{1u}$ state while NES-VMC with the Psiformer achieves 0.088 eV—a ~16× improvement using a different method on the same molecule; (2) the FermiNet and Psiformer produce similar results when both succeed, suggesting neither architecture is the bottleneck; (3) the network size ablation on H $_2$ O (Figure S4) shows that even heavily reduced architectures achieve near-chemical accuracy, implying the default networks are substantially overparameterized for small systems. This redirects research effort: the field should not prioritize designing more expressive neural network architectures for excited states, but should instead focus on (a) adopting variational principles with structural collapse-prevention, and (b) addressing the shared challenge of initialization sensitivity that affects all methods.

The identification of initialization sensitivity as a shared, fundamental challenge is itself a landscape-changing diagnostic. The paper shows that both NES-VMC and the ensemble penalty method miss the same states with the same pretraining strategy (Figure S1), and that the failure mode is convergence to a saddle point (Figure S2) rather than any method-specific pathology. This converts initialization from an implementation detail into a first-class research problem. Prior work had treated initialization heuristically—match Hartree-Fock orbitals and hope. The paper's evidence that saddle points systematically trap the optimization on states with double excitation character (HNO, nitrosomethane, butadiene, glyoxal) suggests a structural feature of the energy landscape: the subspace containing double excitations may be separated from the Hartree-Fock initialization by a barrier that gradient descent struggles to cross. Understanding the topology of this landscape—why saddle points form, how their prevalence scales with system size and the number of computed states $K$ , and whether alternative optimization algorithms (natural gradient, saddle-point-aware methods) can escape them more reliably—becomes a central research direction enabled by the paper's careful documentation of the phenomenon.

Finally, the paper establishes that neural network VMC can be competitive with state-of-the-art quantum chemistry for excited states at scales up to benzene. This was not obvious before this work. The PauliNet penalty method had reached ~0.29 eV error on the benzene $^3B_{1u}$ state (with variance matching corrections), which is far from competitive with coupled-cluster (~0.05–0.1 eV) or even CASPT2 (~0.1–0.2 eV). NES-VMC with the Psiformer reaches ~0.09 eV on the same state and is competitive with DMC-corrected Slater-Jastrow results (Pathak et al., 2021) without any post-processing. On double excitations, NES-VMC resolves discrepancies as large as ~1 eV in the QUEST database (tetrazine $2^1A_g$ ) and produces results within the uncertainty range of the two best competing methods on cyclopentadienone. This establishes neural network VMC as a credible alternative to projector QMC methods (DMC, AFQMC) for excited states—an alternative that provides explicit wavefunction forms enabling downstream analysis (natural orbitals, density matrices) that projector methods cannot provide.

Follow-Up Research This Work Enables

Systematic characterization of the initialization landscape for excited-state optimization. The paper documents that saddle points trap optimization and that different pretraining strategies (ordered vs. random) find different states, but does not map the landscape. A systematic study would: initialize many copies of the same architecture with different random seeds and both pretraining strategies, run NES-VMC to convergence on a fixed set of molecules (including HNO, butadiene, and glyoxal where the paper observed saddle-point behavior), record which states are found and at what energies, and cluster the results to identify how many distinct local minima exist, what states they correspond to, and whether the number of minima grows with system size or the number of computed states $K$ . This would reveal whether saddle points are isolated accidents or a systematic feature of the landscape, and would guide the development of initialization protocols that avoid them without requiring trial-and-error. The paper's Figure S2 (gradient norm increasing transiently near the saddle point) suggests a possible diagnostic: monitor the gradient norm and reinitialize if it spikes after a long plateau. Validating this diagnostic across multiple systems would provide a practical tool.

Training a classifier to predict difficulty of state separation directly from molecular structure. The paper's two most challenging systems—fluorine/neon (many closely-spaced Rydberg states) and cyclopentadienone/tetrazine (dense manifolds of triplet states requiring singlet-targeting)—share a common feature: high density of states near the excitation energies of interest. The paper handles this with ad-hoc solutions (increase $K$ for atoms, add $\lambda \hat{S}^2$ for the double excitations), but a more principled approach would predict a priori how many states $K$ are needed to resolve the spectrum of interest. A specific experiment: for a range of molecules with varying state densities (from atoms through polycyclic aromatics), compute the full spectrum at some level of theory (e.g., TD-DFT or EOM-CCSD), measure the spacing between the $K$ -th and $(K+1)$ -th states, and correlate this with features of the molecular structure (number of electrons, density of virtual orbitals near the HOMO-LUMO gap, presence of heavy atoms). Train a lightweight model to predict the required $K$ from these features. Validate by running NES-VMC with the predicted $K$ and checking whether states of interest are resolved without contamination from higher states.

Combining NES-VMC with geometry optimization for excited-state potential energy surfaces. The paper demonstrates single-point vertical excitation energies and a manually-sampled potential energy curve for carbon dimer and twisted ethylene, but does not perform automated geometry optimization on excited-state surfaces. This is the natural next step for studying photochemical reaction mechanisms: locate excited-state minima, transition states, and conical intersections automatically. The challenge is that NES-VMC provides energies and gradients for $K$ states simultaneously, but the energy of a specific state is only recovered after diagonalizing the local energy matrix, which introduces non-differentiability into the state-assignment step. A specific development: implement a state-tracking algorithm that follows the eigenvectors of the local energy matrix across geometry steps (by maximizing overlap with the previous step's eigenvectors), enabling smooth potential energy surfaces for individual states without manual state assignment. Test on the twisted ethylene conical intersection (where the paper already has reference geometries from Barbatti et al., 2004) by running NES-VMC at each step of an automated minimum-energy conical intersection search and comparing the located intersection geometry to the MR-CI reference.

Stress-testing on transition metal complexes. The paper studies exclusively main-group atoms and molecules (first-row atoms, small organic molecules, benzene). Transition metals introduce qualitatively new challenges: near-degenerate $d$ -orbital manifolds, multiple low-lying spin states (singlet, triplet, quintet often within ~1 eV), strong static correlation requiring many determinants, and relativistic effects that become significant for heavier elements. A specific stress-test: compute the lowest $K$ states of an octahedral $d^6$ complex (e.g., $[\text{Fe}(\text{H}_2\text{O})_6]^{2+}$ ) with NES-VMC using the Psiformer and pseudopotentials for iron. The known spin-state ordering (high-spin $^5T_{2g}$ ground state, low-spin $^1A_{1g}$ at higher energy, with several intermediate-spin states) provides a challenging test of whether the determinant construction can prevent collapse across spin multiplicities. The paper's singlet-targeting trick ( $\lambda \hat{S}^2$ ) would need to be applied separately for each spin multiplicity, or the method would need to compute a very large $K$ to resolve states of all spins simultaneously—a direct test of the $O(K^3)$ scaling barrier documented in Figure S3.

Developing a "natural gradient" optimizer specifically for NES-VMC. The paper uses standard stochastic gradient descent (with the gradient scaling from the Psiformer paper, learning rate 0.05) and notes saddle-point convergence as a challenge. Natural gradient methods (which precondition the gradient by the inverse Fisher information matrix) are known to escape saddle points more effectively in some settings because they account for the geometry of the parameter space. The NES-VMC objective $\mathbb{E}_{x \sim \Psi^2}[\text{Tr}[E_L(x)]]$ has a specific structure: it is the trace of a matrix expectation. The Fisher information matrix for the distribution $\Psi^2$ is computable from the same samples used for the energy gradient, since it involves expectations of outer products of the score function $\nabla_\theta \log|\Psi|$ . A concrete implementation: adapt the Kronecker-Factored Approximate Curvature (KFAC) optimizer (already used in the FermiNet and Psiformer implementations) to the multi-state setting by computing the Fisher information from the total ansatz $\Psi$ rather than from individual $\psi_i$ , and testing whether this reduces the frequency of saddle-point trapping on the systems where the paper observed it (especially nitrosomethane, butadiene, and glyoxal with ordered pretraining).

Combining NES-VMC with diffusion Monte Carlo for excited states. The paper notes that Pathak et al. (2021) used a DMC correction to reduce VMC penalty method errors by ~0.1 eV on benzene. NES-VMC provides a better VMC starting point (Psiformer error ~0.09 eV on benzene $^3B_{1u}$ vs. the penalty method's ~0.29 eV corrected), so a DMC correction on top of NES-VMC might push accuracy below 0.05 eV—competitive with the best coupled-cluster results. The specific challenge is that DMC for excited states typically requires fixed-node approximations where the nodal surface is taken from the VMC wavefunction, and NES-VMC provides $K$ wavefunctions that are not individually orthogonal. A specific experiment: take the NES-VMC-optimized single-state ansätze for a small molecule (e.g., H $_2$ O or C $_2$ H $_4$ ), diagonalize the local energy matrix to recover approximate eigenstates, use each eigenstate's nodal surface as a fixed-node constraint for a separate DMC calculation, and measure the reduction in excitation energy error relative to the VMC result. Compare against the DMC correction magnitude reported in Pathak et al. If the DMC correction is smaller for NES-VMC than for penalty-method VMC (because NES-VMC's VMC starting point is already more accurate), that would suggest NES-VMC is recovering more of the dynamic correlation at the VMC level, leaving less for DMC to correct—an important characterization of the method's completeness.

Practical Applications and Downstream Use Cases

Screening candidate molecules for photocatalytic activity by predicting absorption spectra. The paper demonstrates that NES-VMC can compute oscillator strengths to within ~0.001 of theoretical best estimates for a diverse set of small molecules (Section IV, Figure 2, Table S3), and can handle double excitations that single-reference methods miss entirely (Section VII). For photocatalytic applications—where the key quantity is whether a molecule absorbs light at a specific wavelength with sufficient oscillator strength to drive a reaction—TD-DFT is the standard screening tool but "often gives qualitatively incorrect results" (Section I). NES-VMC provides an alternative that is more accurate than TD-DFT (as demonstrated on benzene, Figure 6, where TD-DFT-PBE0 deviates significantly from the TBE while NES-VMC closely matches) and can be applied without prior knowledge of the excited-state character (single vs. double excitation). The practical workflow: for a candidate photocatalyst with 20–40 electrons (within the scale where the paper demonstrates robust performance), run NES-VMC with $K = 5$ –10 states and random pretraining, compute the first few singlet excitation energies and oscillator strengths, and compare against the target wavelength. The paper's results on butadiene—where NES-VMC correctly resolves the historically controversial $1^1B_u$ / $2^1A_g$ ordering that TD-DFT with standard functionals often gets wrong—directly demonstrate the value for conjugated organic molecules relevant to photovoltaics and photosynthesis. The cost is higher than TD-DFT (hours to days on a GPU vs. seconds on a CPU), but the accuracy gain—avoiding qualitatively incorrect state orderings—may justify the cost for lead candidates before experimental synthesis.

Generating reference data for training machine learning force fields that include excited states. The paper's ability to compute off-diagonal observables (transition dipole moments) and potential energy surfaces (carbon dimer, Figure 3; twisted ethylene, Figure 4) makes it suitable for generating training data for neural network potential energy surfaces that must describe multiple electronic states and non-adiabatic couplings. Current approaches for training such potentials rely on MR-CI or CASPT2 calculations, which become prohibitively expensive for systems with more than ~10–15 heavy atoms. NES-VMC scales as $O(N^3)–O(N^4)$ with system size—favorable for a wavefunction method—and provides explicit wavefunction amplitudes that can be used to compute non-adiabatic coupling vectors (derivatives of one state's wavefunction with respect to nuclear coordinates, projected onto another state). The paper does not compute non-adiabatic couplings, but the framework for off-diagonal observables (Section S2, Eq. 51) extends naturally: the coupling vector involves the matrix element $\langle \psi_i | \nabla_R \psi_j \rangle$ , which can be accumulated as an off-diagonal observable matrix and demixed using the same eigenvectors $U$ as the energy and dipole moment. For a specific system like the twisted ethylene conical intersection (where the paper already has validated geometries and energies), adding non-adiabatic coupling calculations would provide the complete dataset needed for a machine learning model of the $N$ / $V$ state crossing—a standard benchmark for non-adiabatic dynamics methods.

State characterization via natural orbital analysis for understanding photochemical mechanisms. The paper demonstrates on the carbon dimer (Section V, Figures 3d–3e) and benzene (Section VIII, Figure 6b) that the optimized wavefunctions can be analyzed to extract natural orbitals and their occupancies, revealing the character of each excitation (single vs. double, which orbitals are involved). This capability is absent from projector QMC methods (DMC, AFQMC), which provide only energies, and from TD-DFT, which provides only approximate transition densities. For a researcher trying to understand why a particular molecule photoisomerizes or where a conical intersection occurs, knowing that a state has double excitation character (as NES-VMC reveals for the $b^3\Sigma_g^-$ and $B^1\Delta_g$ states of C $_2$ ) directly explains why single-reference methods fail to describe it correctly. The computational cost of this analysis is modest: the density matrix accumulation (Section S13B) requires additional sampling from the same $\Psi^2$ distribution, and the natural orbital extraction involves diagonalizing matrices of size $M \times M$ (where $M$ is the number of basis functions used for projection, typically 50–200 for the def2-TZVPD basis used in the paper). This enables a workflow where NES-VMC serves not just as an energy calculator but as a diagnostic tool for understanding the electronic structure of excited states—complementary to the energetic benchmarks that dominate the paper but potentially more valuable for mechanistic insight.

When to Prefer This Method

The paper positions NES-VMC explicitly against penalty-based VMC methods for excited states (particularly Entwistle et al., 2023 and Wheeler et al., 2024) and implicitly against conventional quantum chemistry methods (TD-DFT, CASPT2, EOM-CCSD) by benchmarking against their results. The following decision criteria are grounded in the paper's specific findings:

Prefer NES-VMC with neural network ansätze when:

The system has ground-state VMC results that are already accurate (indicating the ansatz has sufficient expressivity), but excited-state VMC with penalty methods fails—as on benzene, where the PauliNet penalty method produces 1.4 eV error while NES-VMC with the Psiformer achieves 0.088 eV on the $^3B_{1u}$ state (Section VIII). The determinant construction eliminates the penalty strength hyperparameter entirely.
The excited states of interest include double excitations or states with significant multi-reference character, where TD-DFT with standard functionals gives qualitatively wrong results and even CASPT2 may have large errors. NES-VMC correctly identifies double excitation character in C $_2$ (Section V, Figure 3e) and butadiene (Section VII, where the $2^1A_g$ state has ~30% double excitation character), and resolves QUEST database errors as large as ~1 eV on tetrazine.
Physical insight from the wavefunction (natural orbitals, density matrices, excitation character) is as valuable as the energies themselves. NES-VMC provides explicit wavefunction forms enabling the analyses in Figures 3d–3e and 6b that projector methods cannot replicate. If you need to understand why an excitation occurs (which orbitals are involved, whether it is single or double), NES-VMC is preferable to DMC or AFQMC.
Off-diagonal properties like oscillator strengths are required. NES-VMC computes transition dipole moments from the same accumulated observable matrices as the energies, achieving oscillator strengths within ~0.001 of TBEs (Section IV, Table S3). The paper notes that oscillator strengths are "known to be highly sensitive to the choice of basis set," yet NES-VMC (which is basis-set-free) matches basis-set-extrapolated coupled-cluster results—suggesting the method captures the relevant physics directly.

Prefer conventional quantum chemistry (TD-DFT, EOM-CCSD, CASPT2) when:

The molecule is larger than ~40 electrons or contains many heavy atoms (beyond second row). The paper's largest system is benzene (42 electrons), the bottleneck envelope was required to manage memory, and scaling beyond this point is uncharacterized. EOM-CCSD with a triple-zeta basis can handle 50–100 electrons routinely on CPU clusters; TD-DFT can handle hundreds of atoms. NES-VMC's $O(N^3)–O(N^4)$ scaling with GPU implementation is competitive in principle but unproven above benzene scale.
The system has a very high density of low-lying excited states of different spin multiplicities, requiring the singlet-targeting trick (Section S5) to separate states. The paper used this for cyclopentadienone and tetrazine because "a large number of low-lying triplet states meant that it was impractical to compute enough states." This workaround requires choosing $\lambda$ (the spin shift) and means singlet and triplet states cannot be computed simultaneously. For transition metal complexes where spin-state energetics are the central question, this is a significant limitation.
The geometry or state assignment is known to be sensitive to pretraining, and you cannot afford multiple runs with different initializations to verify convergence. The paper shows that ordered and random pretraining find different states on HNO and HCF (Figure S1), and that on double excitation systems the best result was selected post-hoc by trying both. Without a reliable diagnostic for saddle-point trapping, a single NES-VMC calculation cannot be trusted to have found the correct states unless the expected spectrum is already known—in which case the calculation is confirming rather than predicting.
Absolute energies (not just excitation energies) are needed at high precision. The paper notes that "optimal absolute energies" require 200,000 iterations versus 100,000 for converged relative energies, and that pseudopotentials were needed for heavy atoms because "the energy scale of the core electrons is far greater than the valence electrons." For applications like thermochemistry where total energies must be accurate to ~1 mHa, the convergence requirements and core electron treatment may make NES-VMC less practical than explicitly correlated coupled-cluster methods that handle core-valence correlation systematically.

Prefer projector QMC (DMC, AFQMC) when:

Only the lowest one or two excited states are needed at very high accuracy, and the fixed-node approximation from a cheaper VMC method is acceptable. Pathak et al. (2021) achieved DMC-corrected energies on benzene with errors comparable to NES-VMC-VMC (Section VIII, Table S10), and DMC can be applied after a sufficiently accurate VMC calculation regardless of which VMC method produced the trial function. If NES-VMC provides a better VMC starting point than penalty methods (as the paper argues), then NES-VMC + DMC may be a powerful combination, but this is not demonstrated in the present paper.