Insights from exact exchange-correlation kernels

The iDEA code [44] is used in order to obtain the interacting and Kohn-Sham response functions. This software implements quantum mechanics for finite systems in one dimension interacting with a softened Coulomb electron-electron interaction

where $\alpha$ is the extent of the softening; we use $\alpha=1$ a.u. A delta function basis set, i.e. real-space grid, of dimension $N$ is used to discretize the spatial domain $[-a,a]$ of length $L=2a$ subject to Dirichlet boundary conditions.

Our three prototype systems each consist of two spinless electrons in the external potentials described below. (The use of spinless fermions delivers a richer ground state and excitation spectrum for a given number of electrons.)

For some input external potential $v_{\text{ext}}(x)$, the full set of eigenvectors $\{|\Psi_{i}\rangle\}$ of the interacting Hamiltonian is found using exact diagonalization. The corresponding exact Kohn-Sham potential $v_{\text{KS}}(x)$ is then reverse-engineered by applying preconditioned root-finding techniques to an appropriate fixed-point map [45]. The full set of Kohn-Sham eigenvectors $\{|\phi_{i}\rangle\}$ is also obtained using exact diagonalization. The causal response functions are computed in the frequency domain directly using the Lehmann representation; the interacting response function, for example, is given by

where $\Omega_{n}=E_{n}-E_{0}$ is the $n^{\text{th}}$ excitation energy of the interacting Hamiltonian, and the response function is zero for times $t>t^{\prime}$. Construction of the interacting response function in this fashion is an accurate but demanding procedure, whereas the methods outlined in [42] to construct the response functions are amenable to larger systems, but more prone to error; either method will suffice here. On a finite spatial grid, the response functions at a given $\omega$ become response matrices, denoted $\chi(\omega)$ and $\chi_{0}(\omega)$. The Dyson equation gives an alternate definition of the xc kernel,

where ^-1 is to be understood as the matrix inverse, and the Hartree kernel $f_{\text{H}}$ becomes softened due to Eq. (4). This expression is used as the definition of $f_{\text{xc}}$ in the present context, where the matrix inverses require the careful treatment described in the next section.

Numerical difficulties arise when attempting to construct $f_{\text{xc}}$ as an object in itself using Eq. (5), which is one of a few reasons that has prevented or hindered studies of the exact $f_{\text{xc}}$ along these lines [41, 42, 43]. The xc kernel represents the solution to an inverse problem, i.e. find the $\delta v_{\text{xc}}$ that produces a given $\delta n$, and inverse problems are notoriously sensitive to small error [46], such as those introduced by finite-precision arithmetic. As discussed, $f_{\text{xc}}$ in Eq. (5) requires the matrix inverse of the response matrices at a given $\omega$, and hence a naive inversion procedure introduces numerical error at a given $\omega$ in proportion to the condition of the response matrices, that is, the ratio of the maximum to minimum eigenvalue.

Physical eigenvalues that are close to, or below, machine precision manifest in the response matrices from various sources [47]. One such source is related to the linear response $v$-representability problem. That is to say, there exist density perturbations that oscillate with some non-resonant frequency $\omega$ that cannot be produced by a perturbing potential at linear order [48, 47, 27]. Therefore, at such a frequency, the response function has an eigenvector $|u(\omega)\rangle$ whose eigenvalue is zero, indicating that the density perturbation $\delta n=|u(\omega)\rangle$ does not correspond to some finite $\delta v$, as the response function is non-invertible ²²2For example, the constant perturbation $\delta v=c(\omega)$ oscillating with frequency $\omega$ produces no response in the density $\delta n$ at all orders, including at first order.. This can happen in both the interacting and Kohn-Sham response functions at distinct frequencies. These eigenvalues are of importance for the work to follow, and are discussed in more depth in Section III.2.

