# Polarizable Continuum Models and Green's Function GW Formalism: On the Dynamics of the Solvent Electrons

Ivan Duchemin, David Amblard, Xavier Blase

► To cite this version:

Ivan Duchemin, David Amblard, Xavier Blase. Polarizable Continuum Models and Green's Function GW Formalism: On the Dynamics of the Solvent Electrons. Journal of Chemical Theory and Computation, 2024, 20 (20), pp.9072-9083. 10.1021/acs.jctc.4c00745 . hal-04752298

HAL Id: hal-04752298

https://hal.science/hal-04752298v1

Submitted on 24 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Polarizable Continuum Models and Green's Function GW Formalism: On the Dynamics of the Solvent Electrons

Ivan Duchemin,$^{*,†}$ David Amblard,$^{‡}$ and Xavier Blase$^{*,‡}$

$^{†}$Univ. Grenoble Alpes, CEA, IRIG-MEM-L_Sim, 38054 Grenoble, France

$^{‡}$Univ. Grenoble Alpes, CNRS, Inst NEEL, F-38042 Grenoble, France

E-mail: ivan.duchemin@cea.fr; xavier.blase@neel.cnrs.fr

## Abstract

The many-body GW formalism, for the calculation of ionization potentials or electronic affinities, relies on the frequency-dependent dielectric function built from the electronic degrees of freedom. Considering the case of water as a solvent treated within the polarizable continuum model, we explore the impact of restricting the full frequency-dependence of the solvent electronic dielectric response to a frequency-independent ($\epsilon_\infty$) optical dielectric constant. For solutes presenting small to large highest-occupied to lowest-unoccupied molecular orbital energy gaps, we show that such a restriction induces errors no larger than a few percent on the energy level shifts from the gas to the solvated phase. We further introduce a remarkably accurate single-pole model for mimicking the effect of the full frequency dependence of the water dielectric function in the visible-UV range. This allows a fully dynamical embedded GW calculation with the only knowledge of the cavity reaction field calculated for the $\epsilon_\infty$ optical dielectric constant.
# 1 Introduction

The renormalization of electronic excitations by a polarizable environment, such as an electrode, an organic crystal or a solvent, plays an important role in many fields pertaining to physics, chemistry, and biology. This is particularly true in the case of charged excitations, as described by direct or inverse photoemission, where the reaction field from the environment can stabilize the added charge by several electronvolts. The stabilization energy, labeled a polarization energy, depends on the environment dielectric properties. The complexity of the environment, in many applications of interest, initiated hybrid theoretical strategies, merging the quantum mechanical treatment of a central active subsystem, with a simplified description of the environment response properties described as a continuum or discrete (atomistic) medium.¹,²

Originating in simple image models for the description of a charge facing a dielectric surface or located inside a dielectric cavity,³,⁴ the polarizable continuum model (PCM) in its various implementations has emerged as a very successful approach,²,⁵⁻¹¹ in particular in the description of solvated species for which more accurate but expensive atomistic simulations require averaging over molecular dynamics trajectories. In the PCM approach, the solute is placed in a cavity carved into an homogeneous medium characterized by its macroscopic dielectric function $\epsilon(\omega)$ that depends on general principles on the excitation frequency. Of importance for the following discussions, two specific solvent dielectric constants can be introduced when it comes to study fast electronic excitations, namely the static and optical dielectric constants labeled $\epsilon_0$ and $\epsilon_\infty$, respectively. The first one (e.g., $\epsilon_0=78.35$ for water) includes the slow ionic motions, while the second ($\epsilon_\infty=1.78$ for water) only accounts for the fast electronic degrees of freedom. The use of different dielectric constants in different frequency domains leads to the so-called non-equilibrium formalisms.⁷ In practice, within this very simplified description of the frequency dependence of the dielectric constant, $\epsilon_0$ is employed to renormalize the ground-state properties (charge density, one-body molecular orbitals, etc.) while $\epsilon_\infty$ is adopted for the solvent reaction to the fast electronic excitations from the ground-state to an excited state.

The solvent is chosen to be mostly transparent (non-absorbing) in the energy range of the solute
lowest electronic excitations of interest. As such, the poles of the solvent electronic susceptibility are expected to be located at significantly higher energy. This may certainly justify the use of a frequency-independent $\epsilon_\infty$ optical dielectric constant defined from the static limit of the solvent electronic susceptibility. The subscript $\infty$ really points to a frequency range significantly larger than the typical ionic frequencies, even though smaller than the energy of the solvent electronic susceptibility poles. Restricting the PCM dielectric response to the $\epsilon_0$ and $\epsilon_\infty$ dielectric constants can be improved by considering the full frequency-dependent $\epsilon(\omega)$ macroscopic dielectric constant for frequencies spanning the typical ionic to electronic time scales. This was introduced within a time-dependent non-equilibrium dielectric response PCM approach.<sup>12–18</sup> Concerning the solvent/solute problem, the frequency dependence often relied on the Debye relaxation model for the solvent degrees of freedom dynamics, interpolating between $\epsilon_0$ and $\epsilon_\infty$.<sup>12,15–17</sup> As such, the high-frequency limit of the dielectric response was set to $\epsilon_\infty$. In practice, and as shown in Fig. 1(a) (blue line) for water as the solvent, the macroscopic optical dielectric function, the square of the refractive index, shows significant variations with frequency in the visible-UV range, and $\epsilon_\infty=1.78$ for water corresponds to the low frequency limit response from the electronic degrees of freedom. In this spirit, a fit to experimental measurements of the frequency-dependent dielectric function in the visible-UV range was exploited for the time-dependent study of the reaction field close to metallic nanoparticles treated as a continuous polarizable medium.<sup>18</sup> Besides specific formulations of the frequency-dependent macroscopic dielectric function, a general PCM approach with a fully frequency-dependent solvent response matrix was introduced within an open quantum system theory framework.<sup>19</sup>

As another issue, the comparison between the dynamics of the solvent and solute electronic degrees of freedom leads to two limiting regimes. If the low lying poles of the solute electronic susceptibility are located at much lower energy than that of the solvent, then the solvent electronic response to the solute low-lying excitations can be considered to be instantaneous. This is called the Born-Oppenheimer (BO) regime in the PCM literature.<sup>19–22</sup> The opposite regime, where the solvent electronic degrees of freedom are assumed to be slower than that of the solute, is called the
self-consistent regime. Clearly, in cases where the decoupling of energy between the solute and solvent electron dynamics cannot be made, none of these regime strictly applies. In particular, the actual frequency dependence of the solvent optical dielectric response needs to be considered.

Already merged with the PCM formalism,²³⁻²⁵ the Green’s function GW many-body perturbation theory,²⁶ where *G* stands for the one-body Green’s function and *W* for the screened Coulomb potential, has recently gained much popularity in quantum chemistry. Significant efforts have been devoted to benchmark its accuracy for molecular systems as compared to quantum chemistry techniques such as coupled-cluster approaches,²⁷⁻³⁹ together with much progress in developing low-scaling implementations.⁴⁰⁻⁵⁴ Targeting the calculation of electronic energy levels, including ionization potentials and electronic affinities, properly defined as charging energies, the GW operator can be described as a non-local and dynamical self-energy term that includes exchange and correlation effects. In particular, the screened Coulomb potential *W*(*r*, *r*'; *ω*) is fully dynamical, accounting for the frequency dependence in the visible-UV range of the dielectric function *ε*(*r*, *r*'; *ω*). As such, the GW formalism offers an obvious path to explore fully dynamical models of polarizable environments, with in particular a macroscopic dielectric function originating from electronic degrees of freedom that is not restricted to a single *ε*<sub>∞</sub> dielectric constant, but is instead allowed to depend on the excitation frequency in the visible-UV range.

Replacing the fully dynamical optical dielectric function of Fig. 1(a) by the frequency-independent ε∞ constant, amounts to pushing the poles of the solvent electronic susceptibility to infinity. This leads to the BO regime where the solvent electronic degrees of freedom are assumed to have a much faster dynamics than that of the solute. Such an analysis is consistent with the conclusions provided by Guido and coworkers in Ref. 19 concerning the early merging²³ of the GW formalism with a PCM optical dielectric function fixed to the ε∞ constant. Such an approximation of a much faster solvent electrons dynamics is supposed to be valid when the solvent transition energies from occupied to unoccupied states are located at significantly higher energy than the solute analogs. Reintroducing the full frequency-dependence of the solvent optical dielectric response, as provided in Fig. 1 for water, allows to tackle systems where solute and solvent electron dynamics are not
decoupled, including situations where the poles of the solute electronic susceptibility may be located
at higher energies than that of the solvent. Recently, a fully dynamical PCM has been introduced
and integrated with the GW formalism,²⁵ but no comparison was made to an approach where the
frequency-dependent optical dielectric function is restricted to the ε∞ constant.

In a recent work,⁵⁵ the authors merged the GW formalism with a dynamical treatment of
the environment degrees of freedom in a fully *ab initio* QM/QM' scheme. These developments
allowed in particular to scrutinize the accuracy of the instantaneous environment electrons response
approximation.⁵⁶ Considering a molecular system immersed in its parent molecular crystal, that is
a situation where there is no decoupling of energy between the “solute” and “solvent” electronic
degrees of freedom, the error associated with treating the environment in the BO limit was shown
to induce errors no larger than 10% for the polarization energy associated with frontier orbitals
or the energy gap. Such an *ab initio* and dynamical treatment of the environment stands as a
generalization of previous QM/MMpol implementations²³,⁵⁷⁻⁶¹ based on semi-empirical and low-
frequency descriptions of the environment dielectric properties in the optical range.⁶²⁻⁶⁴

In the present study, we explore explicitly the effect of considering the full frequency dependence
of the optical dielectric constant of water as a solvent described by the PCM. More explicitly, we
calculate the GW ionization potential and electronic affinities of solvated molecules, switching on
and off the frequency-dependence of the optical macroscopic dielectric constant of water. For a
large set of molecules showing very different highest-occupied to lowest-unoccupied molecular
orbital (HOMO-LUMO) energy gaps, we show that treating the solvent in the BO limit does not
induce errors larger than a few percent on the polarization energy. We further introduce a simple
pole-model for the solvent electronic dielectric response that accurately reproduces the effect of
considering the full dynamics of the solvent electronic degrees of freedom, while only requiring the
calculation of the PCM reaction field associated with the ε∞ optical dielectric constant.
## 2 Theory

We start this theory section with a brief outline of the specific GW features and flavors that will be used hereafter to support strategies for merging with a polarizable environment. Broader and extensive descriptions of the Hedin's formalism can be found elsewhere, as we refer the reader to either seminal articles$^{26,65–69}$ or more recent books or reviews.$^{70–76}$ We will keep our demonstration to the quantities relevant to the present work.

### 2.1 The GW formalism