Another source of low eigenvalues is due to extended regions of nearly vanishing ground-state density. Since the aim of this work is, in part, to study the optical response of confined systems far beyond the first excitation, an appropriately large spatial domain $[-a,a]$ is required in order to accommodate the more extended excited states without introducing spurious features due to the boundary conditions. Within such a system, a perturbing potential localized toward the edge of the domain yields a negligible ground-state density response, the effect of which is to introduce near-zero eigenvalues into the response functions. Therefore, ill-conditioning is unavoidable if we are to study the response functions, and thus $f_{\text{xc}}$, at high frequencies. The extent of this ill-conditioning depends on the maximum frequency up to which one wishes to examine matters.

Note that these near-machine-precision eigenvalues of the response functions are problematic only if this difference accounts for some particular physical phenomenon, such as charge transfer, whereby in order to capture excitations of charge-transfer character an $f_{\text{xc}}$ that is divergent in proportion to the increasing separation between the subsystems involved is required ³³3We observed this situation in a double-well system that we explored as background to the present study. [51]. In such cases the construction of numerical xc kernels will be challenging. However, the systems studied in this work do not suffer this issue: the ill-conditioning of the interacting and Kohn-Sham response functions can be assumed to cancel in the definition of $f_{\text{xc}}$, Eq. (5), thus producing a regular $f_{\text{xc}}$. This procedure can be viewed as a form of basis set truncation, i.e. assign $\chi=\chi_{0}$ within some subset of the basis responsible for ill-conditioning and proceed to compute $f_{\text{xc}}$ under this assumption. We now describe two such approaches: a truncation in real space, and a truncation in eigenspace.

As noted in Refs. [28, 27, 57, 37], the following transformation

leaves the output of the Dyson equation unchanged, and thus we are, in principle, free to choose the arbitrary complex-valued functions $g(x,\omega),h(x^{\prime},\omega)$ and $c(\omega)$. All three transformations are a direct manifestation of the invariance of quantum Hamiltonian systems under a constant time-dependent shift of the potential. From the point of view of $f_{\text{xc}}$ approximations, two xc kernels are equivalent if they exist within this family of functions ⁶⁶6Note that the quantities $\langle ij|f_{\text{xc}}(\omega)|kl\rangle$, where $(i,j,k,l)$ label indices of single-particle wavefunctions, are unique [27], and it is these quantities that form the input to the Casida equation [72], for example..

A preferred gauge is defined by Eq. (5), since the objects $\chi$ and $\chi_{0}$ are themselves invariant to a shift in the potential. The unique $f_{\text{xc}}$, modulo a constant shift (see below), defined in Eq. (5) can be considered the physical $f_{\text{xc}}$, and it is this definition of $f_{\text{xc}}$ that is assumed in discussions on its various properties and limits [16, 15, 52, 53, 59]. To modify this $f_{\text{xc}}$ using its gauge freedom changes its underlying structure; for example, setting $g\neq 0$ gives $f_{\text{xc}}$ spurious long-range behavior, and setting $g\neq h$ produces an $f_{\text{xc}}$ that is not symmetric under interchange of $x\leftrightarrow x^{\prime}$. In this work, we illustrate $f_{\text{xc}}$ as it is defined in Eq. (5), which also defines $g=h=0$. The constant shift $c$ has special meaning, as it is itself an eigenvector of $\chi$ and $\chi_{0}$ with eigenvalue zero, i.e. the response functions are non-invertible in this direction. Since the Dyson equation is therefore silent regarding the value of $c$, we anchor $f_{\text{xc}}$ by requiring that, in the long-range limit, far outside the confined density, $f_{\text{xc}}+f_{\text{H}}\rightarrow 0$, and find that this limit is reliably achieved.

Having decided upon a preferred gauge, one can consider the possible consequences of this gauge freedom on matters of practical interest. An approximate $f_{\text{xc}}$ that differs in relevant structure from the exact $f_{\text{xc}}$ largely due to a change of gauge provides at least a partial explanation for the performance of a given approximate $f_{\text{xc}}$; in Section III.4 we consider this line of inquiry.

Our first system consists of two interacting electrons confined in the atom-like potential, $v_{\text{ext}}(x)=-2/(|0.1x|+0.5)$, within the domain $[-8,8]$ a.u. An illustration of this system, and its associated numerical parameters, are given in the supplemental material. The purpose of the atom demonstration is two-fold. First, it constitutes a proof-of-concept, and defines a standard of accuracy to which the remainder of the calculations are held unless stated otherwise. Second, the optical spectrum of the atom is calculated in the range $\omega\in[0,6]$ a.u., which includes many excitations beyond the first, and the efficacy of various approximations to $f_{\text{xc}}$ are examined in relation to the optical spectrum.

The exact $f_{\text{xc}}$, constructed using the real-space truncation method, is shown for the first three visible interacting excitations in the optical spectrum, and at $\omega=0$ a.u., in Fig. 3. The last is sometimes termed the exact adiabatic $f_{\text{xc}}$, and it correctly describes any system in which the response to a perturbation is essentially instantaneous, $f_{\text{xc}}(x,x^{\prime},\omega=0)$.

The real-space truncation parameter is chosen as $b=5.8$ a.u., meaning the inner region is defined as $-5.8<x<5.8$. It is important to stress that this choice of $b$ is not unique, and there exists some feasible range of $b$ within which $f_{\text{xc}}$ itself is insensitive to changes. Moreover, within this feasible range, both the method error and numerical error are acceptable – a discussion on the precise quantification of error here is given in the supplemental material. The mean absolute error between the output of the Dyson equation, which is hereafter defined as

and the interacting response function $\chi(\omega)$ is $\mathcal{O}(10^{-9})$ over the entire grid. The zero-force sum rule [2, 60] is used to further validate the numerics, which is discussed and illustrated in the supplemental material.

In order to extract the optical transition energies and transition rates, given a single-particle response function $\chi_{0}$ and xc kernel $f_{\text{xc}}$, one can construct the entire optical absorption spectrum Eq. (3) using the corresponding output of the Dyson equation, denoted $\chi_{\text{Dyson}}(f_{\text{xc}},v_{\text{xc}})$, where $\chi_{0}$ is specified with some $v_{\text{xc}}$, see Fig. 4.

As established in [61, 62, 63], the exact Kohn-Sham single-particle transitions are in excellent agreement with the interacting transitions, but this agreement becomes increasingly poor at higher energies. An approach to understanding this is to consider the overlap between the final states involved in a given interacting $|\Psi_{0}\rangle\rightarrow|\Psi_{f}\rangle$ and Kohn-Sham $|\Phi_{(0,1)}\rangle\rightarrow|\Phi_{f}\rangle$ transition, where $|\Phi_{(i,j)}\rangle$ denotes the Slater determinant constructed from $i^{\text{th}}$ and $j^{\text{th}}$ Kohn-Sham single-particle states.

The overlap of the ground-state is $\langle\Psi_{0}|\Phi_{(0,1)}\rangle=0.99991$, and the overlap of the final states involved in the first transition at $\omega=0.76$ a.u. is $\langle\Psi_{1}|\Phi_{(0,2)}\rangle=0.9995$. The static correlation in the interacting state here is modest, meaning the interacting state has strong single-particle character, which leads to agreement between the low-energy transitions in the optical spectrum. At higher energies, the overlap decays by multiple orders of magnitude, however this is not the predominant source of error in higher energy transitions. Rather, interacting excitations that correspond to Kohn-Sham single-particle excitations out of the highest occupied state (the second) are much more accurate than single-particle excitations out of the first state. For example, the interacting excitation $|\Psi_{0}\rangle\rightarrow|\Psi_{19}\rangle$ at $\omega=4.41$ a.u. in Fig. 4 is captured well with the Kohn-Sham excitation $|\Phi_{(0,1)}\rangle\rightarrow|\Phi_{(0,12)}\rangle$, whereas this is not true of the preceding interacting excitation. This is not surprising, as the highest occupied Kohn-Sham state has energy equal to minus the exact electron removal energy [64], and thus at least one energy involved in the transition is correct. This might often be the case, and to comment further would require additional many-body calculations.