A standard GW calculation starts from an input time ordered Green's function $G(\mathbf{r}, \mathbf{r}'; \omega)$, built upon $\{\epsilon_n, \phi_n\}$ Kohn-Sham eigenstates, and constructs the dynamical GW exchange-correlation self-energy operator defined as

$$ \Sigma(\mathbf{r}, \mathbf{r}'; E) = \frac{i}{2\pi} \int d\omega \, e^{i\eta\omega} G(\mathbf{r}, \mathbf{r}'; E + \omega) W(\mathbf{r}, \mathbf{r}'; \omega), \quad (1) $$

with $\eta$ a positive infinitesimal. The dynamically screened Coulomb potential $W(\mathbf{r}, \mathbf{r}'; \omega)$ comes instead of the bare Coulomb potential $v(\mathbf{r}, \mathbf{r}')$ that would lead to the Hartree-Fock exact exchange $Gv$ operator. As the self-energy is dynamical, performing a GW calculation in a polarizable or dielectric environment relies straightforwardly on the relation between $W(\mathbf{r}, \mathbf{r}'; \omega)$ and the inverse dielectric function $\epsilon^{-1}(\mathbf{r}, \mathbf{r}'; \omega)$:

$$ \begin{aligned} W(\mathbf{r}, \mathbf{r}'; \omega) &= \int d\mathbf{r}'' \, \epsilon^{-1}(\mathbf{r}, \mathbf{r}''; \omega) v(\mathbf{r}'' , \mathbf{r}') \\ &= \epsilon^{-1}(\omega) \cdot v, \end{aligned} \quad (2) $$

where we introduce the $\cdot$ product to indicate composition over the space variables of two adjacent linear operators, i.e. the result of their successive application as defined by the matrix product of their matrix representations. Since the GW formalism tackles the study of electronic correlation effects, we emphasize that the dielectric function $\epsilon^{-1}(\omega)$ considered here shall only account for the
electronic degrees of freedom, justifying below the notation $\epsilon_{opt}$ where “opt” stands for optical. In contrast, the effect of the slow ionic degrees of freedom was accounted for at the Density Functional Theory (DFT) level, e.g., generating initial DFT@PCM Kohn-Sham states with a dielectric constant of $\epsilon_0 = 78.35$ for water as the solvent.⁷⁷

Alternatively, one can define $W(\mathbf{r}, \mathbf{r}'; \omega)$ through the electronic susceptibility $\chi(\mathbf{r}, \mathbf{r}'; \omega)$:

$$W(\omega) = v + v \cdot \chi(\omega) \cdot v, \qquad (3)$$

which is itself related to the free susceptibility $\chi_0(\mathbf{r}, \mathbf{r}'; \omega)$:

$$\chi(\omega) = \chi_0(\omega) + \chi_0(\omega) \cdot v \cdot \chi(\omega), \qquad (4)$$

adopting the Random Phase Approximation (RPA). All the above quantities are dynamical, and require in particular the knowledge of the full dynamics of $\epsilon_{opt}^{-1}(\mathbf{r}, \mathbf{r}'; \omega)$. As such, merging with low frequency limit models such as PCM cannot be made straightforwardly as taking $\epsilon_\infty = \epsilon_{opt}(\omega \to 0)$ for the environment would lead to the wrong high frequency limit for $W(\mathbf{r}, \mathbf{r}'; \omega)$.

The difficulties brought by the frequency dependence within the GW self-energy of eq 1 were circumvented by Hedin²⁶ through the introduction of the so-called Coulomb-hole (COH) plus screened-exchange (SEX) static approximation:

$$\Sigma^{\text{SEX}}(\mathbf{r}, \mathbf{r}') = - \sum_{i}^{\text{occ}} \phi_i(\mathbf{r}) \phi_i^*(\mathbf{r}') W(\mathbf{r}, \mathbf{r}'; \omega = 0) \qquad (5)$$

$$\Sigma^{\text{COH}}(\mathbf{r}, \mathbf{r}') = \frac{1}{2} \sum_{n} \phi_{n}(\mathbf{r}) \phi_{n}^{*}(\mathbf{r}') [W(\mathbf{r}, \mathbf{r}'; \omega = 0) - v(\mathbf{r}, \mathbf{r}')] , \qquad (6)$$

requiring only the low-frequency $W(\omega \to 0)$ limit of the screened Coulomb potential, or equivalently of the susceptibility $\chi(\omega \to 0)$. An elegant and simple way to recover the static COHSEX approximation is to consider that the system’s electron density reacts immediately to any perturbation
within the GW equations response functions 3 and 4.⁵⁵ In this case, the dynamical response of the system may be described by a simple pole model tailored to reproduce this low frequency limit:

$$ \chi_{\Omega_p}(\mathbf{r}, \mathbf{r}'; \omega) = \chi(\mathbf{r}, \mathbf{r}'; \omega = 0) \times f(\omega; \Omega_p) \quad (7) $$

with

$$ f(\omega; \Omega_p) = \frac{\Omega_p}{2} \left[ \frac{1}{\omega + \Omega_p - i\eta} - \frac{1}{\omega - \Omega_p + i\eta} \right], \quad (8) $$

and $\Omega_p$ a unique pole energy for simplicity. Such an expression has the correct static limit and time-ordering structure in the energy plane. The static COHSEX self energy 5+6, is recovered by taking the pole energy $\Omega_p$ to infinity after performing the energy integration in eq 1. We will exploit here below the same strategy in order to merge the standard PCM model, described by a single $\epsilon_{\text{opt}}(\omega \to 0)$ dielectric constant, within a fully dynamical GW formalism for the central subsystem. We will also further explore the possibility to efficiently account for the solvent dynamics through such a simple pole model but with a finite $\Omega_p$ pole energy.

## 2.2 Embedding the GW equations

### 2.2.1 General considerations

We start with a partition of the system into the solute S and the solvent (i.e., environment) E, and assume independent contributions of the solute and solvent free susceptibilities $\chi_0^S$ and $\chi_0^E$ to the free susceptibility $\chi_0$ of the full system. This assumption is fully justified in the case of non-overlapping wave-functions between subsystems S and E. In such case, the screened Coulomb potential can be obtained through:

$$ \begin{aligned} W^{-1}(\omega) &= v^{-1} - (\chi_0^S(\omega) + \chi_0^E(\omega)) \\ &= \tilde{v}^{-1}(\omega) - \chi_0^S(\omega), \end{aligned} \quad (9) $$
with the frequency dependent screened Coulomb interactions $\tilde{v}(\omega)$ defined through the Woodbury identity:

$$ \begin{aligned} \tilde{v} &= (\nu^{-1} - \chi_0^E(\omega))^{-1} \\ &= \nu + \nu \cdot \chi^E(\omega) \cdot \nu, \end{aligned} \qquad (10) $$

and $\chi^E(\omega)$ is the interacting susceptibility of the environment, or solvent, in the absence of the solute

$$ \chi^E(\omega)^{-1} = \chi_0^E(\omega)^{-1} - \nu. \qquad (11) $$

For clarity in the following equations, we will write

$$ \nu_{\text{reac}}(\omega) = \nu \cdot \chi^E(\omega) \cdot \nu \qquad (12) $$

the reaction potential of the solvent induced by the solute: qualitatively, a change of the solute charge density induces a change in the solvent charge density proportional to $\chi^E \cdot \nu$ that in return exerts a reaction potential on the solute via $\nu$. Finally, the solution to eq 9 can be decomposed in

$$ W(\omega) = \tilde{v}(\omega) + \tilde{v}(\omega) \cdot \tilde{\chi}(\omega) \cdot \tilde{v}(\omega) \qquad (13) $$

and

$$ \tilde{\chi}(\omega) = \chi_0^S(\omega) + \chi_0^S(\omega) \cdot \tilde{v}(\omega) \cdot \tilde{\chi}(\omega), \qquad (14) $$

which follows the ingredients of a standard isolated GW calculation, but where $\tilde{\chi}(\omega)$ is the interacting susceptibility of the solute obtained with the frequency dependent Coulomb interactions that have been screened through the solvent reaction potential.
### 2.2.2 Dynamical PCM models

The dynamical cavity reaction field $v_{\text{reac}}(\omega)$ can be seen as the result of a frequency-dependent PCM model associated with the frequency-dependent $\epsilon_{\text{opt}}(\omega)$ solvent optical dielectric constant. The PCM formalism can thus be straightforwardly generalized to a fully dynamical scheme, recalculating the cavity reaction fields for each frequency involved in the energy-quadrature used to evaluate the self-energy (see eq 1). In our approach combining the contour-deformation scheme with the analytic continuation of the screened Coulomb potential (see Technical details section 2.3),⁷⁸ this means recalculating the cavity charges for typically twelve imaginary frequencies. Since the dielectric constant is real along the imaginary axis, the standard continuity equations on the cavity surface⁹,²³ can be applied without modifications.

The frequency dependence of the water macroscopic optical dielectric constant was recently parametrized under the form of damped harmonic oscillators reproducing faithfully available experimental data.⁷⁹ The resulting frequency-dependent optical dielectric constant $\epsilon_{\text{opt}}(\omega)$ is represented in Fig. 1(a) along the real-axis (blue line) and along the imaginary-axis (green line). We further plot the related $[1/\epsilon_{\text{opt}}(\omega) - 1]$ along the imaginary axis (green line, Fig. 1(b)) and along the real axis (blue line, Fig. 1(c)). As inferred from eqs 2 and 3, this quantity is more closely related to the susceptibility $\chi$ entering the reaction field:

$$ \epsilon_{\text{opt}}(\omega)^{-1} = 1 + \nu \cdot \chi(\omega). \qquad (15) $$

In particular, $\epsilon_{\text{opt}}(\omega)^{-1}$ and $\chi(\omega)$ share the same pole structure.

The fit of Ref. 79 relies on damped oscillators. This is a difficulty since the resulting susceptibility, or related $\epsilon_{\text{opt}}(\omega)^{-1}$ inverse macroscopic dielectric constant, does not offer the proper symmetry along the real and imaginary-frequency axes. We bypass this problem by showing in Fig. 1 that a simple single-pole model with proper symmetry:
<figure><img src="image_1.png" /><figcaption>Figure 1: (a) Real-part of the macroscopic optical dielectric constant of water ($\epsilon_{\text{opt}}$) along the real (blue line) and imaginary (green line) frequency axes [from Ref. 79]. The $\omega = 0$ limit allows to recover the $\epsilon_{\infty} = 1.78$ optical dielectric constant. The dashed purple line represents the dielectric constant at imaginary frequencies, as obtained with the one-pole plasmon-pole model (PPM) described in the main text (eq 16). (b) Plot of $[1/\epsilon_{\text{opt}}(i\omega) - 1]$ from Ref. 79 (green line) and corresponding single-pole functional fit with a pole located at $\Omega_p=21$ eV (dashed purple line). The red dashed line represents the limit case where $\Omega_p$ is pushed to infinity. (c) Real-part of $[1/\epsilon_{\text{opt}}(\omega) - 1]$ from Ref. 79 (blue line) and using the single-pole functional fit with a pole at 21 eV and a 0.5 eV broadening (dashed purple line).</figcaption></figure>
**dynamical PCM(ω):**

$$ \epsilon_{\text{opt}}(\omega)^{-1} = 1 + (\epsilon_{\infty}^{-1} - 1) \times f(\omega; \Omega_p) \quad (16) $$

offers an accurate approximation to the full fit along the imaginary axis, fixing the pole energy at $\Omega_p=21$ eV. The resulting $[1/\epsilon_{\text{opt}}(\omega) - 1]$ function (dashed purple line, Fig. 1(b)) follows closely the original fit along the imaginary axis. A similar agreement is found for the related $\epsilon_{\text{opt}}(\omega)$ function (dashed purple line, Fig. 1(a)). Even though not needed for our analytic continuation scheme applied to frontier orbitals, we show in Fig. 1(c) that the single-pole fit is accurate at low energy along the real-energy axis. As shown in the Supporting Information (SI, Figs. S1 and S2), a more accurate fit, allowing for several poles, does not change the resulting polarization energies by more than a very few meV. Such a scheme, and the related polarization energies, represent our reference dynamical PCM calculations and will be labeled PCM(ω).

Fully dynamical PCM calculations, requiring to recalculate the cavity reaction fields at several frequencies, increase the complexity and computational cost as compared to the BO limit (fast solvent electron dynamics) where only the $\epsilon_{\infty}$ optical dielectric constant needs to be considered. We now introduce a seemingly severe approximation by modeling the dynamics of the non-local susceptibility as follows:

**plasmon-pole model solvent (PCMPP):**

$$ \chi_{\Omega_p}^E(\mathbf{r}, \mathbf{r}'; \omega) = \chi_{\text{PCM}}^E(\mathbf{r}, \mathbf{r}'; \omega = 0) \times f(\omega; \Omega_p). \quad (17) $$

While the fit of the experimental macroscopic $\epsilon_{\text{opt}}(\omega)^{-1}$ function by a single-pole function (eq 16) represents just a properly symmetrized representation of a frequency-dependent scalar, the latter approximation is much more questionable since it introduces a decoupling between the real-space
and frequency degrees of freedom. On formal grounds, the solvent susceptibility reads:

$$ \chi^E(\mathbf{r}, \mathbf{r}'; \omega) = \sum_n^{\text{poles}} \left[ \frac{A_n(\mathbf{r}, \mathbf{r}')}{\omega - \Omega_n + i\eta} - \frac{A_n(\mathbf{r}, \mathbf{r}')}{\omega + \Omega_n - i\eta} \right], $$

with $A_n(\mathbf{r}, \mathbf{r}') = \rho_n(\mathbf{r})\rho_n(\mathbf{r}')$ and where the $\Omega_n$ are the RPA transition energies and $\rho_n(\mathbf{r})$ the corresponding transition densities. Clearly, the spatial and frequency degrees of freedom are entangled. The decoupling of the spatial and frequency degrees of freedom for the solvent susceptibility extends to the reaction field $v_{\text{reac}}(\omega) = v \cdot \chi^E(\omega) \cdot v$. As such, the reaction field at finite frequency can be straightforwardly obtained by multiplying the static reaction field by the dynamical factor $f(\omega, \Omega_p)$. This is a dramatic simplification since only the low-frequency limit of the reaction fields needs to be calculated as in a standard PCM calculation. A similar approach was recently proposed to improve on the static COHSEX approximation for calculating the IP and EA of molecular systems.⁸⁰ Our approach is less ambitious as we only intend to capture the dynamics of the environment susceptibility. The solute electrons correlation energy is in our case treated at the fully dynamical GW level.

The validation of this plasmon-pole approximation for the solvent susceptibility, as compared to the dynamical PCM approach where the reaction fields are explicitly recalculated for each $\epsilon_{\text{opt}}(\omega)$ value, represents one of the main goal of this paper. This approach will be labeled PCMpp in what follows, where PP stands for plasmon-pole.

### 2.2.3 The static COHSEX PCM model: $\Omega_p \to \infty$

To connect to PCM calculations that neglect the frequency dependence of the optical dielectric constant, and in order to explore the consequences of such an approximation as compared to the fully dynamical PCM, we now turn to GW@PCM calculations in the BO limit, namely assuming that the solvent degrees of freedom react instantaneously to any excitation on the solute. As shown recently in a fully ab initio QM/QM' implementation of embedded GW calculations,⁵⁵ this can be very simply obtained from eq 17 by bringing the $\Omega_p$ frequency to infinity after performing the frequency integration that defines the GW self-energy (eq 1). This leads to the so-called static
Coulomb-Hole plus Screened-Exchange (COHSEX) decomposition of the reaction field as briefly described now. More details can be found in Ref. 55.

We focus on the correlation-only $\Sigma_c(E)$ self-energy, leaving aside the bare exchange contribution, with:

$$ \Sigma_c(E) = \frac{i}{2\pi} \int d\omega \, e^{i\eta\omega} G(E + \omega) [W(\omega) - v], \quad (18) $$

where the space variables are omitted for the sake of compactness. This equation can be rewritten in terms of the embedded Coulomb interactions and reaction field as

$$ \Sigma_c(E) = \frac{i}{2\pi} \int d\omega \, e^{i\eta\omega} G(E + \omega) [W(\omega) - \tilde{v}(\omega)] + \frac{i}{2\pi} \int d\omega \, e^{i\eta\omega} G(E + \omega) v_{\text{reac}}(\omega). \quad (19) $$

It turns out that, beyond its similarity with the usual formulation 18, the first term of eq 19 brings a contribution to the self energy that can be computed using only the low frequency limit of $\tilde{v}(\omega)$. This is thanks to the contribution of $W(\omega) - \tilde{v}(\omega)$ poles coming from the solvent dynamics vanishing as $\Omega_p \to \infty$. Thanks to that, and as detailed in appendix A, we can write:

$$ \begin{aligned} & \lim_{\Omega_p \to \infty} \frac{i}{2\pi} \int d\omega \, e^{i\eta\omega} G(E + \omega) [W(\omega) - \tilde{v}(\omega)] \\ &= \frac{i}{2\pi} \int d\omega \, e^{i\eta\omega} G(E + \omega) [\tilde{W}(\omega) - \tilde{v}(0)] \end{aligned} \quad (20) $$

with

$$ \tilde{W}(\omega)^{-1} = \tilde{v}(0)^{-1} - \chi_0^S(\omega). \quad (21) $$

Remains thus only the contribution from the finite frequency poles of the solute, that appears as re-normalized within the static limit of the solvent. The direct contribution of the poles of the
environment response comes from the second term of eq 19:

$$ \begin{aligned} & \frac{i}{2\pi} \int d\omega e^{i\eta\omega} G(\mathbf{r}, \mathbf{r}'; E + \omega) v_{\text{reac}}(\mathbf{r}, \mathbf{r}'; \omega) \\ & = - \sum_{i}^{\text{occ}} \phi_i(\mathbf{r}) \phi_i^*(\mathbf{r}') v_{\text{reac}}(\mathbf{r}, \mathbf{r}'; \omega = 0) \\ & \quad + \frac{1}{2} \sum_{n} \phi_n(\mathbf{r}) \phi_n^*(\mathbf{r}') v_{\text{reac}}(\mathbf{r}, \mathbf{r}'; \omega = 0) \end{aligned} \quad (22) $$

where (i) runs over occupied states only. The first line takes the form of a screened-exchange-like contribution to the reaction field, the second one is the analog of the Coulomb-hole (COH) contribution. Using directly the static reaction field $v_{\text{reac}}(\omega = 0)$ without the pole structure coming from $f(\omega; \Omega_p)$ leads to the SEX contribution only, leaving aside the COH term. As shown below, this second term is crucial when looking at the polarization energy associated with individual energy levels.

We finally have an expression that differentiates between indirect and direct solvent contributions to the embedded GW quasiparticle self energy:

**Static COHSEX PCM (BO limit):**

$$ \begin{aligned} \Sigma_{\text{C}}(\mathbf{r}, \mathbf{r}'; E) = & \frac{i}{2\pi} \int d\omega e^{i\eta\omega} G(\mathbf{r}, \mathbf{r}'; E + \omega) [\tilde{W}(\mathbf{r}, \mathbf{r}'; \omega) - \tilde{v}(\mathbf{r}, \mathbf{r}'; 0)] \\ & - \sum_{i}^{\text{occ}} \phi_i(\mathbf{r}) \phi_i^*(\mathbf{r}') v_{\text{reac}}(\mathbf{r}, \mathbf{r}'; \omega = 0) \\ & + \frac{1}{2} \sum_{n} \phi_n(\mathbf{r}) \phi_n^*(\mathbf{r}') v_{\text{reac}}(\mathbf{r}, \mathbf{r}'; \omega = 0). \end{aligned} \quad (23) $$

The last two lines can be regarded as perturbative corrections to the gas phase self-energy. However, the first line integral also differs from its gas phase analog through the renormalization of the screened Coulomb and bare Coulomb potential by the static reaction field. Each term gives rise to a contribution to the polarization energy through the impact of the PCM on the correlation energy.
Overall, the polarization energy reads:

$$P_n = \epsilon_n^{GW@PCM} - \epsilon_n^{GW@gas}$$

where $\epsilon_n^{GW@PCM}$ and $\epsilon_n^{GW@gas}$ are respectively the embedded (PCM) and gas phase quasiparticle energies. We emphasize again that in a GW@PCM calculation, the input Kohn-Sham eigenstates are calculated at the DFT@PCM($\epsilon_0$) level.

## 2.3 Technical details

Our calculations are performed with the BEDEFT (beyond-DFT) package$^{52,78}$ implementing the GW and Bethe-Salpeter equation (BSE) formalisms using Gaussian basis sets and Coulomb-fitting (RI-V) resolution-of-the-identity.$^{81–83}$ The correlation part of the self-energy is calculated adopting a contour-deformation scheme with an integration performed along the imaginary-frequency axis, completed by residues involving the value of the screened Coulomb potential along the real-axis. The values of W along the real-energy axis can be accurately obtained by analytic continuation, a scheme much more stable than the direct analytic continuation of the full GW self-energy.$^{78}$ GW calculations are performed at the partially-self-consistent evGW level where quasiparticle energies are reinjected self-consistently in the construction of the solute Green’s function G and independent-particle susceptibility $\chi_0$. We adopt the def2-TZVP basis set$^{84}$ together with the corresponding def2-TZVP-RIFIT auxiliary basis set.$^{85}$ Input Kohn-Sham eigenstates are generated using the PBE0 functional.$^{86,87}$

Our integral-equation-formalism (IEF)$^9$ PCM formulation is described in Ref. 23, where the reaction field is obtained through a double layer potential version of the formalism: the reaction field $v_{reac}$ corresponds to the potential experienced by the solute coming from the surface charges and surface dipoles induced on the cavity walls. This double layer potential has the merit of handling correctly solute charge spilling out of the cavity, which turns out to be an important point when considering virtual states. As in the standard (IEF)-PCM framework, the model depends
implicitly on the solvent dielectric constant through the continuity equations at the cavity boundaries.
Moreover, $v_{\text{reac}}$ is calculated in the auxiliary basis $\{P\}$ used to describe charge density variations
in our Coulomb-fitting scheme. Namely, we calculate $v_{\text{reac}}(P, Q; \omega)$ matrix elements which can be
interpreted as the action on $P$ of the reaction field associated with the surface charges and dipoles
generated by the $Q$ charge density in the cavity for a given $\epsilon_{\text{opt}}(\omega)$ macroscopic dielectric constant.
Once the reaction-field matrix is obtained in the auxiliary basis, its action on any Kohn-Sham state
is straightforward. In particular, there is no need to recalculate the reaction field components at each
GW iteration in the case of a self-consistent GW scheme.

Input Kohn-Sham eigenstates are generated with the ORCA package$^{88}$ that implements the C-PCM version$^{10}$ of the PCM approach. Embedded DFT calculations with the water environment are thus performed using the $\epsilon_0=78.355$ dielectric constant. As such, all GW calculations with the PCM environment start with the same DFT@PCM($\epsilon_0$) eigenstates. The differences between the GW@PCM$_{COHSEX}$ and dynamical PCM polarization energies can only stem from the way the dynamical dependence of $\epsilon_{\text{opt}}(\omega)$ is treated.

**List of molecules:** We study a large set of solvated molecules, starting with the adenine, cytosine, thymine, and uracil nucleobases studied in a first implementation of the GW@PCM formalism in the standard static PCM limit.$^{23}$ Furthermore, as studied in the first merging of the Bethe-Salpeter equation (BSE) with the PCM,$^{89}$ we select acrolein and indigo, the push-pull p-nitro-aniline (PNA) molecule in its planar and rotated (perpendicular) conformation, the donor-acceptor benzene/TCNE complex, and the 4-Nitropyridine N-oxide organic probe commonly used to assess solvent polarities. Geometries are taken from Refs. 23,89. Finally, acetaldehyde, ethanol, formaldehyde, and water molecule and hexamer clusters are selected to offer a large range of HOMO-LUMO gaps, from 5.66 eV for indigo, 13.14 eV for formaldehyde, to 16.21 eV for the water monomer at the evGW@PBE0/def2-TZVP (gas phase) level. Selecting such a large spread of HOMO-LUMO gaps is expected to serve as a test of the BO approximation (faster solvent electron dynamics) since such an approximation may be expected to be more robust if the gap of the solute is much smaller than that of the solvent. The water cluster geometries are taken from Refs. 90,91
while the acetaldehyde, ethanol, and formaldehyde molecules are relaxed at the B3LYP/def2-TZVP level.

# 3 Results and discussions

## 3.1 Accuracy of the PCM<sub>PP</sub> simplified dynamical scheme

We first compare in Fig. 2 the error associated with the simplified plasmon-pole GW@PCM<sub>PP</sub> approach where only the low-frequency reaction field is needed [see eq 17]. The error is calculated with respect to the fully dynamical GW@PCM(ω) calculations [see eq 16] that requires the recalculation of the reaction field at each needed frequency.

<figure><img src="image_2.png" /></figure>

Figure 2: Error (meV) on the GW@PCM quasiparticle energies associated with the simplified (plasmon-pole) dynamical PCM<sub>PP</sub> approach (eq 17 with Ω<sub>p</sub>=21 eV). The error is taken with respect to the fully dynamical PCM(ω) approach, for which the PCM reaction fields are recalculated at each needed frequency. Errors for occupied and virtual energy levels are represented with blue-filled and red-empty diamonds, respectively.

Spanning the occupied and virtual energy levels over a large energy window for the full set of molecules, we observe that decoupling the spatial and frequency degrees of freedom in the PCM<sub>PP</sub> approach leads to remarkably small errors of the order of a very few meVs. Such an error is seen to increase for virtual states with positive energy (above the vacuum level). However, the increase in error is shown in the SI (Fig. S3) to correlate with the spilling of the associated one-body
<table><caption>Table 1: HOMO and LUMO PCM polarization energies P@PCM(ω) as given by the fully dynamical GW@PCM(ω) approach. The polarization energies associated with the GW@PCM<sub>PP</sub> and GW@PCM<sub>COHSEX</sub> are given as differences (ΔP@PCM<sub>PP</sub> and ΔP@PCM<sub>COHSEX</sub>) with respect to the P@PCM(ω) reference. All energies are given in eV. Spilling values (spill.) express the amount of each state density leaking outside of the PCM cavity (in %). Mean values are provided for each quantity. We separate the case of water monomer and hexamers when taking the averages in relation with the large spilling of the associated LUMOs.</caption><thead><tr><th rowspan="2"></th><th colspan="4">HOMO</th><th colspan="4">LUMO</th></tr><tr><th>spill.</th><th>P@PCM(ω)</th><th>ΔP@PCM<sub>PP</sub></th><th>ΔP@PCM<sub>COHSEX</sub></th><th>spill.</th><th>P@PCM(ω)</th><th>ΔP@PCM<sub>PP</sub></th><th>ΔP@PCM<sub>COHSEX</sub></th></tr></thead><tbody><tr><td>Acetaldehyde</td><td>0.8</td><td>1.137</td><td>-0.0013</td><td>-0.059</td><td>1.6</td><td>-1.194</td><td>-0.0018</td><td>-0.053</td></tr><tr><td>Acroleine</td><td>0.7</td><td>1.058</td><td>-0.0012</td><td>-0.054</td><td>1.6</td><td>-1.035</td><td>-0.0022</td><td>-0.056</td></tr><tr><td>Adenine</td><td>0.6</td><td>0.741</td><td>-0.0018</td><td>-0.060</td><td>1.2</td><td>-1.061</td><td>-0.0017</td><td>-0.045</td></tr><tr><td>Benzene-TCNE</td><td>0.6</td><td>1.124</td><td>-0.0020</td><td>-0.054</td><td>0.5</td><td>-0.570</td><td>-0.0010</td><td>-0.028</td></tr><tr><td>Cytosine</td><td>0.7</td><td>0.839</td><td>-0.0019</td><td>-0.063</td><td>1.2</td><td>-0.875</td><td>-0.0018</td><td>-0.059</td></tr><tr><td>Ethanol</td><td>2.1</td><td>0.916</td><td>-0.0026</td><td>-0.087</td><td>2.2</td><td>-1.361</td><td>-0.0019</td><td>-0.057</td></tr><tr><td>Formaldehyde</td><td>0.8</td><td>1.344</td><td>-0.0014</td><td>-0.070</td><td>2.7</td><td>-1.205</td><td>-0.0024</td><td>-0.070</td></tr><tr><td>Indigo</td><td>0.6</td><td>0.677</td><td>-0.0017</td><td>-0.051</td><td>0.7</td><td>-0.840</td><td>-0.0014</td><td>-0.035</td></tr><tr><td>PNA<sub>perp</sub></td><td>0.7</td><td>1.400</td><td>-0.0017</td><td>-0.057</td><td>0.6</td><td>-1.104</td><td>-0.0009</td><td>-0.029</td></tr><tr><td>PNA<sub>plan</sub></td><td>0.6</td><td>1.120</td><td>-0.0019</td><td>-0.064</td><td>0.6</td><td>-1.159</td><td>-0.0008</td><td>-0.027</td></tr><tr><td>Probe</td><td>0.4</td><td>1.043</td><td>-0.0014</td><td>-0.052</td><td>0.6</td><td>-0.773</td><td>-0.0008</td><td>-0.027</td></tr><tr><td>Thymine</td><td>0.6</td><td>1.102</td><td>-0.0016</td><td>-0.055</td><td>1.1</td><td>-0.762</td><td>-0.0017</td><td>-0.042</td></tr><tr><td>Uracil</td><td>0.7</td><td>1.177</td><td>-0.0019</td><td>-0.061</td><td>1.2</td><td>-0.766</td><td>-0.0018</td><td>-0.045</td></tr><tr><td>mean val.</td><td>0.8</td><td>1.052</td><td>-0.0017</td><td>-0.061</td><td>1.2</td><td>-0.977</td><td>-0.0016</td><td>-0.044</td></tr><tr><td>Water</td><td>0.6</td><td>1.491</td><td>-0.0015</td><td>-0.074</td><td>29.2</td><td>-1.338</td><td>-0.0067</td><td>-0.211</td></tr><tr><td>(H<sub>2</sub>O)<sub>6</sub> book</td><td>0.5</td><td>0.787</td><td>-0.0012</td><td>-0.050</td><td>22.5</td><td>-1.187</td><td>-0.0082</td><td>-0.190</td></tr><tr><td>(H<sub>2</sub>O)<sub>6</sub> cage</td><td>0.5</td><td>0.916</td><td>-0.0011</td><td>-0.053</td><td>20.6</td><td>-1.106</td><td>-0.0076</td><td>-0.178</td></tr><tr><td>(H<sub>2</sub>O)<sub>6</sub> prism</td><td>0.6</td><td>0.854</td><td>-0.0013</td><td>-0.050</td><td>19.1</td><td>-1.062</td><td>-0.0077</td><td>-0.173</td></tr><tr><td>(H<sub>2</sub>O)<sub>6</sub> ring</td><td>0.5</td><td>1.079</td><td>-0.0013</td><td>-0.054</td><td>24.6</td><td>-1.445</td><td>-0.0086</td><td>-0.205</td></tr><tr><td>mean val.</td><td>0.5</td><td>1.025</td><td>-0.0013</td><td>-0.056</td><td>23.2</td><td>-1.228</td><td>-0.0078</td><td>-0.192</td></tr></tbody></table>
wavefunction outside the cavity, namely a situation where the PCM approximation becomes more
questionable.

The fully dynamical GW@PCM(ω) polarization energies for the HOMO and LUMO energies
are reported in Table 1. The polarization energy per level is rather stable around 1 eV in absolute
value. This leads to a gap closing of ~2 eV from the gas phase to the solvated phase. The errors
induced by the PCMpp scheme are shown to amount to about 1-2 meV, except for the water monomer
and hexamer LUMO levels where it increases to ~7-8 meV, in relation with the very large associated
spilling.

The accuracy of the GW@PCM<sub>PP</sub> approximation is very remarkable. This indicates that
standard PCM calculations, relying only on the knowledge of the low-frequency reaction field, can
be easily extended to include the effect of the dynamics of the solvent electronic degrees of freedom.
Furthermore, and as shown in the SI (Fig. S2), the polarization energy varies rather weakly with the
value of the pole energy, with a shift of the order of the meV for a shift of Ω<sub>p</sub> of the order of the eV.
This indicates that a very precise determination of Ω<sub>p</sub> is not a central issue, provided that the low
and high frequency limits of [1/ε<sub>opt</sub>(ω) - 1], or equivalently χ<sup>E</sup>(ω), are satisfied.

Along that line, we note that the water pole energy value (Ωₚ = 21 eV) is close to the classical
plasma frequency Ω<sub>plasmon</sub> = √(4πn) (a.u.) with n the electron density. For water, Ω<sub>plasmon</sub> amounts
indeed to 19.2 eV, accounting for the 8 valence electrons per water molecule, or 21.5 eV including
the oxygen 1s electrons. The rather large stability of the polarization energy with respect to the
plasmon frequency suggests that selecting the classical plasma frequency may be a simple and
accurate strategy in the case of solvents for which no sufficient experimental data for ε<sub>opt</sub>(ω) in the
full visible-UV range are available. Such a strategy can also be retrieved by plugging a one-pole
model into the f-sum rule.<sup>66</sup>

To provide some understanding of the accuracy of the single-plasmon-pole model for χ^E(r, r'; ω) (eq 17), we note that in the simple Born model of a unit charge at the center of a spherical cavity of
radius R, the polarization energy is directly related to $[1/\epsilon_{opt} - 1]$ through:

$$P = \frac{e^2}{2R} \left[ \frac{1}{\epsilon_{opt}} - 1 \right]$$

In the case of a non-spherical cavity, a similar relation holds for the integral of the induced surface charge on the cavity:

$$\int_{\Gamma} \sigma(\mathbf{x}) \, d\mathbf{x} = \left[ \frac{1}{\epsilon_{\text{opt}}} - 1 \right] Q$$

with Q the charge in the cavity inducing the distribution $\sigma(\mathbf{x})$ of surface charges in the simplest single-layer IEF approach to PCM.⁹,²³ These two observations shed light on the fact that the effect of the solvent on the monopole of the solute charge redistribution, associated with the photoemission charging process, is directly proportional to $f(\omega, \Omega_p)$. As such, factoring out the frequency dependence of the response in eq 17 mostly amounts to neglecting the form factor of the cavity and the higher order multipolar (dipole, etc.) contributions from the reaction field.

## 3.2 Accuracy of the PCM<sub>COHSEX</sub> scheme

Even though including the full dynamics of the solvent electrons can be performed at no cost as compared to standard PCM calculations in the BO limit, we now explore the accuracy of the GW@PCM<sub>COHSEX</sub> approach. The reference is again the fully dynamical GW@PCM(ω) calculation. The error is represented in Fig. 3 for occupied and unoccupied energy levels over a large energy window for the full set of solutes. The error is further provided in Table 1 for the HOMO and LUMO levels.

Analyzing Fig. 3, the error is found to be of the order of -50 meV for occupied states and unoccupied states below the vacuum level. With an absolute polarization energy of the order of the electronvolt (in absolute value), this represents a deviation of about 5%. This is significantly larger than the error associated with the dynamical PCM<sub>PP</sub> scheme. Whether such an error should be considered as moderate or large clearly depends on the required accuracy for the problem of interest. For states above the vacuum level, the error can be much larger, with again a clear correlation of the
<figure><img src="image_3.png" /></figure>

Figure 3: Error (meV) on the GW@PCM quasiparticle energies associated with the PCM<sub>COHSEX</sub> approach. The error is taken with respect to the fully dynamical PCM(ω) approach for which the PCM reaction fields are recalculated at each needed frequency. Errors for occupied and virtual energy levels are represented with blue-filled and red-empty diamonds, respectively.

error with the spilling outside the cavity (see SI Fig. S3). For such states, the PCM approximation becomes questionable independently of the treatment of environment electronic degrees of freedom dynamics.

The systematically negative error indicates that the BO approximation underestimates the polarization energy for the occupied levels but overestimates the correction for the virtual energy levels. In other words, the static approximation leads to ionization potential and electronic affinities in solution that are too large as compared to a fully dynamical PCM calculation. This leads to a partial cancellation of error when considering the polarization energy associated with the HOMO-LUMO gap, namely the closing in energy of the gap from the gas phase to the solvated phase.

It is interesting to analyze the decomposition reported in eq 23 partitioning the polarization energy in a SEX, COH and integral renormalization terms. This is represented in Fig. 4. For the occupied states, the renormalization $P_{Dyn}$ term is seen to be significantly smaller than the perturbative $P_{SEX}$ and $P_{COH}$ terms. In the limit of a reaction field very slowly varying over the solute, it was further shown<sup>92</sup> that the relation $P_{SEX} = -2P_{COH}$ should hold. This relation is only
qualitatively verified in the present case.

Concerning the unoccupied levels, the screened-exchange $P_{SEX}$ contribution is much smaller, as again justified in the limit of a smooth reaction field.⁵⁵,⁹² The renormalization $P_{Dyn}$ term is again significantly smaller than the perturbative contribution $P_{COH}$ that hardly changes from occupied to unoccupied energy levels. As such, the COH-like term does not contribute much to the renormalization of the HOMO-LUMO gap in solution as compared to the gas phase. Neglecting the COH contribution can thus be qualitatively correct when considering the evolution of energy differences between virtual and unoccupied levels, but fails for the absolute value of the ionization potential and electronic affinities.

<figure><img src="image_4.png" /><figcaption>Figure 4: Decomposition of the polarization energy for (a) occupied and (b) unoccupied energy levels. The $P_{Dyn}$, $P_{SEX}$, $P_{COH}$ correspond to the first, second and third-line contributions from eq 23, respectively. Dashed lines indicate occupied state manifold average values.</figcaption></figure>
To conclude this exploration of the solvent electron dynamics, we observe that the error as-
sociated with the PCM<sub>COHSEX</sub> scheme does not seem to correlate with the gap of the molecules
(given in the Table S1 of the SI). For sake of illustration, the errors associated, e.g., with the HOMO
polarization energies for indigo and benzene-TCNE, the two smallest gap molecules in the test set,
are not smaller than the ones associated with the water hexamers that display much larger gaps. For
the water solute in a water solvent, one may expect the BO approximation to become questionable
since there is no decoupling of energy between the solute and solvent electronic degrees of freedom.
However, clearly, the error does not seem to increase significantly.

Our results are in line with a previous QM/QM' exploration of the embedded GW scheme,⁵⁵
with an equivalent COHSEX-like treatment of the reaction field, showing that for fullerenes in a
fullerene environment, the approximation that the electrons in the environment react instantaneously
to an electronic excitation in the central subsystem would induce errors limited to a few percent
as well. Nevertheless, in the case of a metallic environment (a metallic nanotube), the error was
shown to increase to a few tenths of an eV for the polarization energy associated with the HOMO
and LUMO levels of a water molecule inside the tube. As such, the approximation of faster solvent
electrons seems rather robust, except in extreme situations where the solvent electronic degrees
of freedom are characterized by much smaller energies as compared to the solute. We emphasize
however that in any case, a fully dynamical model can be easily set up, with a cost similar to that
associated with a PCM approach where the reaction field is only required in the low-frequency
limit of the electronic dielectric response.

## 4 Conclusions

We have studied the impact on the electronic energy level polarization energies of considering
explicitly the frequency dependence of the optical dielectric constant of water as a solvent. By po-
larization energies, we mean the energy levels shifts from the gas phase to the solvated environment.
This study is performed at the GW@PCM level, with the GW operator depending explicitly on the
$\epsilon_{\text{opt}}^{-1}(\mathbf{r}, \mathbf{r}'; \omega)$ inverse dielectric function through the screened Coulomb potential $W$. Such an explicit frequency dependence offers an ideal playground to explore the impact of the dynamical nature of the solvent optical dielectric response, relying further on a recent parametrization of $\epsilon_{\text{opt}}(\omega)$ in the visible-UV range for water.⁷⁹ Accounting for the frequency dependence of the solvent optical dielectric response goes beyond the approximation restricting $\epsilon_{\text{opt}}(\omega)$ to its low frequency $\epsilon_\infty$ value. This restriction is expected to be valid only in the limit where the solvent electrons are assumed to be much faster than that of the solute. Such a limit is labeled the Born-Oppeinheimer (BO) limit in the PCM literature.

In the framework of GW@PCM calculations, the limit of instantaneous solvent electrons response can be properly obtained by assuming that the solvent susceptibility $\chi^E(\omega)$ present a pole structure with poles brought to infinity. This yields the so-called static Screened-Exchange and Coulomb-Hole contributions to the solvent reaction field. This static limit is compared to reference dynamical PCM($\omega$) calculations where the cavity reaction fields are recalculated explicitly for each needed frequency when building the GW self-energy integral. We observe that the BO limit leads to errors of the order of -50 meV for the polarization energy associated with individual energy levels. Due to cancellations of error between occupied and unoccupied levels polarization energies, the error on the HOMO-LUMO gap of the solute is found to be smaller. We emphasize that using directly the low frequency limit of the reaction field, without accounting for the pole structure of the solvent susceptibility, leads to neglecting the Coulomb-hole contribution, with erroneous polarization energies for individual energy levels.

Besides appraising the error associated with the instantaneous solvent electronic response limit, we have introduced a very simple and accurate approach to dynamical solvents within a single plasmon-pole model for the solvent electronic susceptibility. Such a formulation is approximate as it decouples real-space and frequency degrees of freedom but is strictly valid for the monopolar component of the reaction field. As a great advantage, only the low-frequency cavity reaction field associated with the solvent electronic degrees of freedom needs to be calculated. Finite-frequency reaction fields can further be obtained by simple rescaling. Remarkably, this approach, that we
label GW@PCM<sub>PP</sub> offers an accuracy at the few meV level as compared to the reference PCM(ω)
calculations.

The water single-pole frequency (Ω<sub>p</sub>=21 eV) was obtained from the best fit of the available
experimental [1/ε<sub>opt</sub> − 1] function. Such a value is close to the classical plasma frequency Ω<sub>plasmon</sub> =
√4πn (a.u.) associated with water. Further, the GW@PCM<sub>PP</sub> polarization energies are rather
insensitive to the exact Ω<sub>p</sub> value, with a shift of about a meV for an eV change of Ω<sub>p</sub>. This suggests
that in cases where the full visible/UV ε<sub>opt</sub>(ω) function is not known experimentally for a given
solvent, accurate fully dynamical calculations can still be straightforwardly conducted.

Finally, the merging of the GW formalism with the PCM opens the way to combining the Bethe-Salpeter formalism (BSE)<sup>93</sup> with the PCM for the study of optical excitations in systems immersed in a polarizable environment.<sup>89,94</sup> The BSE/PCM combination was shown in particular to account simultaneously for linear-response and state-specific contributions to the solvatochromic shifts.<sup>89</sup> This is an important feature for excitations presenting a hybrid local (Frenkel) and charge-transfer character. Even though the most common BSE implementations rely on the low frequency limit W(ω → 0) of the screened Coulomb potential, dynamical implementations are being explored.<sup>93,95-102</sup> The present study may facilitate the implementation of a dynamical BSE formalism combined with a dynamical PCM for the environment.

**Supporting Information Available**

Supporting Information: (I) a multiple-pole fit of ε⁻¹<sub>opt</sub>(ω) as an improvement over equation 16 and in
(II) the error for selected polarization energies associated with the single-pole PCM(ω) approach
with respect to the multiple-fit values, considering several single-pole Ω<sub>p</sub> energies. In (III), we
evidence the correlation between the errors associated with the PCM<sub>PP</sub> and PCM<sub>COHSEX</sub> schemes
and the spilling of the orbitals outside the cavity. In (IV), a Table compiling the Kohn-Sham and
various evGW@PCM gaps are provided for all molecules of the test set.
Acknowledgement

XB acknowledges numerous discussions with Gabriele D'Avino. DA is indebted to ENS Paris-Saclay for his PhD fellowship. This project was provided with computer and storage ressources by GENCI@TGCC thanks to the grants A0130910016 and A0150910016 on the Joliot-Curie supercomputer (SKL and Rome partitions). XB and ID acknowledge support from the French Agence Nationale de la Recherche (ANR) under contract ANR-20-CE29-0005.

A Appendix: Demonstration of equation 20

In order to tackle the limit $\Omega_p \to \infty$ within eq 20, we first remark that the integral along a closed contour $C$ enclosing the two poles $\Omega_1$ and $\Omega_2$ of the following product:

$$ \int_C d\omega \frac{1}{\omega - \Omega_1} \frac{1}{\omega - \Omega_2} = \frac{2i\pi}{\Omega_2 - \Omega_1} + \frac{2i\pi}{\Omega_1 - \Omega_2} = 0 \qquad (24) $$

is always null. We can use this property to rewrite the GW integral by considering only the residues taken at the poles of $G$. Once we have such an expression, it becomes possible to take the limit value of these residues when $\Omega_p \to \infty$. Finally, we can revert the expression by casting the sum over limit values of the residue as an integral of the form $G(E + \omega)[\tilde{W}(\omega) - \tilde{v}(0)]$.

To start our demonstration, we first rewrite the term $[W - \tilde{v}]$ as $\Delta^{(+)} + \Delta^{(-)}$, where $\Delta^{(+)}$ regroup all the poles of $W - \tilde{v}$ within the complex upper-plane while $\Delta^{(-)}$ holds all the poles of the complex lower-plane. Similarly, we split $G$ into $G_{occ}$ and $G_{vir}$ with poles occupying also respectively the
upper-plane and the lower-plane.

$$ 
\begin{aligned}
& \frac{i}{2\pi} \int_{-\infty}^{\infty} d\omega \, e^{i\omega\eta} G(E + \omega) [W(\omega) - \tilde{v}(\omega)] \\
& = \frac{i}{2\pi} \int_{-\infty}^{\infty} d\omega \, e^{i\omega\eta} \left( G_{\text{occ}}(E + \omega) + G_{\text{vir}}(E + \omega) \right) \\
& \qquad \times \left( \Delta^{(+)}(\omega) + \Delta^{(-)}(\omega) \right).
\end{aligned} 
\qquad (25) $$

We insist on the fact that all these sub-quantities have a $1/\omega$ asymptotic behavior, with no constant terms, as they are only composed of simple poles. Thus, expression 25 is clearly an integral over products of simple poles. Using eq 24 as selection rule for which of these product residues to keep after integration, we can now express the equation 20 as:

$$ 
\begin{aligned}
& \frac{i}{2\pi} \int_{-\infty}^{\infty} d\omega \, e^{i\omega\eta} G(E + \omega) [W(\omega) - \tilde{v}(\omega)] \\
& = -\sum_{i} \text{Res}(G_{\text{occ}}, \varepsilon_i) \Delta^{(-)}(\varepsilon_i - E) \\
& \qquad -\sum_{k} G_{\text{vir}}(E + \Omega_k) \text{Res}(\Delta^{(+)}, \Omega_k).
\end{aligned} 
\qquad (26) $$

$\text{Res}(f, x)$ corresponds to the residue of the function $f$, associated to its pole $x$. We can then use eq 24 once again to rewrite the sum over $k$ as a sum over the empty states $a$ as:

$$ 
\begin{aligned}
& \sum_{k} G_{\text{vir}}(E + \Omega_k) \text{Res}(\Delta^{(+)}, \Omega_k) \\
& = -\sum_{a} \text{Res}(G_{\text{vir}}, \varepsilon_a) \Delta^{(+)}(\varepsilon_a - E),
\end{aligned} 
\qquad (27) $$
leading to an expression of eq 20 that involves only residues taken at the poles of G:

$$ \begin{aligned} & \frac{i}{2\pi} \int_{-\infty}^{\infty} d\omega \, e^{i\omega\eta} G(E + \omega) [W(\omega) - \tilde{v}(\omega)] \\ &= -\sum_{i} \mathrm{Res}(G_{\mathrm{occ}}, \varepsilon_i) \Delta^{(-)}(\varepsilon_i - E) \\ &\quad + \sum_{a} \mathrm{Res}(G_{\mathrm{vir}}, \varepsilon_a) \Delta^{(+)}(\varepsilon_a - E). \end{aligned} \quad (28) $$

Now that we have an expression that requires only evaluation of $\Delta^{(+)}$ and $\Delta^{(-)}$ at finite value of $\omega$, we can explore their limit when $\Omega_p \to \infty$. In order to do so, we remark that for any finite value $\omega$, we can write:

$$ \begin{aligned} \lim_{\Omega_p \to \infty} W(\omega) - \tilde{v}(\omega) &= \lim_{\Omega_p \to \infty} \tilde{v}(\omega) \cdot \tilde{\chi}(\omega) \cdot \tilde{v}(\omega) \\ &= \tilde{v}(0) \cdot \tilde{\chi}(\omega) \cdot \tilde{v}(0), \end{aligned} \quad (29) $$

with $\tilde{\chi}(\omega)$ also taken at the limit:

$$ \tilde{\chi}(\omega) = \chi_0^S(\omega) + \chi_0^S(\omega) \cdot \tilde{v}(0) \cdot \tilde{\chi}(\omega). \quad (30) $$

Thus, we can write:

$$ \lim_{\Omega_p \to \infty} W(\omega) - \tilde{v}(\omega) = \tilde{W}(\omega) - \tilde{v}(0), \quad (31) $$

with

$$ \tilde{W}(\omega) = \tilde{v}(0) + \tilde{v}(0) \cdot \chi_0^S(\omega) \cdot \tilde{W}(\omega). \quad (32) $$

The screened Coulomb integrals $\tilde{W}$ correspond thus to the screened PCM Coulomb integrals obtained within the static limit of the solvent response, and renormalized by the dynamical solute free susceptibility. We can then proceed as previously by separating $\tilde{W}(\omega) - \tilde{v}(0) = \tilde{\Delta}^{(+)}(\omega) + \tilde{\Delta}^{(-)}(\omega)$
into contributions of the upper and lower-plane poles. The equation 31 being true for all finite $\omega$, this implies that we now can make the identification:

$$ \lim_{\Omega_p \to \infty} \Delta^{(+/-)}(\omega) = \tilde{\Delta}^{(+/-)}(\omega). \qquad (33) $$

The limit values of $\Delta^{(+/-)}(\omega)$ can be re-injected into eq 28:

$$ \begin{aligned} & \lim_{\Omega_p \to \infty} \frac{i}{2\pi} \int_{-\infty}^{\infty} d\omega \, e^{i\omega\eta} G(E + \omega) [W(\omega) - \tilde{v}(\omega)] \\ & = -\sum_i \mathrm{Res}(G_{\mathrm{occ}}, \epsilon_i) \tilde{\Delta}^{(-)}(\epsilon_i - E) \\ & \quad + \sum_a \mathrm{Res}(G_{\mathrm{vir}}, \epsilon_a) \tilde{\Delta}^{(+)}(\epsilon_a - E). \end{aligned} \qquad (34) $$

It is now possible to revert the derivation starting from eq 34 and rewinding the steps 28 back to 25 to finally obtain that:

$$ \begin{aligned} & \lim_{\Omega_p \to \infty} \frac{i}{2\pi} \int_{-\infty}^{\infty} d\omega \, e^{i\omega\eta} G(E + \omega) [W(\omega) - \tilde{v}(\omega)] \\ & = \frac{i}{2\pi} \int_{-\infty}^{\infty} d\omega \, e^{i\omega\eta} G(E + \omega) [\tilde{W}(\omega) - \tilde{v}(0)]. \end{aligned} \qquad (35) $$

## References

(1) Warshel, A.; Levitt, M. Theoretical studies of enzymic reactions: Dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme. *J. Mol. Bio.* **1976**, *103*, 227–249.

(2) Miertuš, S.; Scrocco, E.; Tomasi, J. Electrostatic interaction of a solute with a continuum. A direct utilization of *ab initio* molecular potentials for the prevision of solvent effects. *Chem. Phys.* **1981**, *55*, 117–129.
(3) Born, M. Volumen und Hydratationswärme der Ionen. Z. Phys. **1920**, 1, 45.

(4) Jackson, J. *Classical Electrodynamics*; Wiley; 2nd edition, 1975; Chapter Multipoles, Electrostatics of Macroscopic Media, Dielectrics.

(5) Klamt, A.; Schüürmann, G. COSMO: a new approach to dielectric screening in solvents with explicit expressions for the screening energy and its gradient. *J. Chem. Soc., Perkin Trans. 2* **1993**, 799–805.

(6) Cammi, R.; Tomasi, J. Analytical derivatives for molecular solutes. I. Hartree-Fock energy first derivatives with respect to external parameters in the polarizable continuum model. *J. Chem. Phys.* **1994**, *100*, 7495–7502.

(7) Cammi, R.; Tomasi, J. Nonequilibrium solvation theory for the polarizable continuum model: A new formulation at the SCF level with application to the case of the frequency-dependent linear electric response function. *Int. J. Quant. Chem.* **1995**, *56*, 465–474.

(8) Tomasi, J.; Persico, M. Molecular Interactions in Solution: An Overview of Methods Based on Continuous Distributions of the Solvent. *Chem. Rev.* **1994**, *94*, 2027–2094.

(9) Cancès, E.; Mennucci, B.; Tomasi, J. A new integral equation formalism for the polarizable continuum model: Theoretical background and applications to isotropic and anisotropic dielectrics. *J. Chem. Phys.* **1997**, *107*, 3032–3041.

(10) Barone, V.; Cossi, M. Quantum Calculation of Molecular Energies and Energy Gradients in Solution by a Conductor Solvent Model. *J. Phys. Chem. A* **1998**, *102*, 1995–2001.

(11) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation Models. *Chem. Rev.* **2005**, *105*, 2999–3094.

(12) Caricato, M.; Ingrosso, F.; Mennucci, B.; Tomasi, J. A time-dependent polarizable continuum model: Theory and application. *J. Chem. Phys.* **2005**, *122*, 154501.
(13) Nguyen, P. D.; Ding, F.; Fischer, S. A.; Liang, W.; Li, X. Solvated First-Principles Excited-State Charge-Transfer Dynamics with Time-Dependent Polarizable Continuum Model and Solvent Dielectric Relaxation. *J. Phys. Chem. Lett* **2012**, *3*, 2898–2904.

(14) Liang, W.; Chapman, C. T.; Ding, F.; Li, X. Modeling Ultrafast Solvated Electronic Dynamics Using Time-Dependent Density Functional Theory and Polarizable Continuum Model. *J. Phys. Chem. A* **2012**, *116*, 1884–1890.

(15) Corni, S.; Pipolo, S.; Cammi, R. Equation of Motion for the Solvent Polarization Apparent Charges in the Polarizable Continuum Model: Application to Real-Time TDDFT. *J. Phys. Chem. A* **2015**, *119*, 5405–5416.

(16) Ding, F.; Lingerfelt, D. B.; Mennucci, B.; Li, X. Time-dependent non-equilibrium dielectric response in QM/continuum approaches. *J. Chem. Phys.* **2015**, *142*, 034120.

(17) Gil, G.; Pipolo, S.; Delgado, A.; Rozzi, C. A.; Corni, S. Nonequilibrium Solvent Polarization Effects in Real-Time Electronic Dynamics of Solute Molecules Subject to Time-Dependent Electric Fields: A New Feature of the Polarizable Continuum Model. *J. Chem. Theory Comput.* **2019**, *15*, 2306–2319.

(18) Dall'Osto, G.; Gil, G.; Pipolo, S.; Corni, S. Real-time dynamics of plasmonic resonances in nanoparticles described by a boundary element method with generic dielectric function. *J. Chem. Phys.* **2020**, *153*, 184114.

(19) Guido, C. A.; Rosa, M.; Cammi, R.; Corni, S. An open quantum system theory for polarizable continuum models. *J. Chem. Phys.* **2020**, *152*, 174114.

(20) Kim, H. J.; Hynes, J. T. Equilibrium and nonequilibrium solvation and solute electronic structure. III. Quantum theory. *J. Chem. Phys.* **1992**, *96*, 5088–5110.

(21) Kuznetsov, A. M. Role of high-frequency and low-frequency polarization of the medium in
the kinetics of electron transfer and thermodynamics of solvation. *J. Chem. Phys.* **1992**, *96*, 3337–3345.

(22) Marcus, R. A. Schroedinger equation for strongly interacting electron-transfer systems. *J. Phys. Chem.* **1992**, *96*, 1753–1757.

(23) Duchemin, I.; Jacquemin, D.; Blase, X. Combining the GW formalism with the polarizable continuum model: A state-specific non-equilibrium approach. *J. Chem. Phys.* **2016**, *144*, 164106.

(24) Kim, S.-J.; Lebègue, S.; Ringe, S.; Kim, H. GW Quasiparticle Energies and Bandgaps of Two-Dimensional Materials Immersed in Water. *J. Phys. Chem. Lett.* **2022**, *13*, 7574–7582.

(25) Clary, J. M.; Del Ben, M.; Sundararaman, R.; Vigil-Fowler, D. Impact of solvation on the GW quasiparticle spectra of molecules. *J. App. Phys.* **2023**, *134*, 085001.

(26) Hedin, L. New Method for Calculating the One-Particle Green's Function with Application to the Electron-Gas Problem. *Phys. Rev.* **1965**, *139*, A796–A823.

(27) Rostgaard, C.; Jacobsen, K. W.; Thygesen, K. S. Fully self-consistent GW calculations for molecules. *Phys. Rev. B* **2010**, *81*, 085103.

(28) Blase, X.; Attaccalite, C.; Olevano, V. First-principles GW calculations for fullerenes, porphyrins, phtalocyanine, and other molecules of interest for organic photovoltaic applications. *Phys. Rev. B* **2011**, *83*, 115103.

(29) Marom, N.; Caruso, F.; Ren, X.; Hofmann, O. T.; Körzdörfer, T.; Chelikowsky, J. R.; Rubio, A.; Scheffler, M.; Rinke, P. Benchmark of GW methods for azabenzenes. *Phys. Rev. B* **2012**, *86*, 245127.

(30) Bruneval, F.; Marques, M. A. L. Benchmarking the Starting Points of the GW Approximation for Molecules. *J. Chem. Theory Comput.* **2013**, *9*, 324–329.
(31) van Setten, M. J.; Caruso, F.; Sharifzadeh, S.; Ren, X.; Scheffler, M.; Liu, F.; Lischner, J.; Lin, L.; Deslippe, J. R.; Louie, S. G.; Yang, C.; Weigend, F.; Neaton, J. B.; Evers, F.; Rinke, P. GW100: Benchmarking G0W0 for Molecular Systems. *J. Chem. Theory Comput.* **2015**, *11*, 5665–5687.

(32) Katharina Krause, M. E. H.; Klopper, W. Coupled-cluster reference values for the GW27 and GW100 test sets for the assessment of GW methods. *Mol. Phys.* **2015**, *113*, 1952–1960.

(33) Knight, J. W.; Wang, X.; Gallandi, L.; Dolgounitcheva, O.; Ren, X.; Ortiz, J. V.; Rinke, P.; Körzdörfer, T.; Marom, N. Accurate Ionization Potentials and Electron Affinities of Acceptor Molecules III: A Benchmark of GW Methods. *J. Chem. Theory Comput.* **2016**, *12*, 615–626.

(34) Kaplan, F.; Harding, M. E.; Seiler, C.; Weigend, F.; Evers, F.; van Setten, M. J. Quasi-Particle Self-Consistent GW for Molecules. *J. Chem. Theory Comput.* **2016**, *12*, 2528–2541.

(35) Rangel, T.; Hamed, S. M.; Bruneval, F.; Neaton, J. B. Evaluating the GW Approximation with CCSD(T) for Charged Excitations Across the Oligoacenes. *J. Chem. Theory Comput.* **2016**, *12*, 2834–2842.

(36) Caruso, F.; Dauth, M.; van Setten, M. J.; Rinke, P. Benchmark of GW Approaches for the GW100 Test Set. *J. Chem. Theory Comput.* **2016**, *12*, 5076–5087.

(37) Maggio, E.; Liu, P.; van Setten, M. J.; Kresse, G. GW100: A Plane Wave Perspective for Small Molecules. *J. Chem. Theory Comput.* **2017**, *13*, 635–648.

(38) Govoni, M.; Galli, G. GW100: Comparison of Methods and Accuracy of Results Obtained with the WEST Code. *J. Chem. Theory Comput.* **2018**, *14*, 1895–1909.

(39) Förster, A.; Visscher, L. GW100: A Slater-Type Orbital Perspective. *J. Chem. Theory Comput.* **2021**, *17*, 5080–5097.

(40) Rojas, H. N.; Godby, R. W.; Needs, R. J. Space-Time Method for Ab Initio Calculations
of Self-Energies and Dielectric Response Functions of Solids. *Phys. Rev. Lett.* **1995**, *74*, 1827–1830.

(41) Foerster, D.; Koval, P.; Sánchez-Portal, D. An $O(N^3)$ Implementation of Hedins GW Approximation for Molecules. *J. Chem. Phys* **2011**, *135*, 074105.

(42) Neuhauser, D.; Gao, Y.; Arntsen, C.; Karshenas, C.; Rabani, E.; Baer, R. Breaking the Theoretical Scaling Limit for Predicting Quasiparticle Energies: The Stochastic GW Approach. *Phys. Rev. Lett.* **2014**, *113*, 076402.

(43) Liu, P.; Kaltak, M.; Klimeš, J.; Kresse, G. Cubic scaling GW: Towards fast quasiparticle calculations. *Phys. Rev. B* **2016**, *94*, 165109.

(44) Vlček, V.; Rabani, E.; Neuhauser, D.; Baer, R. Stochastic GW Calculations for Molecules. *J. Chem. Theory Comput.* **2017**, *13*, 4997–5003.

(45) Vlček, V.; Li, W.; Baer, R.; Rabani, E.; Neuhauser, D. Swift GW beyond 10,000 electrons using sparse stochastic compression. *Phys. Rev. B* **2018**, *98*, 075107.

(46) Wilhelm, J.; Golze, D.; Talirz, L.; Hutter, J.; Pignedoli, C. A. Toward GW Calculations on Thousands of Atoms. *J. Phys. Chem. Lett.* **2018**, *9*, 306–312.

(47) Förster, A.; Visscher, L. Low-Order Scaling $G_0W_0$ by Pair Atomic Density Fitting. *J. Chem. Theory Comput.* **2020**, *16*, 7381–7399.

(48) Kim, M.; Martyna, G. J.; Ismail-Beigi, S. Complex-time shredded propagator method for large-scale GW calculations. *Phys. Rev. B* **2020**, *101*, 035139.

(49) Kutepov, A. Self-consistent GW method: $O(N)$ algorithm for polarizability and self energy. *Comput. Phys. Commun.* **2020**, *257*, 107502.

(50) Gao, W.; Chelikowsky, J. R. Accelerating Time-Dependent Density Functional Theory and GW Calculations for Molecules and Nanoclusters with Symmetry Adapted Interpolative Separable Density Fitting. *J. Chem. Theory Comput.* **2020**, *16*, 2216–2223.
(51) Wilhelm, J.; Seewald, P.; Golze, D. Low-Scaling GW with Benchmark Accuracy and Application to Phosphorene Nanosheets. *J. Chem. Theory Comput.* **2021**, *17*, 1662–1677.

(52) Duchemin, I.; Blase, X. Cubic-Scaling All-Electron GW Calculations with a Separable Density-Fitting Space–Time Approach. *J. Chem. Theory Comput.* **2021**, *17*, 2383–2393.

(53) Förster, A.; van Lenthe, E.; Spadetto, E.; Visscher, L. Two-Component GW Calculations: Cubic Scaling Implementation and Comparison of Vertex-Corrected and Partially Self-Consistent GW Variants. *J. Chem. Theory Comput.* **2023**, *19*, 5958–5976.

(54) Scott, C. J. C.; Backhouse, O. J.; Booth, G. H. A “moment-conserving” reformulation of GW theory. *J. Chem. Phys.* **2023**, *158*, 124102.

(55) Amblard, D.; Blase, X.; Duchemin, I. Static versus dynamically polarizable environments within the many-body GW formalism. *J. Chem. Phys.* **2024**, *160*, 154104.

(56) In Ref. 55, the assumption of an instantaneous response from the environment electronic degrees of freedom was called an adiabatic approximation. While the association of the Born-Oppenheimer (BO) limit with the wording adiabatic was used in Ref. 20, the wording antiadiabatic was used in Ref. 103 for the same regime of faster solvent electronic response. Since the adiabatic limit for the solvent electrons can be seen as an antiadiabatic limit for the solute electrons, we stick here to the BO wording.

(57) Baumeier, B.; Rohlfing, M.; Andrienko, D. Electronic Excitations in Push-Pull Oligomers and Their Complexes with Fullerene from Many-Body Green's Functions Theory with Polarizable Embedding. *J. Chem. Theory Comput.* **2014**, *10*, 3104–3110.

(58) Li, J.; D'Avino, G.; Duchemin, I.; Beljonne, D.; Blase, X. Combining the Many-Body GW Formalism with Classical Polarizable Models: Insights on the Electronic Structure of Molecular Solids. *J. Phys. Chem. Lett.* **2016**, *2814–2820*.
(59) Li, J.; D'Avino, G.; Duchemin, I.; Beljonne, D.; Blase, X. Accurate description of charged excitations in molecular solids from embedded many-body perturbation theory. *Phys. Rev. B* **2018**, *97*, 035108.

(60) Wehner, J.; Brombacher, L.; Brown, J.; Junghans, C.; Çaylak, O.; Khalak, Y.; Madhikar, P.; Tirimbò, G.; Baumeier, B. Electronic Excitations in Complex Molecular Environments: Many-Body Green's Functions Theory in VOTCA-XTP. *J. Chem. Theory Comput.* **2018**, *14*, 6253–6268.

(61) Tirimbò, G.; Sundaram, V.; Çaylak, O.; Scharpach, W.; Sijen, J.; Junghans, C.; Brown, J.; Ruiz, F. Z.; Renaud, N.; Wehner, J.; Baumeier, B. Excited-state electronic structure of molecules using many-body Green's functions: Quasiparticles and electron–hole excitations with VOTCA-XTP. *J. Chem. Phys.* **2020**, *152*, 114103.

(62) Thole, B. Molecular polarizabilities calculated with a modified dipole interaction. *Chem. Phys.* **1981**, *59*, 341–350.

(63) D'Avino, G.; Muccioli, L.; Castet, F.; Poelking, C.; Andrienko, D.; Soos, Z. G.; Cornil, J.; Beljonne, D. Electrostatic phenomena in organic semiconductors: fundamentals and implications for photovoltaics. *J. Phys. Condens. Matter* **2016**, *28*, 433002.

(64) Bondanza, M.; Nottoli, M.; Cupellini, L.; Lipparini, F.; Mennucci, B. Polarizable embedding QM/MM: the future gold standard for complex (bio)systems? *Phys. Chem. Chem. Phys.* **2020**, *22*, 14433–14448.

(65) Strinati, G.; Mattausch, H. J.; Hanke, W. Dynamical Correlation Effects on the Quasiparticle Bloch States of a Covalent Crystal. *Phys. Rev. Lett.* **1980**, *45*, 290–294.

(66) Hybertsen, M. S.; Louie, S. G. Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies. *Phys. Rev. B* **1986**, *34*, 5390–5413.
(67) Godby, R. W.; Schlüter, M.; Sham, L. J. Self-energy operators and exchange-correlation potentials in semiconductors. *Phys. Rev. B* **1988**, *37*, 10159–10175.

(68) von der Linden, W.; Horsch, P. Precise quasiparticle energies and Hartree-Fock bands of semiconductors and insulators. *Phys. Rev. B* **1988**, *37*, 8351–8362.

(69) Farid, B.; Daling, R.; Lenstra, D.; van Haeringen, W. GW approach to the calculation of electron self-energies in semiconductors. *Phys. Rev. B* **1988**, *38*, 7530–7534.

(70) Aryasetiawan, F.; Gunnarsson, O. The GW method. *Rep. Prog. Phys.* **1998**, *61*, 237–312.

(71) Farid, B. In *Electron Correlation in the Solid State - Chapter 3*; March, N., Ed.; Imperial College Press, London, 1999.

(72) Onida, G.; Reining, L.; Rubio, A. Electronic excitations: density-functional versus many-body Green's-function approaches. *Rev. Mod. Phys.* **2002**, *74*, 601–659.

(73) Ping, Y.; Rocca, D.; Galli, G. Electronic excitations in light absorbers for photoelectrochemical energy conversion: first principles calculations based on many body perturbation theory. *Chem. Soc. Rev.* **2013**, *42*, 2437–2469.

(74) Leng, X.; Jin, F.; Wei, M.; Ma, Y. GW method and Bethe–Salpeter equation for calculating electronic excitations. *WIREs Comput. Mol. Sci.* **2016**, *6*, 532–550.

(75) Martin, R.; Reining, L.; Ceperley, D. *Interacting Electrons: Theory and Computational Approaches*; Cambridge University Press, 2016.

(76) Golze, D.; Dvorak, M.; Rinke, P. The GW Compendium: A Practical Guide to Theoretical Photoemission Spectroscopy. *Front. Chem.* **2019**, *7*, 377.

(77) It may seem somehow awkward to use the notation $\epsilon_\infty$ for the low frequency limit of the macroscopic dielectric constant in the optical range. However, in the PCM language, the frequencies associated with the electronic degrees of freedom remain always much larger than that associated with the slow degrees of freedom contributing to $\epsilon_0$.
(78) Duchemin, I.; Blase, X. Robust Analytic-Continuation Approach to Many-Body GW Calculations. *J. Chem. Theory Comput.* **2020**, *16*, 1742–1756.

(79) Fiedler, J.; Boström, M.; Persson, C.; Brevik, I.; Corkery, R.; Buhmann, S. Y.; Parsons, D. F. Full-Spectrum High-Resolution Modeling of the Dielectric Function of Water. *J. Phys. Chem. B* **2020**, *124*, 3103–3113.

(80) Tyagi, R.; Voora, V. K. Single-Pole Polarization Models: Rapid Evaluation of Electron Affinities of Solvated-Electron and Superatomic Molecular Anionic States. *J. Phys. Chem. Lett.* **2024**, *15*, 1218–1226.

(81) Vahtras, O.; Almlöf, J.; Feyereisen, M. Integral approximations for LCAO-SCF calculations. *Chem. Phys. Lett.* **1993**, *213*, 514–518.

(82) Ren, X.; Rinke, P.; Blum, V.; Wieferink, J.; Tkatchenko, A.; Sanfilippo, A.; Reuter, K.; Scheffler, M. Resolution-of-identity approach to Hartree-Fock, hybrid density functionals, RPA, MP2 and GW with numeric atom-centered orbital basis functions. *New J. Phys.* **2012**, *14*, 053020.

(83) Duchemin, I.; Li, J.; Blase, X. Hybrid and Constrained Resolution-of-Identity Techniques for Coulomb Integrals. *J. Chem. Theory Comput.* **2017**, *13*, 1199–1208.

(84) Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. *Phys. Chem. Chem. Phys.* **2005**, *7*, 3297–3305.

(85) Weigend, F.; Häser, M.; Patzelt, H.; Ahlrichs, R. RI-MP2: optimized auxiliary basis sets and demonstration of efficiency. *Chem. Phys. Lett.* **1998**, *294*, 143–152.

(86) Perdew, J. P.; Ernzerhof, M.; Burke, K. Rationale for mixing exact exchange with density functional approximations. *J. Chem. Phys.* **1996**, *105*, 9982–9985.
(87) Adamo, C.; Barone, V. Toward reliable density functional methods without adjustable parameters: The PBE0 model. *J. Chem. Phys.* **1999**, *110*, 6158–6170.

(88) Neese, F. Software update: The ORCA program system—Version 5.0. *WIREs Comput. Mol. Sci* **2022**, *12*, e1606.

(89) Duchemin, I.; Guido, C. A.; Jacquemin, D.; Blase, X. The Bethe–Salpeter formalism with polarisable continuum embedding: reconciling linear-response and state-specific features. *Chem. Sci.* **2018**, *9*, 4430–4443.

(90) Segarra-Martí, J.; Merchán, M.; Roca-Sanjuán, D. Ab initio determination of the ionization potentials of water clusters ($H_2O)_n$ ($n$ = 2-6). *J. Chem. Phys.* **2012**, *136*, 244306.

(91) Blase, X.; Boulanger, P.; Bruneval, F.; Fernandez-Serra, M.; Duchemin, I. GW and Bethe-Salpeter study of small water clusters. *J. Chem. Phys.* **2016**, *144*, 034109.

(92) Neaton, J. B.; Hybertsen, M. S.; Louie, S. G. Renormalization of Molecular Electronic Levels at Metal-Molecule Interfaces. *Phys. Rev. Lett.* **2006**, *97*, 216405.

(93) Strinati, G. Application of the Green's functions method to the study of the optical properties of semiconductors. *La Rivista del Nuovo Cimento* **1988**, *11*, 1–86.

(94) Kim, S.-J.; Lebègue, S.; Ringe, S.; Kim, H. Elucidating Solvatochromic Shifts in Two-Dimensional Photocatalysts by Solving the Bethe–Salpeter Equation Coupled with Implicit Solvation Method. *J. Phys. Chem. Lett.* **2024**, *15*, 4575–4580.

(95) Strinati, G. Dynamical Shift and Broadening of Core Excitons in Semiconductors. *Phys. Rev. Lett.* **1982**, *49*, 1519–1522.

(96) Ma, Y.; Rohlfing, M.; Molteni, C. Excited states of biological chromophores studied using many-body perturbation theory: Effects of resonant-antiresonant coupling and dynamical screening. *Phys. Rev. B* **2009**, *80*, 241405.
(97) Loos, P.-F.; Blase, X. Dynamical correction to the Bethe–Salpeter equation beyond the plasmon-pole approximation. *J. Chem. Phys.* **2020**, *153*, 114120.

(98) Bintrim, S. J.; Berkelbach, T. C. Full-frequency dynamical Bethe–Salpeter equation without frequency and a study of double excitations. *J. Chem. Phys.* **2022**, *156*, 044114.

(99) Yamada, S.; Noguchi, Y.; Ishii, K.; Hirose, D.; Sugino, O.; Ohno, K. Development of the Bethe-Salpeter method considering second-order corrections for a GW electron-hole interaction kernel. *Phys. Rev. B* **2022**, *106*, 045113.

(100) Loos, P.-F.; Romaniello, P. Static and dynamic Bethe–Salpeter equations in the T-matrix approximation. *J. Chem. Phys.* **2022**, *156*, 164101.

(101) Zhang, X.; Leveillee, J. A.; Schleife, A. Effect of dynamical screening in the Bethe-Salpeter framework: Excitons in crystalline naphthalene. *Phys. Rev. B* **2023**, *107*, 235205.

(102) Monino, E.; Loos, P.-F. Connections and performances of Green's function methods for charged and neutral excitations. *J. Chem. Phys.* **2023**, *159*, 034105.

(103) Phan Huu, D. K. A.; Dhali, R.; Pieroni, C.; Di Maiolo, F.; Sissa, C.; Terenziani, F.; Painelli, A. Antiadiabatic View of Fast Environmental Effects on Optical Spectra. *Phys. Rev. Lett.* **2020**, *124*, 107401.
TOC Graphic

<figure><img src="image_5.png" /></figure>