The following $f_{\text{xc}}$ approximations are now considered: the RPA $f^{\text{RPA}}_{\text{xc}}=0$, an adiabatic LDA $f^{\text{ALDA}}_{\text{xc}}[n](x,x^{\prime},\omega=0)\propto\delta(x-x^{\prime})$ parameterized with reference to the homogeneous electron gas in [65], and the exact adiabatic xc kernel $f_{\text{xc}}(x,x^{\prime},\omega=0)$, Fig. 3(a). These $f_{\text{xc}}$ approximations are used to solve the Dyson equation in conjunction with the exact Kohn-Sham response function. The atomic optical spectrum, using the aforementioned series of approximations, is shown in Fig. 5 for the first transition and a chosen higher energy transition.

The exact adiabatic $f_{\text{xc}}$ reproduces the first peak well, which reflects the fact that $f_{\text{xc}}$ at the first excitation, Fig. 3(b), is visually indistinguishable from the exact adiabatic $f_{\text{xc}}$. This agreement demonstrates that, not only is the adiabatic approximation valid for low-energy transitions, but also the non-local spatial structure in the exact adiabatic $f_{\text{xc}}$ is required in order to reproduce the low-energy peaks in the optical spectrum – in lacking this structure, both the RPA and adiabatic LDA significantly over-correct the underestimation of the non-interacting transition energy. At higher energies, i.e. beyond the third peak in the optical spectrum (not shown), all three approximations perform similarly, and do not improve matters significantly beyond the corresponding non-interacting peak. This behavior appears to be generic for all peaks observed up to $\omega=6$ a.u., namely, the transition energies output from the Dyson equation with these approximate xc kernels are bound to the quality of the non-interacting transition energies. Furthermore, the inset of Fig. 5 demonstrates that the exact adiabatic approximation gives an excitation energy that is worse than the other adiabatic approximations. This suggests that, in order to capture higher energy transitions, a frequency dependence is required, and in particular ‘improving’ the spatial structure of adiabatic approximations toward the exact adiabatic structure does not assist matters here.

As alluded to above, the exact adiabatic approximation ceases to out-perform the adiabatic LDA and RPA beyond the third transition in the optical spectrum, i.e. the same transition for which the corresponding exact $f_{\text{xc}}$ departs from its adiabatic structure in a serious manner, Fig. 3. In fact, as the subsection to follow demonstrates, this violent departure from adiabaticity has a particular and fairly simple form, and its origin is understood in the context of eigenvalues of $\chi(\omega)$ or $\chi_{0}(\omega)$ that cross zero. The eigenvector corresponding to the eigenvalue that touches zero is a non-$v$-representable density perturbation at linear order.

The infinite potential well is defined with the external potential $v_{\text{ext}}=0$ inside the domain $[-5,5]$ a.u. (See supplemental material for the corresponding Kohn-Sham potential, density, and numerical parameters.) The numerical response functions for this system are well-conditioned and valid up to arbitrary $\omega$, there are no regions of nearly vanishing density, and thus there is no significant numerical error present here.

The non-adiabatic character of $f_{\text{xc}}$ is illustrated up to $\omega\approx 6$ a.u. in Fig. 6. As in the atom, $f_{\text{xc}}$ exhibits little frequency dependence, until after the second transition. This behavior is related to the poles that occur in $f_{\text{xc}}$ infinitesimally below the $\omega$-axis [47, 27]. The eigenvalues of $f_{\text{xc}}$ as a function of frequency are given in Fig. 7, and in particular, three divergences are shown (the singularities are tempered slightly by evaluating the response functions just above the real $\omega$-axis). The lower panels of Fig. 7 demonstrate that these singularities coincide with a single eigenvalue of either the interacting or non-interacting response function crossing zero. Thus, the visual character of $f_{\text{xc}}$ is dominated by the outer product of the eigenvector whose eigenvalue is either beginning to diverge, or recovering from a divergence. Moreover, nothing in principle is preventing these singularities in $f_{\text{xc}}$ coming arbitrarily close to an interacting excitation – either $\chi_{0}$ can cross zero close to an interacting excitation, or $\chi$ itself can cross zero close to an interacting excitation ⁷⁷7The interacting response function $\chi$ diverges at an interacting excitation, but this does not, in a finite basis, prevent an eigenvalue from crossing zero at this energy under certain circumstances that are set out in the supplemental material.. Since most $f_{\text{xc}}$ approximations lack these divergences, one can question the importance of these divergences in relation to optical properties.

In order to determine the impact of the divergences in $f_{\text{xc}}$ on the optical spectrum, we examine how the optical spectrum is affected after projecting out the divergent eigendirection in $f_{\text{xc}}$, but otherwise keeping $f_{\text{xc}}$ identical (the details of this procedure are given in the supplemental material). This is tantamount to setting $\chi|v\rangle\coloneqq\chi_{0}|v\rangle$ for some vector $|v\rangle$ within the Hilbert space that is associated with the divergence.

The interacting excitation at $\omega=1.09$ a.u., see Fig. 7, is visible in the optical spectrum, and furthermore the character of $f_{\text{xc}}$ at this energy, see Fig. 6, is dominated by the outer product of an eigenvector whose eigenvalue is much larger in magnitude than the rest, and between the two divergences at $\omega=1.01,1.17$ a.u. The removal of these divergences from $f_{\text{xc}}$ across the energy range of interest, and subsequently solving the Dyson equation with the projected xc kernel $f_{\text{xc}}^{\text{projected}}$, shifts the interacting optical peak back toward the non-interacting peak, Fig. 8. Conversely, the much weaker divergence at $\omega=0.99$ a.u., as seen in Fig. 7, has a tail that also yields an eigenvalue much larger in magnitude than the rest at the interacting excitation, and removal of this divergence produces no change in the optical spectrum, i.e. this eigenvector is not relevant for capturing the transition in question. In the inset of Fig. 8, the visual character of $f_{\text{xc}}$ is shown at the interacting excitation $\omega=1.09$ a.u. after the divergences visible in Fig. 7 have been removed. Underneath these divergences, $f_{\text{xc}}$ is indistinguishable from the adiabatic $f_{\text{xc}}$, meaning the frequency dependence of $f_{\text{xc}}$ in this system is largely due to its pole structure. These results suggest that functional approximations that do not capture the non-adiabatic pole structure of $f_{\text{xc}}$ will struggle to improve matters beyond the non-interacting peaks for certain transitions.

A further point of note is that there are many more zeros in the interacting response function than the non-interacting response function, as a simple counting argument is sufficient to demonstrate. All $N$ eigenvalues of the response functions begin negative [47], and each excitation brings a negative eigenvalue to a positive eigenvalue across a divergence, which otherwise evolves as a continuous function of $\omega$. For two electrons discretized with a basis set of dimension $N$, there are $\frac{1}{2}N(N-1)-1$ interacting excitations, and $2N-4$ non-interacting excitations, the difference being made up of double (triple, etc. with more than two electrons) excitations that are notoriously not present in the Kohn-Sham response function [21]. Therefore, there are a great deal more eigenvalues that must pass through zero in the interacting response function, and moreover, these eigenvalues that cross zero are connected to the excitations that require them to do so – an account of the precise conditions under which this occurs is given in the supplemental material using a two-state model, see also [65].

In Fig. 7, there are three divergences, two of which are paired at $\omega=1.01,1.17$ a.u., meaning the zero in the non-interacting response function is a shifted version of the zero in the interacting response function; such behavior can be related to single excitations. There also exists an unpaired divergence at $\omega=0.99$ a.u. related to the excitation at $\omega=0.98$ a.u. which has double excitation character. That is, the overlap between the final state involved in this transition $|\Psi_{0}\rangle\rightarrow|\Psi_{5}\rangle$ and the Slater determinant $|\Phi_{(2,3)}\rangle$ is $\langle\Phi_{(2,3)}|\Psi_{5}\rangle=0.98$.

It is perhaps, then, no surprise that removal of the eigenvector and eigenvalue whose source is the unpaired divergence did not alter the visible (single excitation) transition in Fig. 8. In fact, the pole in the interacting response function relating to the double excitation disappears with removal of the unpaired divergence, and so this divergence, and its surrounding character, is important in order to capture the transition for which it is relevant. The xc kernel exhibits unpaired divergences after all double excitations up to $\omega=6$ a.u. at frequencies slightly higher than the interacting double excitation energy. Refs. [21, 23] derive, under certain assumptions, the necessarily divergent character of $f_{\text{xc}}$ around a double excitation. The above analysis reveals that divergences in $f_{\text{xc}}$ are, in fact, common and associated with the $\omega$ neighborhood containing multiple and single excitations alike.

The quantum harmonic oscillator is defined in the domain $[-8,8]$ a.u. with the potential $v_{\text{ext}}(x)=\frac{1}{2}\nu^{2}x^{2}$ where $\nu\coloneqq 0.45$ a.u. (See supplemental material for the corresponding exact Kohn-Sham potential, the interacting ground-state density, and the numerical parameters.)

First, the spatial structure, and in particular the long-range behavior [16], of the exact $f_{\text{xc}}$ is examined. The delicate nature of the numerics involved in computing the exact $f_{\text{xc}}$ is brought to the fore when attempting to capture its long-range limit. The ill-conditioning in the response matrices can be identified with regions of nearly vanishing ground-state density, see Section II.2, and it is precisely in these regions that we expect to observe the long-range character of $f_{\text{xc}}$. However, $f_{\text{xc}}$, when evaluated outside some central region where we can be confident there is little-to-no numerical error, diverges in a manner that is not consistent with the known long-range limit of $f_{\text{xc}}$. This spurious divergence in $f_{\text{xc}}$ does not much affect the accuracy of the output of the Dyson equation, $\chi_{\text{Dyson}}$, because the Dyson equation Eq. (10) involves the matrix product $\chi_{0}f_{\text{xc}}$, and the divergent regions of $f_{\text{xc}}$ operate on the nearly vanishing regions of $\chi_{0}$.

If we are to observe the long-range limit of $f_{\text{xc}}$ in the present context, the region where the numerical error is low must overlap with the region where the long-range limit is observed. This is the case for the exact adiabatic $f_{\text{xc}}$, which is shown, together with slices of $f_{\text{xc}}$ along a particular axis, in Fig. 9 and Fig. 10.

The center-most region of the domain is where exchange and correlation effects are most important, and in this region $f_{\text{xc}}$ exhibits a strong local response that quickly decays as $x\neq x^{\prime}$. Since this region can be interpreted as the most crucial for recovering accurate observable properties from the Dyson equation, this observation supports, at least in part, local approximations to $f_{\text{xc}}$, such as the adiabatic LDA.

The non-local structure of $f_{\text{xc}}$ at $\omega=0$ a.u. can be seen most clearly in the lower panel of Fig. 10, where perturbations in the density outside the center-most region cause a significant change in the exchange-correlation potential inside the center-most region. The failure of the adiabatic LDA to capture this non-local response leads to fairly poor agreement between the exact and approximate transition energies, which is seen for the atom in Fig. 5, and illustrated for the quantum harmonic oscillator below ⁸⁸8It is possible that the particular form of the non-local $f_{\text{xc}}$ functionals considered in [73, 13], in which $f_{\text{xc}}(x,x^{\prime})=f_{\text{xc}}[(n(x)+n(x^{\prime}))/2]$, can capture certain features of the adiabatic non-local response depicted in Fig. 9 and Fig. 10..

The long-range limit $f_{\text{xc}}(x,x^{\prime},\omega=0)\rightarrow-f_{\text{H}}(|x-x^{\prime}|)$ is explored in Fig. 10 along both $x^{\prime}=0$ a.u. and $x^{\prime}=-2$ a.u. Due to the rapid decay of the ground-state density in the quantum harmonic oscillator, convergence of $f_{\text{xc}}$ toward $-f_{\text{H}}$ along $x=-2$ a.u. is not directly observed in the $x\rightarrow-\infty$ limit – the atom, whose ground-state density does not decay as quickly, converges toward $-f_{\text{H}}$ in both the positive and negative limits (see supplemental material). It is known that the long-range limit of $f_{\text{xc}}$, in both finite and periodic systems, satisfies $f_{\text{xc}}(\omega)\rightarrow-\alpha(\omega)f_{\text{H}}$ [16, 15, 52, 53], where $\alpha(\omega)$ is a frequency-dependent constant that reflects dielectric screening in the system. Therefore, one expects $\alpha=1$ in finite systems, and at lower frequencies, where the numerical methodology provides a robust long-range limit, our observations are consistent with this, see Fig. 10.

To conclude matters for the quantum harmonic oscillator, we examine its optical absorption spectrum, and in particular highlight a failing of the $f_{\text{xc}}$ approximations considered in this work when compared to the exact and exact adiabatic xc kernels. The optical absorption spectrum, using the same range of approximations discussed in the case of the atom, is shown in Fig. 11. In the non-interacting quantum harmonic oscillator, all transitions but the first are disallowed, otherwise known as its selection rules. The inclusion of the Coulomb interaction lifts these special symmetries of the non-interacting quantum harmonic oscillator, but not enough for the previously disallowed transitions to be observed in the optical spectrum – the dipole matrix elements for these allowed transitions are $\mathcal{O}(10^{-8})$. On the other hand, the exact Kohn-Sham potential differs significantly from the harmonic form, which creates a series of visible peaks beyond the first in the optical spectrum; the transition rates for these transitions are vastly overestimated.

In this situation, the RPA and adiabatic LDA do not achieve the required strong suppression of the optical peaks. Interestingly, the exact adiabatic kernel, as seen in Fig. 9, is able to reproduce the exact optical spectrum across the entire frequency range considered $\omega\in[0,6]$ a.u., presumably due to the correct oscillator strengths used in its construction. This suggests that, in cases where the exact system possesses heavily suppressed transitions, perhaps due to symmetries that are not shared by the Kohn-Sham system, the typical $f_{\text{xc}}$ approximations are insufficient to recover the exact state of affairs, but improvement of their non-local spatial structure toward the exact adiabatic kernel can assist matters.

Two xc kernels that differ in structure that is predominantly captured with the gauge transform defined in Section II.3 provides an explanation for the approximate agreement between the derived properties of the two xc kernels in question. This possibility is now considered, and in fact we shall demonstrate that the gauge freedom of $f_{\text{xc}}$ is not sufficient to explain the similarity observed between, for example, the optical properties calculated using the adiabatic LDA and RPA in Section III. Moreover, the particular form of the non-local spatial structure within the exact $f_{\text{xc}}(x,x^{\prime},\omega)$ is in general not possible to capture with this gauge freedom, and it is unlikely that this line of reasoning is able to explain the efficacy of approximations of any kind.

In order to demonstrate this, an optimal gauge is defined that transforms, in as much as is possible, one xc kernel into another using the gauge freedom, i.e. it brings one $f_{\text{xc}}$ toward another $f_{\text{xc}}$ in a particular matrix norm. The definition and derivation of the optimal gauge is given in the relevant section of the supplemental material. Furthermore, an illustration of the optimal gauge transform for each of the examples to follow is also provided in the supplemental material. The atom of Section III is considered, and in particular the optimal gauge is found in order to match the RPA, $f_{\text{xc}}^{\text{RPA}}=0$, with the adiabatic LDA used in this work, $f_{\text{xc}}^{\text{ALDA}}$ [65]. As discussed in Section II.3, the vectors $g$ and $h$ introduce an unavoidable degree of non-locality, and the local spatial structure of the adiabatic LDA simply cannot be reproduced with a gauge transform of this kind applied to the RPA.

The story remains much the same when an attempt is made to match the adiabatic LDA xc kernel to the exact adiabatic xc kernel for the atom, $f_{\text{xc}}(x,x^{\prime},\omega=0)$. The particular form of the spatial non-locality in the exact adiabatic xc kernel does not lend itself to the fairly inflexible structural freedom afforded by the functions $g$ and $h$. This conclusion holds true for almost all the xc kernels examined in this work, namely, the spatial profile of the exact $f_{\text{xc}}$ at any $\omega$ is not possible to reproduce with a gauge transform applied to the adiabatic LDA or RPA, despite, at high frequencies, the similarity of the derived optical spectra from these approximations.

A final remark is given in relation to the divergences in $f_{\text{xc}}$ studied in the previous section. The xc kernel around a divergence is approximately of the form of an outer product $|u\rangle\langle u|$, where $|u\rangle$ is the eigenvector of $f_{\text{xc}}$ that is diverging. The xc kernel is therefore mostly composed of $N$ degrees of freedom, a state of affairs that is a priori much more amenable to the gauge freedom. In fact, the first pole in the $f_{\text{xc}}$ of the atom, i.e. the beginning of the non-adiabatic behavior, see panel (d) of Fig. 3, is non-divergent after the optimal gauge is applied. This divergence therefore does not affect observable properties such as its associated optical peak, and indeed the exact adiabatic kernel is able to capture this peak, despite it existing squarely within the divergence. However, such a situation is not detected again for the remainder of the divergences.

In finite systems, it is the conventional wisdom that an accurate ground-state xc potential is more important for capturing the interacting excitation spectrum than a sophisticated $f_{\text{xc}}$ [16, 15, 68, 69, 70]. The effect of an improved treatment of exchange and correlation in the ground state is shown in Fig. 12, which demonstrates the optical spectrum for the atomic system calculated from the non-interacting response function $\chi_{0}$ at various levels of approximation. The transition energies and rates calculated from the exact Kohn-Sham system are considerably closer to the exact transitions than, for example, Hartree theory, which is increasingly poor at higher energies.

These non-interacting response functions can be used, in conjunction with their corresponding $f_{\text{xc}}$, to solve the Dyson equation, thus shifting the position and weight of the peaks. An xc kernel functional corresponds to an xc potential functional if it is the second functional derivative of the xc potential with respect to the density. Thus, the RPA xc kernel corresponds to Hartree theory, and the adiabatic LDA xc kernel corresponds to the ground-state LDA from which it came. To use an xc kernel that does not correspond to the ground-state xc potential can violate various exact conditions [71] – this is evidently the case here for the zero-force sum rule.

The aforementioned conventional wisdom is exhibited clearly in Fig. 12, namely, use of a corresponding $f_{\text{xc}}$ is only able to improve matters slightly beyond the non-interacting peaks, which is most visible at higher frequencies. For this reason, even when using an incompatible $f_{\text{xc}}$, the calculation benefits from an improved treatment of ground-state exchange and correlation. This is evidenced in the case of the atom in Section III, where using the RPA and adiabatic LDA xc kernels in conjunction with the exact Kohn-Sham ground-state yields a more accurate optical spectrum than if one were to attempt to keep the xc potential and xc kernel compatible.