DOI:
10.1039/C8SC00529J
(Edge Article)
Chem. Sci., 2018,
9, 44304443
The Bethe–Salpeter formalism with polarisable continuum embedding: reconciling linearresponse and statespecific features†
Received
1st February 2018
, Accepted 2nd April 2018
First published on 5th April 2018
The Bethe–Salpeter equation (BSE) formalism has been recently shown to be a valuable alternative to timedependent density functional theory (TDDFT) with the same computing time scaling with system size. In particular, problematic transitions for TDDFT such as chargetransfer, Rydberg and cyaninelike excitations were shown to be accurately described with BSE. We demonstrate here that combining the BSE formalism with the polarisable continuum model (PCM) allows us to include simultaneously linearresponse and statespecific contributions to solvatochromism. This is confirmed by exploring transitions of various natures (local, chargetransfer, etc.) in a series of solvated molecules (acrolein, indigo, pnitroaniline, donor–acceptor complexes, etc.) for which we compare BSE solvatochromic shifts to those obtained by linearresponse and statespecific TDDFT implementations. Such a remarkable and unique feature is particularly valuable for the study of solvent effects on excitations presenting a hybrid localised/chargetransfer character.
1 Introduction
The exploration of the excitedstate (ES) properties of chemical systems certainly stands as a central question in theoretical chemistry. Indeed, ES phenomena govern many applications such as solar energy conversion, photocatalysis, lightemission or optical information storage. Further, while experimental characterisations can provide reference absorption and/or emission spectra, they are less suited to obtain some key information, e.g., ES geometries, nature of the excitation (localised, chargetransfer, etc.) or time evolution of hot electrons. Such a need for quantum mechanical formalisms allowing us to study realistic systems certainly explains the formidable popularity of timedependent density functional theory (TDDFT)^{1,2} that can be used to study the optical properties of systems comprising up to a few hundred atoms, thanks to a (formal) scaling with system size. Further, the availability of analytical TDDFT derivatives^{3–7} together with the extension of efficient continuum models, such as the polarisable continuum model (PCM),^{8,9} to TDDFT^{10–15} has dramatically helped in bridging the gap between quantum simulations and realistic systems, by respectively allowing us to explore ES potential energy surfaces and to take into account the impact of the surroundings. In TDDFT, the coupling with the PCM was initially performed within a linearresponse (LR) formalism,^{10,11} that is, using the electronic transition density for including solvent effects. While such a LR model is generally accurate for describing the local ES, it is less suited for the chargetransfer (CT) ES, in which a large reorganisation of the electron density occurs. To tackle such an ES, statespecific (SS) PCMTDDFT models, in which the solvatochromic effects depend on the total electronic density of the ES, have been developed.^{7,13–15} At this stage, let us point out that, in using these SS models, one can encounter some cases for which the exact details of the chosen SS approach as well as the selected exchangecorrelation functional have a very large impact on the results, especially when selfconsistent iterative methods are selected.^{7,16–18} In addition, as first pointed out by Corni et al.,^{19} who used a simple formal model explicitly including two states for the solute and two solvent macrostates, there is a need to simultaneously account for both LR and SS effects. However, to date, only the ad hoc sum of both LR and SS terms, determined in the context of a corrected linear response (cLR) approach,^{13} was proposed in a TDDFT context.^{20} Alternatively one can turn towards singlereference electroncorrelated wavefunction approaches, such as ADC(2), CC2, CCSD or SACCI, that have all been coupled to continuum models,^{21–30} but these models imply a significantly increased computational effort compared to TDDFT. In this framework, we underline that the importance of the inclusion of both LR and SS effects was also clearly underlined by Lunkenheimer and Köhn in their work describing the coupling of the ADC(2) theory to a continuum approach of solvation effects.^{26}
As another alternative to TDDFT, the Bethe–Salpeter equation (BSE) formalism^{31–35} has been recently experiencing a growing interest in the study of molecular systems due to its ability to overcome some of the problems that TDDFT is facing, including chargetransfer^{36–44} and cyaninelike^{45,46} excitations, while preserving the same scaling in its standard implementations. Extensive benchmark studies on diverse molecular families have been performed,^{47–53} demonstrating that excellent agreement with higherlevel manybody wavefunction techniques, such as coupledcluster (CC3) or CASPT2, could be obtained for all types of transitions, provided that they do not present a strong multipleexcitation character. Singlet–triplet transitions constitute the only notable exception as they may present the same instability problems with BSE and TDDFT.^{51,52} We have recently reviewed the differences between BSE and TDDFT formalisms in a chemical context, and we refer the interested reader to that original work for more details.^{35}
As compared to TDDFT, the BSE formalism relies on transition matrix elements between occupied and virtual energy levels calculated at the GW level, where G and W stand for the onebody Green's function and the screenedCoulomb potential. These GW energy levels, including HOMO and LUMO frontier orbital energies, were shown to be in much better agreement with reference wavefunction calculations as compared to standard DFT Kohn–Sham (KS) eigenvalues. Namely, GW HOMO and LUMO energies can be directly associated with the ionisation potential (IP) and electronic affinity (AE). Further, the TDDFT exchangecorrelation kernel matrix elements in the occupiedtovirtual transition space are replaced within the BSE formalism by matrix elements involving the screened Coulomb potential interaction between the hole and the electron.
The coupling of the GW and PCM formalisms was recently achieved,^{54} within an integral equation formalism^{9} (IEFPCM) implementation including the socalled nonequilibrium (neq) effects related to the separation of the solvent response into its “fast” electronic and “slow” nuclear contributions. The electron and hole polarisation energies, namely the shifts of the ionisation potential and electronic affinity from the gas phase to solution, were very well reproduced with GW taking as a reference standard ΔSCF + PCM calculations where total energy calculations of the solvated ions and neutral species were performed at the DFT or CCSD levels. Similar studies were also performed adopting discrete polarisable models to study organic semiconductors and complexes of interest for optoelectronic applications.^{55–58} As central (GW + PCM) formalism features, we underline that the polarisation energies for all occupied (P_{i}^{+}) and virtual (P_{a}^{−}) energy levels can be calculated, not only those associated with frontier orbitals, and that these polarisation energies are SS. Such polarisation energies, of the order of an eV in the case of water solvated nucleobases,^{54} dramatically reduce the HOMO–LUMO gap as compared to its gas phase value (see Fig. 1).

 Fig. 1 Schematic representation of the statespecific contribution to solvatochromic shifts within the (BSE/GW + PCM) formalism. The GW energy levels are first renormalised by the SS solventinduced polarisation energies P_{a}^{−}(ε_{∞}) and P_{i}^{+}(ε_{∞}) for electrons (occupied) and holes (virtual), respectively. Further, the strength of the screened electron–hole interaction W_{eh} is also renormalised, namely reduced by a quantity that we label here P^{eh}(ε_{∞}), thanks to the additional screening provided by the PCM. Groundstate polarisation effects associated with the slow (ε_{0}) degrees of freedom are incorporated at the initial (DFT + PCM) level (not represented here).  
In the present study, we demonstrate that the BSE formalism combined with the PCM in a nonequilibrium formulation intrinsically combines the LR and SS solvatochromic contributions associated with the effect of the polarisable environment. Such a remarkable feature hinges in particular on the proper inclusion of dynamic polarisation energies not only for the occupied and virtual energy levels, calculated within the GW formalism, but also for the screened Coulomb electron–hole interaction. This renormalisation by polarisation of energy levels and electron–hole interactions leads to the SS contribution to solvatochromic shifts, in addition to the familiar transition matrix elements stemming from the LR contributions. This simultaneous account of both LR and SS effects was, to the best of our knowledge, only achieved using an a posteriori sum of the two terms in the context of the computationally more expensive ADC(2)/COSMO approach,^{26} and more recently with a TDDFT/discrete polarizable scheme.^{20} Here, the obtained (BSE + PCM) solvatochromic shifts are computed for paradigmatic transitions in acrolein, indigo, pnitroaniline (PNA), a small donor–acceptor complex and a solvatochromic probe (see Fig. 2), and are compared to the sum of the shifts obtained at the TDDFT level within the standard LR and SS implementations of the PCM, respectively.

 Fig. 2 Schematic representation of the molecules and complexes studied: (a) acrolein (A), (b) indigo (I), (c) paranitroaniline (PNA), (d) “twisted” paranitroaniline (PNA_{perp}), (e) the donor–acceptor benzene–tetracyanoethylene (BTCNE) complex and (f) the 4nitropyridine Noxide solvatochromic probe.  
2 Formalism
We will not detail here the BSE formalism^{31–34,59–61} for which several reviews are available.^{35,62,63} We restrict ourselves to highlight the important points that allow understanding the merging with the PCM formalism. Since the BSE approach requires as an input the screened Coulomb potential W and accurate occupied/unoccupied {ε^{GW}_{i/a}} energy levels calculated at the GW level, we start by a short introduction to the Green's function GW formalism^{62,64–67} from which the BSE approach derives.
2.1
GW formalism
The GW approach belongs to the family of manybody perturbation theories (MBPT) that takes as central quantity the Green's function (G) instead of, e.g., the charge density within DFT. In its timeordered formulation, the Green's function G reads: 
 (1) 
where {ε_{i},ε_{a}} are the “true” electronic energy levels as defined experimentally in a direct/inverse photoemission experiment. The small positive infinitesimal η controls the proper analytic properties of G in the complex energy plane.^{67} More precisely, we define: ε_{a} = E(N + 1, a) − E(N, 0) for “virtual” energy levels and ε_{i} = E(N, 0) − E(N − 1, 0) for “occupied” energy levels, where E(N + 1, a) is the total energy of the (N + 1)electron system in its a^{th} quantum state, E(N − 1, i) the total energy of the (N − 1)electron system in its i^{th} quantum state, and E(N, 0) the Nelectron system groundstate energy. {ϕ_{i}, ϕ_{a}} are called “Lehman amplitudes”. It can be formally shown that G verifies a simple Dyson equation: 
 (2) 
with, e.g., 1 = (r_{1},t_{1}) a spacetime coordinate. G^{0} is the independentelectron Green's function and the “selfenergy” operator Σ^{HXC} contains all electron–electron interactions (Hartree plus exchange and correlation). Plugging the expression for G into the Dyson equation leads to a familiar eigenvalue equation for the (photoemission) excitation energies defined here above, namely: 
 (3) 
where the h^{^}_{0} Hamiltonian contains the kinetic, ionic and classical Hartree operators. Such an equation is formally equivalent to the DFT KS equation, but with an exchangecorrelation “potential” that is both nonlocal and energydependent.
While eqn (3) is exact, an expression for Σ^{XC} should be defined. The GW formalism provides an approximation for Σ^{XC} to first order in the screened Coulomb potential W, with:

 (4) 

 (5) 

 (6) 

 (7) 
where
v(
r,
r′) is the bare Coulomb potential,
χ_{0} the independentelectron susceptibility and
W the screened Coulomb potential.
E_{F} is the energy of the Fermi level. {
f_{i/j}} are occupation numbers and the input {
ϕ_{n},
ε_{n}} eigenstates are typically KS eigenstates that will be corrected within the
GW formalism.
The GW eigenvalues have been shown in many benchmark studies on gas phase molecular systems to be significantly more accurate than KS or Hartree–Fock eigenstates, providing, e.g., frontier orbital energies within a few tenths of an eV with respect to reference CCSD(T) calculations.^{68–72} For the sake of illustration, our gas phase ionization potential (IP) and electronic affinity (EA) are 10.35 eV (−ε^{GW}_{HOMO}) and 0.68 eV (+ε^{GW}_{LUMO}) respectively for acrolein within our GW approach, comparable to 10.01 eV and 0.70 eV within CCSD(T) for the same atomic basis (ccpVTZ) and the same geometry (see the ESI†). In practice, GW calculations proceed traditionally by correcting input KS eigenstates, substituting the GW selfenergy contribution to the DFT exchangecorrelation potential:

 (8) 
The obtained ε^{GW}_{a/i} eigenvalues and the screened Coulomb potential W serve as input quantities for the BSE excitation energy calculation.
2.2 Bethe–Salpeter equation
While TDDFT starts from the evaluation of the density–density susceptibility χ(12) = ∂n(1)/∂U^{ext}(2) that measures the variation of the chargedensity with respect to an external local perturbation, the excitation energies within BSE are obtained through the poles of a generalized susceptibility L(1234) = ∂G(12)/∂U^{ext}(34), where U^{ext}(34) is a nonlocal external perturbation. Deriving the Dyson equation for G (eqn (2)) leads in particular to the introduction of the (∂Σ^{XC}/∂G) derivative, the analog to the exchangecorrelation kernel within DFT. In the standard BSE/GW formalism, we consistently assume the GW approximation for Σ^{XC}. Expressing then the 4point susceptibility L in the transition space between occupied and unoccupied onebody eigenstates yields a standard linear algebra eigenvalue representation: 
 (9) 
where Ω^{BSE}_{λ} represents the BSE excitation energies. X_{λ} represents the components of the twobody electron–hole ψ_{λ}(r_{e},r_{h}) eigenstates over the ϕ_{i}(r_{h})ϕ_{a}(r_{e}) transition basis, where (i,a) index (occupied, virtual) eigenstates, while Y_{λ} gathers the ϕ_{i}(r_{e}) ϕ_{a}(r_{h}) deexcitation components. Such a linear algebra representation in the transition space resembles the standard TDDFT formalism within the socalled Casida's formulation^{1} but with modified matrix elements. For the resonant block, one obtains: 
A^{BSE}_{ai,bj} = δ_{ab}δ_{ij}(ε^{GW}_{a} − ε^{GW}_{i}) − 〈abWij〉 + 〈aibj〉  (10) 
with 
〈abWij〉 = 〈ϕ_{a}(r)ϕ_{b}(r)W(r,r′)ϕ_{i}(r′)ϕ_{j}(r′)〉  (11) 

〈aibj〉 = 〈ϕ_{a}(r)ϕ_{i}(r)v(r,r′)ϕ_{b}(r′)ϕ_{j}(r′)〉  (12) 
where ε^{GW}_{a/i} are the GW unoccupied/occupied energy levels and W the screened Coulomb potential. The (ε^{GW}_{a} − ε^{GW}_{i}) energy differences replace the TDDFT KS (ε^{KS}_{a} − ε^{KS}_{i}) energy differences between virtual and occupied states, while the W matrix elements can be interpreted as electron–hole interaction terms through the screened Coulomb potential. {ϕ_{i/a}} are the starting KS onebody eigenstates that are not corrected in the GW implementation selected here (see below).
2.3 Merging with the PCM formalism
The merging of the GW formalism with the PCM is described in detail in ref. 54. The seminal point is that the solvent reaction field can be straightforwardly incorporated into the screened Coulomb potential W that accounts for both the polarisability of the solute and that of the solvent when combined with the PCM. Under the assumption that the solute and solvent eigenstates do not spatially overlap, it can be shown that the screened Coulomb potential W is changed from its gas phase value: 
 (13) 
to its condensed phase expression: 
ṽ = v + vχ^{PCM}_{0}ṽ  (15) 
where we have dropped the frequency and space variables in the second equation for conciseness. χ^{QM}_{0} is the gas phase independentelectron susceptibility of the quantum subsystem to be solvated, while χ^{PCM}_{0} is that of the PCM solvent. The modified ṽ potential is thus the bare Coulomb potential renormalized by the solvent reaction field, equivalent to:where χ^{PCM} is the full (interacting) susceptibility of the solvent. The quantity [vχ^{PCM}v] (r_{1},r_{2}) represents the reaction field (v^{reac}) generated in r_{2} by the PCM surface charges developed in response to a unity pointcharge added in r_{1}, with (r_{1},r_{2}) located in the cavity carved in the solvent to accommodate the solute. Details of our implementation using localised Gaussian bases and the Coulombfitting resolutionofidentity (RIV) approach can be found in ref. 54 for the merging of GW with the PCM and ref. 56 and 58 for the merging with a discrete polarisable model.
Importantly, in the present nonequilibrium formulation of the (BSE/GW + PCM) implementation, the reaction field that modifies W is associated with the “fast” electronic excitations, namely with the ε_{∞} dielectric constant (e.g., ε_{∞} = 1.78 in water) that is the lowfrequency optical dielectric response. The inclusion of the “slow” response of the solute (given by, e.g., ε_{0} = 78.35 for water) is accounted for in the preliminary groundstate DFT + PCM(ε_{0}) run that serves as a starting point for GW and BSE calculations, namely in the construction of χ_{0} and G in eqn (5) and (7). The flow of calculations can then be summarised as follows:
(1) Calculation of input KS {ε_{n}, ϕ_{n}} eigenstates within a DFT + PCM(ε_{0}) scheme. These eigenstates contain groundstate solvation effects,
(2) Calculation of the “fast” reaction field v^{reac}(ε_{∞}) in an auxiliary basis representation (see ref. 54, 56 and 58) that is incorporated inside the screened Coulomb potential W following eqn (5) and (7),
(3) Correction of the {ε_{n}} KS eigenstates, that include “slow” polarisation effects, with the GW selfenergy operator that contains the “fast” v^{reac}(ε_{∞}) reaction field in order to yield the GW + PCM neq {ε_{n}^{GW+PCM}} energy levels,
(4) Resolution of the BSE excitation energy eigenvalue problem with the {ε_{n}^{GW+PCM}} eigenstates and the W and ṽ potentials that include v^{reac}(ε_{∞}) = vχ^{PCM}(ε_{∞})v.
2.4 LR and SS solvatochromic contributions
The renormalisation of v into v + v^{reac}(ε_{∞}) and the resulting change ΔW, accounting for switching on the fast ε_{∞} solvent response, leads to modifying eqn (10) as follows: 
 (17a) 

+δ_{ab}δ_{ij}Δ(ε^{GW/PCM}_{a} − ε^{GW/PCM}_{i}) − 〈abΔWij〉  (17b) 

+〈aivχ^{PCM}(ε_{∞})vbj〉,  (17c) 
where A^{BSE/PCM0} is the BSE Hamiltonian that includes the groundstate charges only, namely built with ε_{∞} = 1. Correspondingly, Δε^{GW/PCM}_{a/i} corresponds to the shift of the GW energy levels when switching on the fast solvent response, that is, . Such a decomposition allows a direct correspondence with the solvatochromic shifts calculated within TDDFT in both the LR and SS formulations.
As the easiest identification, the 3^{rd} line contribution (eqn (17c)) straightforwardly corresponds to LR reaction field matrix elements, i.e., to transition density polarisation effects. More explicitly, 〈aivχ^{PCM}vbj〉 describes the action on ϕ_{b}(r′)ϕ_{j}(r′) of the reaction field , where is the surface charge generated by the PCM susceptibility χ^{PCM}(ε_{∞}) in response to the field generated by the transition density ϕ_{a}(r)ϕ_{i}(r). In the original notations of ref. 10, such matrix elements are strictly equivalent to the B^{f}_{ai,bj} = 〈aiK^{f}bj〉 linear response terms (eqn (31) in ref. 10) where K^{f} is the “fast” reaction potential integral operator. As such, the BSE + PCM formalism includes the LR solvent contributions.
Let us now demonstrate that the second line (eqn (17b)) recovers the SS contribution. Part of the demonstration relies on the specificity of the GW + PCM formalism that captures accurately SS dynamic polarisation energies, namely:

 (18a) 

 (18b) 
with positive (negative) shifts for occupied (virtual) levels since Δ
W corresponds to a reduced interaction. This was demonstrated numerically in
ref. 54 and justified in Appendix A. As such, the term of
eqn (17b), that we label Δ
A^{SS}_{ai,bj}, reads:

 (19) 
In the case of a “pure” transition between levels (i) and (a), one straightforwardly obtains:

 (20) 
where (
aa −
ii) represents the unrelaxed variation Δ
ρ_{i→a} of the charge density between the excited and the ground states. Identifying Δ
W to the reaction field
vχ^{PCM}(
ε_{∞})
v as justified in Appendix B, we recognise the action on Δ
ρ_{i→a} of the reaction field
, where
Q_{i→a} (
r) is the surface charge generated by the PCM susceptibility
χ^{PCM}(
ε_{∞}) in response to the field generated by Δ
ρ_{i→a} itself, with the (1/2) factor indicating that it is a selfinteraction term. This corresponds to the SS expression within TDDFT for the solvatochromic shift, namely with the notations of
ref. 12:

 (21) 
where here
i represents the
i^{th} ES. Physically, the terms in
eqn (17b) represent the quantum mechanical version of the classical solvation energy associated with the interaction of a charge redistribution (dipolar, quadrupolar,
etc.) with the polarisable medium, where this charge redistribution corresponds to the charge variation between the ground and the excited states.
We now turn to the case of general BSE electron–hole eigenstates assuming for simplicity the TDA approximation. Hole and electron densities are now described by a correlated codensity ρ_{λ}(r_{e},r_{h}) = ψ_{λ}(r_{e},r_{h})^{2} that cannot be expressed in terms of individual eigenstate densities as in the simple twolevel model. This however does not affect the central result that the screened Coulomb potential W, and thus the electron–hole interaction:

 (22) 
are properly renormalized by the reaction field. Further, in the presence of the PCM, the BSE electron–hole Hamiltonian of
eqn (10) is first dressed with the reaction field operator
v^{reac} and then fully rediagonalised. As such, beyond firstorder perturbation theory, the (BSE + PCM) eigenstates and codensities are relaxed with respect to the presence of the polarisable medium and the variation from gas to solvent of the electron–hole interaction results from the variation of both
ρ_{λ} and
W. It should be emphasized however that in the present implementation, the molecular orbitals and corresponding codensities are not relaxed with respect to the change in QM charge density upon excitation. Namely, the socalled
Zmatrix relaxation effects are here neglected. From this point of view, the present BSE/PCM methodology is at that stage in analogy with the unrelaxed density (UD) approximation used in the Vertical Excitation Method (VEM) approach (VEMUD).
^{7,15}
Concerning now the contribution from the GW electron and hole quasiparticle energies, we obtain, considering e.g. the unoccupied (electron) energy levels ε^{GW/PCM}_{a}:

 (23) 
with
. The variation of the electron and hole energies from the gas to the solvent thus stems from the variation of the
X_{ia} coefficients (relaxation effects) and from that of the individual
ε^{GW/PCM}_{i/a} energy levels, namely their related polarization energy Δ
ε^{GW}_{i/a}. In the correlated electron–hole BSE description, the variation of the electronic energy levels contributes to the excitation energy solvatochromic shift through a weighted average of the individual level polarization energy.
3 Computational details
The geometries of all compounds have been obtained at the MP2/6311G(2d,2p) level in the gasphase, but for the benzene/TCNE complex the geometry was taken from ref. 73 (see the ESI†). Our GW and BSE calculations are performed with the Fiesta package^{44,47,53} that implements these formalisms with Gaussian basis sets and the Coulombfitting RIV approach.^{74} Input KS eigenstates are generated with the NWChem package^{75} at the PBE0/ccpVTZ level.^{76–78} The corresponding ccpVTZRI auxiliary basis set^{79} is used in the manybody RIV calculations. For improved accuracy and removal of most of the dependency on the selected DFT functional, our GW calculations are performed at the partially selfconsistent evGW level^{68,70,72,80,81} with GWupdated eigenvalues and frozen eigenvectors.
Within the present neq approach, the starting DFT calculations are performed in combination with the COSMO formalism^{82} as implemented in NWChem, using the (ε_{0}) dielectric constant associated with the (slow) nuclear and electronic degrees of freedom (e.g., ε_{0} = 78.39 for water). As such, the KS eigenstates used to build our electronic excitations carry information regarding the groundstate polarisation effects. In a secondstep, the screened Coulomb potential used in the GW and BSE calculations, that describes fast electronic excitations out of the groundstate, is “dressed” with the reaction field generated with the (fast) electronic dielectric response in the lowfrequency limit (e.g., ε_{∞} = 1.78 for water). In GW and BSE, the reaction field is described at the full IEFPCM level, following the implementation detailed in ref. 54.
Our BSE calculations are compared to TDDFT calculations performed with the Gaussian16 package^{83} using the IEFPCM nonequilibrium implementation and the same ccpVTZ atomic basis set. The LR shifts have been obtained with the default Gaussian16 implementation,^{10,11} whereas the SS shifts have been determined with the socalled corrected LR (cLR) formalism,^{13} that is a perturbative approach. For the TDDFT calculations we selected the same PBE0 functional, complemented in the case of CT transitions with calculations performed with the rangeseparated hybrid CAMB3LYP functional,^{84} known to more accurately describe CT transitions in the TDDFT context. The character of the electronic transitions in TDDFT has been determined by estimating the effective electronic displacement induced by the excitation using the Γ index.^{85,86}
In order to differentiate static and dynamic contributions within our approach, we make use of a ground state frozen polarization excitation energy (Ω_{0}),^{15,19} obtained by considering that the polarisable medium does not respond to the fast electronic excitation, namely by setting ε_{∞} = 1 while keeping the correct static dielectric constant of the solvent ε_{0}. As such, labelling Ω the final BSE excitation energies, accounting for both static and dynamic PCM responses, the quantity (Ω − Ω_{0}) quantifies the impact of switching on the fast PCM(ε_{∞}) response. We underline that our BSE calculations are performed beyond the Tamm–Dancoff approximation (TDA), that is, include the full BSE matrix. However, the inclusion of contributions from deexcitation processes complicates the simple analysis provided above concerning the LR and SS contributions with a clear distinction between 〈abWij〉 terms, namely the coupling between density terms, and 〈aiv^{reac}bj〉 contributions, that is the coupling between transition dipoles. Beyond TDA, we can decompose the expectation value 〈Ψ_{BSE}H_{BSE}Ψ_{BSE}〉 into

 (24a) 

 (24b) 

 (24c) 
where we have dropped the
λ excitation index for simplicity and where the factor 2 in front of the bare Coulomb potential indicates singlet transition energies. In such a decomposition, all terms involving occupiedvirtual orbital products, including the (boxed) 〈
aj
W
bi〉 term, contribute to the LR response, whereas the remaining matrix elements contribute to the SS response.
4 Results
To illustrate the methodology and to confirm that the PCMBSE formalism provides both LR and SS shifts, we perform calculations on a series of standard test molecules (Fig. 2), providing comparison with TDDFT/PCM calculations conducted with the standard LR and SS (cLR) response formalisms. We start with transitions that have a rather local character and hence are expected to be better described within TDDFT by the LR formalism (Subsection 4.1), before focusing in a second Subsection 4.2 on transitions having a strong CT character that requires the SS formalism. Eventually, in Subsection 4.3, we evidence that some systems may require both LR and SS contributions, highlighting the merits of the (BSE + PCM) formalism that can treat all systems on an equal footing.
4.1 Local ES: dominating LR contributions
We start our analysis with the wellknown cases of acrolein and indigo (Fig. 2a and b). Acrolein is characterised by a negative solvatochromism for the lowest n–π* (A′′) transition but a positive solvatochromism of the lowest π–π* (A′) excitation. Experimental data are available in both the gas phase and water,^{87,88} and acrolein has been used as a benchmark system for studying continuum solvation models in conjunction with several levels of theory, e.g., TDDFT,^{12,13,15} ADC(2),^{26} SACCI,^{21,23} and CCSD.^{22,25,27} Hybrid twolayer QM/MM^{88,89} or threelayer QM/MM/PCM^{90} calculations, combining TDDFT with an atomistic force field description of the explicit solvent, have also been used to understand the pros and cons of the PCM model for water, a protic medium wellknown to be challenging for continuum approaches.
Our data are compiled in Table 1. The gas phase (Ω_{g}) and solvated (Ω) theoretical and experimental excitation energies are provided, together with the overall shift (Ω − Ω_{g}). In our nonequilibrium formalism, the shift is decomposed into a groundstate static contribution (Ω_{0} − Ω_{g}) and a dynamic one (Ω − Ω_{0}) that is itself partitioned into LR and SS contributions. Within TDDFT, the LR and SS (cLR) shifts are obtained as two separate calculations, whereas within BSE, this decomposition is obtained by partitioning (see above) the nonequilibrium BSE/GW + PCM shift. Wavefunction approaches, such as ADC(2) or CCSD, can also be used to provide both contributions simultaneously^{26} or separately^{25} and we give some literature examples in Table 1.
Table 1 Lowest singlet excitation energies for acrolein and indigo in the gas phase (Ω_{g}) and in water (Ω), combining TDDFT or BSE with neq PCM (ε_{0} = 78.39, ε_{∞} = 1.78 for water). The Ω_{0} energies are obtained by setting ε_{∞} = 1, namely accounting only for groundstate polarisation effects with no additional surface charges induced by the excitation. The TDPBE0 (Ω − Ω_{0}) shifts are provided using either the LR or SS (cLR) formalisms. The (Ω − Ω_{0}) BSE shifts are also partitioned into LR and SS contributions following the analysis in Section 3. All values are in eV

Ω
_{g}

Ω(Ω − Ω_{g}) 
Ω
_{0} − Ω_{g} 
Ω − Ω_{0} 
Ref. 
LR 
SS 
In the breakdown approach used in that work, the SS contribution is +0.17 eV and the LR contribution is negligible.
In the breakdown approach used in that work, the SS contribution is −0.21 eV and the LR contribution is −0.18 eV.
In ethanol, the most polar protic solvent in which indigo is soluble experimentally.

Acrolein n–π* in water

TDPBE0 (LR) 
3.599 
3.785(+0.186) 
+0.189 
−0.003 

This work 
TDPBE0 (cLR) 
3.599 
3.736(+0.137) 
+0.189 

−0.052 
This work 
BSE 
3.736 
3.988(+0.252) 
+0.232 
−0.011 
+0.031 
This work 
CC3 
3.74 




This work 
CCSDR(3)/MM 
3.81 
4.08(+0.27) 



91

SACCI 
3.85 
3.95(+0.10) 



21

CCSD 
3.94 
4.14(+0.20) 



22

ADC(2) 
3.69 
3.86(+0.17)^{a} 



26

CCSD (LR) 
3.88 
4.10(+0.22) 



25

CCSD (SS) 
3.88 
4.05(+0.17) 



25

Exp. 
3.69 
3.94(+0.25) 



87 and 88


Acrolein π–π* in water

TDPBE0 (LR) 
6.383 
6.174(−0.209) 
−0.073 
−0.136 

This work 
TDPBE0 (cLR) 
6.383 
6.281(−0.102) 
−0.073 

−0.029 
This work 
BSE 
6.498 
6.214(−0.284) 
−0.112 
−0.163 
−0.004 
This work 
CC3 
6.82 




This work 
CCSDR(3)/MM 
6.73 
6.22(−0.51) 



91

SACCI 
6.97 
6.75(−0.22) 



21

CCSD 
6.89 
6.54(−0.35) 



22

ADC(2) 
6.79 
6.40(−0.39)^{b} 



26

CCSD (LR) 
6.80 
6.39(−0.41) 



25

CCSD (SS) 
6.80 
6.54(−0.26) 



25

Exp. 
6.42 
5.89(−0.53)^{a} 



87 and 88


Indigo in water

TDPBE0 (LR) 
2.304 
2.160(−0.144) 
−0.068 
−0.076 

This work 
TDPBE0 (cLR) 
2.304 
2.229(−0.075) 
−0.068 

−0.007 
This work 
BSE 
2.259 
2.047(−0.212) 
−0.122 
−0.082 
−0.008 
This work 
Exp. 
2.32 
2.04(−0.19)^{c} 



92, 93, 94, 95, 96, 97 and 98

The acrolein n–π* transition in the gas phase (Ω_{g}) is found to be located at 3.60 eV and 3.74 eV within TDPBE0 and BSE respectively, in good agreement with the CC3 value of 3.74 eV, as well as with previous wavefunction estimates and experiment. The analysis of the intermediate (Ω_{0} − Ω_{g}) and total (Ω − Ω_{g}) solvatochromic shift indicates that the positive solvatochromism for this transition is entirely dominated by groundstate effects and that the additional shift associated with the fast optical excitation is negligible. TDPBE0 and BSE calculations performed on top of the DFT + PCM(ε_{0}) groundstate yield similar Ω_{0} − Ω_{g} shifts, namely +0.19 eV and +0.23 eV, respectively. Concerning the effect of switching ε_{∞} (1.78 in water), both TDPBE0+PCM, within LR or cLR, and BSE + PCM yield very small additional shifts, ranging from −0.05 eV (cLR) to +0.02 eV (BSE). The BSE computed shift of +0.25 eV turns out to be in close agreement with the experimental values (+0.25 eV), as well as with previous wavefunction approaches (see Table 1).
While the n–π* transition does not allow us to clearly discriminate between LR and SS responses, a more interesting test of the effect of the fast response (ε_{∞}) polarisable environment comes with the higherlying π–π* transition. For this transition, the reaction field associated with the optical excitation, namely v^{reac}(ε_{∞}), leads to a shift that is larger than the one associated with groundstate solvation charges. In the gas phase, the BSE transition energy (6.50 eV) is reasonably close to the CC3 (6.82 eV) and experimental (6.42 eV) values. The most salient feature is that within TDDFT, only the LR scheme can significantly contribute to the redshift, while the cLR approach fails to deliver any sizeable solvation effect, with (Ω − Ω_{0}) being equal to −0.14 eV for the former model, and −0.03 eV for the latter. This effect was expected for a local π–π* transition not involving a strong density reorganisation between the two states: the LRPCMTDDFT is more suited as it captures the dominating contributions originating from the transition densities, that can be viewed as “dispersionlike” terms.^{19} At the CCSD level, the results of Caricato^{25} also demonstrated that the LR contribution is dominant; in that work the decomposition of the total response into various contributions was not performed, but rather two different models have been applied as in TDDFT. We also underline that in their ADC(2) study, Lunkenheimer and Köhn also found that the LR term, negligible for the n–π* case,^{26} becomes large for the π–π* transition, though the approach used to compute the relative contributions is not straightforwardly comparable to ours. In any case, the BSE + PCM formalism, with a −0.17 eV (Ω − Ω_{0}) shift, captures the correct bathochromic effect. Further, consistent with the TDDFT calculations, we observe that the BSE shift is dominated by the LR contribution, while the SS term provides a negligible shift. This shows, consistent with our analysis of eqn (17c), that the BSE formalism correctly captures the LR response. When compared to experiment, the BSE + PCM shift remains too small, but we recall that we neglect here, as in any continuum approach, the explicit solvent–solute interactions that are known to be significant in the present case.^{26,91}
To confirm the present observations, we consider the case of the lowest transition in indigo, a hallmark centrosymmetric dye presenting a lowlying dipoleallowed π–π* transition. This compound was studied previously at the TDDFT level with the LR PCM model, and it was shown that this approach nicely reproduces the experimental solvatochromic shifts.^{99} With TDDFT and the LR formalism, we found that the groundstate and dynamic polarisation effects, as measured respectively by (Ω_{0} − Ω_{g}) and (Ω − Ω_{0}), have the same order of magnitude, whereas the cLR correction does not lead to any significant (Ω − Ω_{0}) effect, as expected for a dye in which both the groundstate and excitedstate total dipoles are strictly null. As in the case of the π–π* transition in acrolein, the BSE + PCM scheme leads to a clear redshift, with in this case a good agreement with the experimental trends as well. As a matter of fact, the (Ω − Ω_{0}) BSE LR (SS) shift is very close to the corresponding TDDFT LR (SS) shift, demonstrating the relevance of the analysis and partitioning of the BSE overall (Ω − Ω_{0}) difference discussed in Section 3.
4.2 Chargetransfer ES: dominating SS contributions
We now turn to the study of chargetransfer (CT) excitations for which the SS contribution should increase with increasing CT character. We start with the pnitroaniline (PNA) molecule (Fig. 2c), a typical push–pull system studied both theoretically^{5,100} and experimentally,^{101} and characterised by a lowlying CT excitation from the amine (–NH_{2}) to the nitro (–NO_{2}) group. Inspired by ref. 5, we first consider an artificial configuration obtained by rotating the amino group perpendicularly to the conjugated plane (PNA_{perp} in Fig. 2d) so as to produce a transition having a pure dark CT character. We focus on the lowlying dark CT excitation showing a dominant HOMO–LUMO character. The “standard” planar PNA molecule is studied below. Further, we consider an intermolecular CT excitation in the donor–acceptor benzene/TCNE complex (Fig. 2e),^{38,73} for which gas phase experimental data are available.^{102} The lowest energy excitation is a clear HOMO–LUMO transition with a strong CT character even though a small delocalization of the HOMO (LUMO) on TCNE (Benzene) leads to a nonzero oscillator strength (see e.g.ref. 38). Since CT transitions are poorly described by global hybrid functionals with a moderate amount of exact exchange such as PBE0, TDDFT calculations are also performed with the CAMB3LYP rangeseparated hybrid functional. Our data are compiled in Table 2.
Table 2 Data as in Table 1 but for the lowest CT transitions in PNA_{perp} and in the benzene–TCNE complex. The experimental value for the latter case is taken from ref. 102

Ω
_{g}

Ω(Ω − Ω_{g}) 
Ω
_{0} − Ω_{g} 
Ω − Ω_{0} 
LR 
SS 
PNA
_{
perp
}
in water (“dark” CT excitation)

TDPBE0 (LR) 
3.686 
3.316(−0.370) 
−0.370 
0.000 

TDPBE0 (cLR) 
3.686 
2.861(−0.795) 
−0.370 

−0.425 
TDCAMB3LYP (LR) 
4.621 
4.308(−0.313) 
−0.312 
−0.001 

TDCAMB3LYP (cLR) 
4.621 
4.038(−0.583) 
−0.312 

−0.271 
BSE 
5.112 
4.399(−0.713) 
−0.423 
+0.015 
−0.304 

Benzene–TCNE in water (“bright” CT excitation)

TDPBE0 (LR) 
2.157 
2.081(−0.076) 
−0.065 
−0.011 

TDPBE0 (cLR) 
2.157 
1.747(−0.410) 
−0.065 

−0.345 
TDCAMB3LYP (LR) 
2.944 
2.876(−0.068) 
−0.061 
−0.007 

TDCAMB3LYP (cLR) 
2.944 
2.492(−0.452) 
−0.061 

−0.291 
BSE 
3.503 
3.121(−0.382) 
−0.086 
−0.009 
−0.287 
Exp. 
3.59 




While the contribution of the fast optical dielectric response (ε_{∞}) to the solvatochromic shifts, as measured by (Ω − Ω_{0}), mainly originates from the LR contribution in both acrolein and indigo, the solventinduced dynamic shift associated with CT transitions can only be described by adopting a SS (cLR here) formalism in the TDDFT context. This is clearly illustrated by the PNA_{perp} system where the TDDFT LR shift is trifling, a logical consequence of the dark nature of the considered transition, whereas the SS contribution is very large, as a result of the large density reorganisation associated with excitation. Consistent with previous studies,^{17} the magnitude of the TDDFT SS shift strongly depends on the chosen functional, and it goes from −0.42 eV to −0.27 eV upon replacing the PBE0 functional by CAMB3LYP, which is more suited for such an excitedstate. The BSE (Ω − Ω_{0}) shift originates mainly from its SS component as well and lies in between the PBE0 and CAMB3LYP values, though much closer to the latter, as expected.
The same conclusions are reached when considering the intermolecular CT transition in the benzene–TCNE complex. First, we observe that the BSE formalism very nicely reproduces the gas phase excitation energy that is available experimentally. The TDCAMB3LYP calculation also provide a reasonable value, while the TDPBE0 approach yields a much too small Ω_{g}, a logical consequence of its lack of longrange corrections. As expected for TDDFT, the LR approach again provides a negligible (Ω − Ω_{0}) dynamic shift, while the SS formalism yields a large redshift. Again, the BSE + PCM formalism captures the correct physics, with a negligible LR contribution and a large SS contribution. Such calculation clearly demonstrates that the proposed BSE + PCM approach can also, following eqn (17b), capture SS dynamic shifts.
4.3 Hybrid ES
We now turn to difficult cases where neither the LR nor SS contributions can be neglected, and our results are collected in Table 3. As such, TDDFT calculations adopting one or the other response formalisms cannot capture the entire solvatochromic shifts, though, as discussed in ref. 26 both contributions should in principle be accounted for. This is first evidenced by considering the planar (standard) PNA molecule. The much smaller disagreement between TDPBE0 and TDCAMB3LYP, as compared to the PNA_{perp} case, is indeed a first indication that the CT character is here significantly reduced. As a result (see Table 3), both SS and LR contributions to the (Ω − Ω_{0}) shift are large: they are of equal magnitude within TDPBE0 while the SS contribution is larger with TDCAMB3LYP. Clearly, the BSE formalism simultaneously accounts for both contributions, the SS contribution being slightly larger than the LR one, consistent with the TDDFT results. The total (Ω − Ω_{0}) dynamic shift amounts to −0.18 eV and −0.27 eV when adding the LR and cLR contribution determined at TDPBE0 and TDCAMB3LYP, respectively, comparable to the −0.19 eV solvatochromic shift obtained at the BSE level.
Table 3 Data as in Table 1 but for the transitions of mixed character

Ω
_{g}

Ω(Ω − Ω_{g}) 
Ω
_{0} − Ω_{g} 
Ω − Ω_{0} 
LR 
SS 
PNA in water (partial CT excitation)

TDPBE0 (LR) 
4.202 
3.802(−0.400) 
−0.314 
−0.086 

TDPBE0 (cLR) 
4.202 
3.796(−0.406) 
−0.314 

−0.092 
TDCAMB3LYP (LR) 
4.513 
4.105(−0.408) 
−0.321 
−0.087 

TDCAMB3LYP (cLR) 
4.513 
4.005(−0.508) 
−0.321 

−0.187 
BSE 
4.527 
3.864(−0.663) 
−0.470 
−0.090 
−0.102 

4Nitropyridine Noxide in benzene (mixed excitation)

TDPBE0 (LR) 
3.989 
3.815(−0.174) 
−0.048 
−0.126 

TDPBE0 (cLR) 
3.989 
3.797(−0.192) 
−0.048 

−0.144 
TDCAMB3LYP (LR) 
4.196 
4.023(−0.173) 
−0.033 
−0.140 

TDCAMB3LYP (cLR) 
4.196 
4.109(−0.087) 
−0.033 

−0.054 
BSE 
3.966 
3.658(−0.267) 
−0.001 
−0.193 
−0.073 

4Nitropyridine Noxide in water (mixed excitation)

TDPBE0 (LR) 
3.989 
3.797(−0.192) 
−0.097 
−0.095 

TDPBE0 (cLR) 
3.989 
3.827(−0.162) 
−0.097 

−0.065 
TDCAMB3LYP (LR) 
4.196 
4.028(−0.168) 
−0.065 
−0.103 

TDCAMB3LYP (cLR) 
4.196 
4.067(−0.129) 
−0.065 

−0.064 
BSE 
3.966 
3.687(−0.279) 
−0.070 
−0.141 
−0.068 
4Nitropyridine Noxide is an organic probe used to assess the nature of solvents following a Kamel–Taft type of analysis.^{103} A previous throughout theoretical analysis of the solvatochromism of this probe is available,^{104} and demonstrates that, none of the available LR or SS PCM model is able to describe the solvent effects in a TDDFT context, as the observed spectral changes come from a fine interplay of several effects. For this compound, the obtained conclusions are similar to the PNA case: both the LR and SS contributions are nonnegligible even though the SS dynamic shift is somehow larger. Adding the LR and SS TDDFT contributions, the overall (Ω − Ω_{0}) shift amounts to −0.16 eV and −0.17 eV with PBE0 and CAMB3LYP in water, respectively, a shift that can be compared to the −0.21 eV value obtained within BSE. The experimental gasphase excitation energy obtained through extrapolation in ref. 104 is 3.80 eV, and one notices that the BSE value is reasonably close to that estimate. In benzene, the measurement gives 3.52 eV,^{103} corresponding to a solvatochromic shift of −0.28 eV, a value that BSE can reproduce (−0.27 eV), whereas TDDFT in unable to do so. In water, the hydrogen bonds with the negatively charged oxygen atom of the probe play a crucial role, and a blueshift is observed (the excitation maximum takes place at 3.94 eV),^{103} an effect that all continuum approaches logically fail to capture.
5 Further discussion
We summarise in Fig. 3 the principal findings of the present study, namely that the BSE/GW formalism accounts for both LR and SS shifts, a fact demonstrated for a large variety of transitions, including excitations where one of the two terms is dominant and excitations for which both are needed. While the partitioning between BSE LR and SS terms is provided in eqn (17b)–(17c), it is interesting to analyse now the contributions to the BSE dynamic shift originating from the solventinduced renormalisation of the ionisation potential and electronic affinity, namely the variation of the GW HOMO–LUMO gap in the solvent, and from the reduction of the electron–hole interaction contained in the 〈abWij〉 terms.

 Fig. 3 Schematic representation of the (Ω − Ω_{0}) dynamic solvent induced shifts for the π–π* excitation in acrolein, the CT excitation in the benzene–TCNE complex, and the (planar) PNA mixed excitation. By dynamic shift, we mean the effect of switching the fast v^{reac}(ε_{∞}) PCM reaction field on top of the (slow) groundstate v^{reac}(ε_{0}) PCM response. The represented data correspond to the (Ω − Ω_{0}) shifts given in Tables 1–3.  
For the π–π* transition in acrolein, with a shift captured by the LR scheme only at the TDDFT level, the decomposition of the (Ω − Ω_{0}) energy difference confirms that the variation of the BSE SSlike 〈Δε − W〉 component (eqn (17b)) does not contribute to the shift. Such a result may seem surprising since, as expected, the 〈Δε^{GW}〉 HOMO–LUMO gap becomes smaller by 2.34 eV upon “switching” v^{reac}(ε_{∞}), consistent with an (absolute) polarisation energy of ca. 1.1–1.2 eV that affects the IP and AE of hydrated compounds compared to the gas phase.^{54} However, very remarkably, the 〈W〉 electron–hole binding energy is also decreased by 2.35 eV. Namely, the screening by the fast reaction field reduces similarly the occupiedtovirtual (ε_{a}^{GW} − ε_{i}^{GW}) GW gap and the electron–hole binding energy (see Fig. 1). As such, only the 〈v^{reac}(ε_{∞})〉 LR terms (eqn (17c)) explain the solvatochromic effect, consistent with the TDDFT results.
We now turn to the opposite situation of a CT excitation for which the solvatochromic shift can only be described by SS implementations, such as cLR, within TDDFT. In the case of the benzene–TCNE complex, the LR 〈iav^{reac}(ε_{∞})bj〉 contributing matrix elements become very small, as a logical consequence of the small overlap between the donor and acceptor wavefunctions. As such the only contribution to the (Ω − Ω_{0}) shift originates from the variation of the SSlike 〈Δε − W〉 term. The impact of v^{reac}(ε_{∞}) leads to the reduction of the 〈Δε〉 average gap by 1.60 eV, while the electron–hole binding energy 〈W〉 reduces by a smaller 1.31 eV amount, accounting for most of the (Ω − Ω_{0}) = −0.29 eV redshift reported in Table 2.
An important consequence of the present analysis, showing that, in BSE/GW theory, SS shifts result from the competition between the reduction of the GW occupiedtovirtual (ε_{a}^{GW} − ε_{i}^{GW}) energy gaps and the electron–hole 〈abWij〉 binding energies, is that both electron–electron and electron–hole interactions must be treated on the same footing, namely here through the screened Coulomb potential W. In particular, BSE + PCM calculations starting from KS eigenstates generated with exchangecorrelation functionals optimallytuned, so as to generate the correct gas phase HOMO–LUMO gap, might not deliver the correct SS contribution to optical excitation energy shifts.
Finally, while Table 2 indicates a large variability of the SS (cLR) TDDFT (Ω − Ω_{0}) dynamic shift as a function of the chosen XC functional, it is interesting to emphasize the stability of the BSE/evGW data with respect to the input KS eigenstates used to build the evGW electronic energy levels and the screened Coulomb potential W. This is illustrated in Fig. 4a where we plot the PNA_{perp} excitation energy in water with the cLRTDPBE(α) (open triangles) and BSE/evGW@PBE(α) (open circles) methods. Here α indicates the percentage of exact exchange (EEX) used in the hybrid functional and going from 0% (PBE) to 100% (HFPBE) in Fig. 4. While the TDDFT excitation energy is shown to increase steeply as a function of α, the BSE data are much more stable, with a variation of ca. 0.12 eV from (α = 0%) to (α = 90%). The stability of gas phase BSE/evGW excitation energies with respect to the starting KS starting point has been documented in previous benchmark studies,^{47,49} and is therefore shown here to pertain for condensedphase calculations. In fact, plotting now in Fig. 4b the dynamic (Ω − Ω_{0}) solvatochromic shifts, we observe again a large variability of the TDPBE(α) cLR shifts (the LR contribution is vanishingly small), while again the BSE shifts are extremely stable, evolving from −0.29 eV to −0.25 eV with α. To rationalize this result, we recall that BSE calculations are performed on top of partially selfconsistent evGW calculations where the corrected electronic energy levels are selfconsistently reinjected during the construction of G and W. This selfconsistent treatment of the quasiparticle energies {ε_{n}^{GW}} and screened Coulomb potential W leads to a large stability of the BSE Hamiltonian and reaction field, with a small residual dependency on the starting KS eigenstates due to the fact that the {ϕ_{n}} eigenstates are unchanged. Such a stability of the BSE/evGW excitation energies and solvent induced shifts, that allows us to alleviate the standard problem of the proper choice of the XC functional central in TDDFT, is one of the interesting features of the present scheme. As a fair tribute to TDDFT, we observe however that the cLR shift matches closely that of BSE for a large exact exchange ratio. It is indeed known that large amounts of exact exchange are required in TDDFT to obtain an accurate description of the groundtoexcited density variation for CT excitations.^{17}

 Fig. 4 (a) Evolution of the PNA_{perp} excitation energy Ω in solution as a function of the percentage of exact exchange α in the PBE(α) functional. Triangles correspond to the cLRTDDFT results while circles indicate BSE/evGW@PBE(α) data. (b) Dynamic shift (Ω − Ω_{0}) obtained with cLRTDDFT (triangles) and BSE/evGW@PBE(α) (circles). The inconsistent value appearing with TDDFT for α = 40% is due to the fact that two excited states of mixed character are nearly degenerate for that EEX ratio. Energies are in eV.  
6 Conclusions
We have demonstrated that the Bethe–Salpeter equation (BSE) formalism combined with the PCM approach straightforwardly includes, on the same footing, both linearresponse (LR) and statespecific (SS) contributions to the solvatochromic shifts. This result was established on analytic grounds as well as by comparing the numerical results obtained for several molecules or complexes using BSE + PCM or TDDFT + PCM approaches, the latter being performed with both the LR and cLR formalisms. To the best of our knowledge, such a consistent theoretical approach simultaneously providing both LR and SS solvent shifts remains unavailable in the TDDFT world, and has been developed previously, only for theories presenting a less favourable scaling with system size, e.g., ADC(2).^{26}
Beyond the standard LR 〈abv^{reac}(ε_{∞})ij〉 matrix element contributions, where v^{reac}(ε_{∞}) is the PCM reaction field to the electronic excitation, the proper inclusion of SS shifts hinges on the incorporation of v^{reac}(ε_{∞}) in the screened Coulomb potential W. This dressed screened Coulomb potential renormalises both electron–electron and electron–hole interactions, namely both the GW quasiparticle energies for occupied/virtual electronic levels and the BSE electron–hole (excitonic) binding energy. Groundstate polarisation effects are accounted for, in the present nonequilibrium scheme, by starting our BSE/GW calculations with KS eigenstates generated at the DFT + PCM(ε_{0}) level.
Following previous studies related to merging the GW formalism with discrete polarisable models,^{56–58} the same analysis can be straightforwardly applied to BSE calculations combined with other polarisation models.^{55,57,105} Specifically, the use of an explicit (molecular) description of the environment, combined with, e.g., empirical force fields, should allow us to account for groundstate electrostatic field effects induced by the solvent molecule static multipoles, an important contribution in the case of polar solvents that cannot be captured by the PCM model. As shown in the case of the π–π* transition in aqueous acrolein,^{89,91} such contributions can explain the difference between the experimental shift (−0.53 eV) and the ∼−0.25 eV redshifts obtained with the present BSE or TDDFTLR formalisms combined with the PCM.
A From Δε^{GW} to the Born solvation energy model
In this Appendix, we justify eqn (18a). As shown in ref. 54, the variation of the ε_{n}^{GW} quasiparticle energies can be simply expressed within the socalled static Coulombhole (COH) plus screenedexchange (SEX) approximation to the selfenergy Σ^{GW}, namely 
ΔΣ^{GW} ≃ ΔΣ^{SEX} + ΔΣ^{COH},  (25) 
with the following matrix elements: 
 (26) 

 (27) 
where W is taken in the lowfrequency optical limit and = W − v. We therefore get for the variation Δε_{n}^{GW}: 
 (28) 
where {ϕ_{n}} represents the frozen groundstate eigenstates generated at the DFT + PCM(ε_{0}) level. Here, one can expect the principal contribution to eqn (28) to come from the reaction field induced by the monopolar terms associated with the unit charge density ϕ^{*}_{n}(r)ϕ_{n}(r) = ϕ_{n}(r)^{2}, i.e., 
 (29) 
with the + or − sign depending on whether n describes an empty (a) or an occupied (i) state. All other contributions are expected to be smaller, implying most dipole–dipole interactions between neutral codensities since, by orthogonalization, and . A similar argument was provided and validated by Neaton and coworkers in their study of the screening of molecular energy levels by metallic electrodes,^{106} in which they assume that ΔW is slowly varying over the solute cavity, that is, the monopolar term is dominant.
B Relation between ΔW and v^{reac}
Working out eqn (14) and (16) to lower order in v^{reac} = vχ^{PCM}v leads to:
with χ_{0}, ε = 1 − vχ_{0}, and W being the solute freeelectron susceptibility, dielectric matrix and screened Coulomb potential in the absence of the PCM. Neglecting the internal response of the QM solute in front of the PCM “reservoir”, namely taking χ_{0} → 0, one obtains that indeed ΔW is equal to the reaction field.
Conflicts of interest
There are no conflicts to declare.
Acknowledgements
This project has received funding from the European Union Horizon 2020 research and innovation programme under grant agreement no. 646176 (EXTMOS). Calculations were performed using resources from the GENCICINES/IDRIS French national supercomputing centers, the CCIPL (Centre de Calcul Intensif des Pays de la Loire), and Grenoble CIMENT (Calcul Intensif, Modélisation, Expérimentation Numérique et Technologique) computational mesocenters. C. A. G and D. J. are indebted to the LumoMat consortium for support through the FCPolResp project.
References
 M. E. Casida, J. Mol. Struct.: THEOCHEM, 2009, 914, 3–18 CrossRef CAS .

C. Ullrich, TimeDependent DensityFunctional Theory: Concepts and Applications, Oxford University Press, New York, 2012 Search PubMed .
 C. van Caillie and R. D. Amos, Chem. Phys. Lett., 1999, 308, 249–255 CrossRef CAS .
 F. Furche and R. Ahlrichs, J. Chem. Phys., 2002, 117, 7433–7447 CrossRef CAS .
 G. Scalmani, M. J. Frisch, B. Mennucci, J. Tomasi, R. Cammi and V. Barone, J. Chem. Phys., 2006, 124, 094107 CrossRef PubMed .
 J. Liu and W. Z. Liang, J. Chem. Phys., 2011, 135, 184111 CrossRef PubMed .
 C. A. Guido, G. Scalmani, B. Mennucci and D. Jacquemin, J. Chem. Phys., 2017, 146, 204106 CrossRef PubMed .
 S. Miertǔs, E. Scrocco and J. Tomasi, Chem. Phys., 1981, 55, 117–129 CrossRef .
 E. Cancès, B. Mennucci and J. Tomasi, J. Chem. Phys., 1997, 107, 3032–3041 CrossRef .
 R. Cammi and B. Mennucci, J. Chem. Phys., 1999, 110, 9877–9886 CrossRef CAS .
 M. Cossi and V. Barone, J. Chem. Phys., 2001, 115, 4708–4717 CrossRef CAS .
 R. Cammi, S. Corni, B. Mennucci and J. Tomasi, J. Chem. Phys., 2005, 122, 104513 CrossRef CAS PubMed .
 M. Caricato, B. Mennucci, J. Tomasi, F. Ingrosso, R. Cammi, S. Corni and G. Scalmani, J. Chem. Phys., 2006, 124, 124520 CrossRef PubMed .
 R. Improta, V. Barone, G. Scalmani and M. J. Frisch, J. Chem. Phys., 2006, 125, 054103 CrossRef PubMed .
 A. V. Marenich, C. J. Cramer, D. G. Truhlar, C. A. Guido, B. Mennucci, G. Scalmani and M. J. Frisch, Chem. Sci., 2011, 2, 2143–2161 RSC .
 A. Pedone, J. Chem. Theory Comput., 2013, 9, 4087–4096 CrossRef CAS PubMed .
 C. A. Guido, D. Jacquemin, C. Adamo and B. Mennucci, J. Chem. Theory Comput., 2015, 11, 5782–5790 CrossRef CAS PubMed .
 C. A. Guido, B. Mennucci, G. Scalmani and D. Jacquemin, J. Chem. Theory Comput., 2018, 14, 1544–1553 CrossRef CAS PubMed .
 S. Corni, R. Cammi, B. Mennucci and J. Tomasi, J. Chem. Phys., 2005, 123, 134512 CrossRef CAS PubMed .
 R. Guareschi, O. Valsson, C. Curutchet, B. Mennucci and C. Filippi, J. Phys. Chem. Lett., 2016, 7, 4547–4553 CrossRef CAS PubMed .
 R. Cammi, R. Fukuda, M. Ehara and H. Nakatsuji, J. Chem. Phys., 2010, 133, 024104 CrossRef PubMed .
 M. Caricato, B. Mennucci, G. Scalmani, G. W. Trucks and M. J. Frisch, J. Chem. Phys., 2010, 132, 084102 CrossRef PubMed .
 R. Fukuda, M. Ehara, H. Nakatsuji and R. Cammi, J. Chem. Phys., 2011, 134, 104109 CrossRef PubMed .
 M. Caricato, J. Chem. Theory Comput., 2012, 8, 4494–4502 CrossRef CAS PubMed .
 M. Caricato, J. Chem. Phys., 2013, 139, 044116 CrossRef PubMed .
 B. Lunkenheimer and A. Köhn, J. Chem. Theory Comput., 2013, 9, 977–994 CrossRef CAS PubMed .
 M. Caricato, Comput. Theor. Chem., 2014, 1040, 99–105 CrossRef .
 R. Fukuda and M. Ehara, J. Chem. Phys., 2014, 141, 154104 CrossRef PubMed .
 J.M. Mewes, Z.Q. You, M. Wormit, T. Kriesche, J. M. Herbert and A. Dreuw, J. Phys. Chem. A, 2015, 119, 5446–5464 CrossRef CAS PubMed .
 J.M. Mewes, J. M. Herbert and A. Dreuw, Phys. Chem. Chem. Phys., 2017, 19, 1644–1654 RSC .
 E. E. Salpeter and H. A. Bethe, Phys. Rev., 1951, 84, 1232–1242 CrossRef .
 L. J. Sham and T. M. Rice, Phys. Rev., 1966, 144, 708–714 CrossRef CAS .
 W. Hanke and L. J. Sham, Phys. Rev. Lett., 1979, 43, 387–390 CrossRef CAS .
 G. Strinati, Phys. Rev. Lett., 1982, 49, 1519–1522 CrossRef CAS .
 X. Blase, I. Duchemin and D. Jacquemin, Chem. Soc. Rev., 2018, 47, 1022–1043 RSC .
 D. Rocca, D. Lu and G. Galli, J. Chem. Phys., 2010, 133, 164109 CrossRef PubMed .
 J. M. GarciaLastra and K. S. Thygesen, Phys. Rev. Lett., 2011, 106, 187402 CrossRef CAS PubMed .
 X. Blase and C. Attaccalite, Appl. Phys. Lett., 2011, 99, 171909 CrossRef .
 I. Duchemin, T. Deutsch and X. Blase, Phys. Rev. Lett., 2012, 109, 167801 CrossRef CAS PubMed .
 B. Baumeier, D. Andrienko and M. Rohlfing, J. Chem. Theory Comput., 2012, 8, 2790–2795 CrossRef CAS PubMed .
 S. Sharifzadeh, P. Darancet, L. Kronik and J. B. Neaton, J. Phys. Chem. Lett., 2013, 4, 2197–2201 CrossRef CAS .
 P. Cudazzo, M. Gatti, A. Rubio and F. Sottile, Phys. Rev. B, 2013, 88, 195152 CrossRef .
 V. Ziaei and T. Bredow, J. Chem. Phys., 2016, 145, 174305 CrossRef PubMed .
 D. Escudero, I. Duchemin, X. Blase and D. Jacquemin, J. Phys. Chem. Lett., 2017, 8, 936–940 CrossRef CAS PubMed .
 P. Boulanger, D. Jacquemin, I. Duchemin and X. Blase, J. Chem. Theory Comput., 2014, 10, 1212–1218 CrossRef CAS PubMed .
 C. Azarias, I. Duchemin, X. Blase and D. Jacquemin, J. Chem. Phys., 2017, 146, 034301 CrossRef PubMed .
 D. Jacquemin, I. Duchemin and X. Blase, J. Chem. Theory Comput., 2015, 11, 3290–3304 CrossRef CAS PubMed .
 F. Bruneval, S. M. Hamed and J. B. Neaton, J. Chem. Phys., 2015, 142, 244101 CrossRef PubMed .
 D. Jacquemin, I. Duchemin and X. Blase, J. Chem. Theory Comput., 2015, 11, 5340–5359 CrossRef CAS PubMed .
 D. Jacquemin, I. Duchemin, A. Blondel and X. Blase, J. Chem. Theory Comput., 2016, 12, 3969–3981 CrossRef CAS PubMed .
 D. Jacquemin, I. Duchemin, A. Blondel and X. Blase, J. Chem. Theory Comput., 2017, 13, 767–783 CrossRef CAS PubMed .
 T. Rangel, S. M. Hamed, F. Bruneval and J. B. Neaton, J. Chem. Phys., 2017, 146, 194108 CrossRef PubMed .
 D. Jacquemin, I. Duchemin and X. Blase, J. Phys. Chem. Lett., 2017, 8, 1524–1529 CrossRef CAS PubMed .
 I. Duchemin, D. Jacquemin and X. Blase, J. Chem. Phys., 2016, 144, 164106 CrossRef PubMed .
 B. Baumeier, M. Rohlfing and D. Andrienko, J. Chem. Theory Comput., 2014, 10, 3104–3110 CrossRef CAS PubMed .
 J. Li, G. D'Avino, I. Duchemin, D. Beljonne and X. Blase, J. Phys. Chem. Lett., 2016, 7, 2814–2820 CrossRef CAS PubMed .
 J. Li, G. D'Avino, A. Pershin, D. Jacquemin, I. Duchemin, D. Beljonne and X. Blase, Physical Review Materials, 2017, 1, 025602 CrossRef .
 J. Li, G. D'Avino, I. Duchemin, D. Beljonne and X. Blase, Phys. Rev. B, 2018, 97, 035108 CrossRef .
 M. Rohlfing and S. G. Louie, Phys. Rev. Lett., 1998, 80, 3320–3323 CrossRef CAS .
 L. X. Benedict, E. L. Shirley and R. B. Bohn, Phys. Rev. Lett., 1998, 80, 4514–4517 CrossRef CAS .
 S. Albrecht, L. Reining, R. Del Sole and G. Onida, Phys. Rev. Lett., 1998, 80, 4510–4513 CrossRef CAS .
 G. Onida, L. Reining and A. Rubio, Rev. Mod. Phys., 2002, 74, 601–659 CrossRef CAS .
 G. Bussi, Phys. Scr., 2004, 2004, 141 CrossRef .
 L. Hedin, Phys. Rev., 1965, 139, A796–A823 CrossRef .
 M. S. Hybertsen and S. G. Louie, Phys. Rev. B, 1986, 34, 5390–5413 CrossRef CAS .
 R. W. Godby, M. Schlüter and L. J. Sham, Phys. Rev. B, 1988, 37, 10159–10175 CrossRef .
 B. Farid, R. Daling, D. Lenstra and W. van Haeringen, Phys. Rev. B, 1988, 38, 7530–7534 CrossRef .
 C. Faber, C. Attaccalite, V. Olevano, E. Runge and X. Blase, Phys. Rev. B, 2011, 83, 115123 CrossRef .
 K. Krause, M. E. Harding and W. Klopper, Mol. Phys., 2015, 113, 1952–1960 CrossRef CAS .
 F. Kaplan, M. E. Harding, C. Seiler, F. Weigend, F. Evers and M. J. van Setten, J. Chem. Theory Comput., 2016, 12, 2528–2541 CrossRef CAS PubMed .
 J. W. Knight, X. Wang, L. Gallandi, O. Dolgounitcheva, X. Ren, J. V. Ortiz, P. Rinke, T. Körzdörfer and N. Marom, J. Chem. Theory Comput., 2016, 12, 615–626 CrossRef CAS PubMed .
 T. Rangel, S. M. Hamed, F. Bruneval and J. B. Neaton, J. Chem. Theory Comput., 2016, 12, 2834–2842 CrossRef CAS PubMed .
 T. Stein, L. Kronik and R. Baer, J. Am. Chem. Soc., 2009, 131, 2818–2820 CrossRef CAS PubMed .
 I. Duchemin, J. Li and X. Blase, J. Chem. Theory Comput., 2017, 13, 1199–1208 CrossRef CAS PubMed .
 M. Valiev, E. Bylaska, N. Govind, K. Kowalski, T. Straatsma, H. V. Dam, D. Wang, J. Nieplocha, E. Apra, T. Windus and W. de Jong, Comput. Phys. Commun., 2010, 181, 1477–1489 CrossRef CAS .
 T. H. Dunning Jr, J. Chem. Phys., 1989, 90, 1007–1023 CrossRef .
 J. P. Perdew, M. Ernzerhof and K. Burke, J. Chem. Phys., 1996, 105, 9982–9985 CrossRef CAS .
 C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170 CrossRef CAS .
 F. Weigend, J. Chem. Phys., 2002, 116, 3175–3183 CrossRef CAS .
 X. Blase, C. Attaccalite and V. Olevano, Phys. Rev. B, 2011, 83, 115103 CrossRef .
 C. Faber, J. L. Janssen, M. Côté, E. Runge and X. Blase, Phys. Rev. B, 2011, 84, 155104 CrossRef .
 A. Klamt, C. Moya and J. Palomar, J. Chem. Theory Comput., 2015, 11, 4220–4225 CrossRef CAS PubMed .

M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. WilliamsYoung, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr, J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian16 Revision A.03, Gaussian Inc., Wallingford CT, 2016 Search PubMed .
 T. Yanai, D. P. Tew and N. C. Handy, Chem. Phys. Lett., 2004, 393, 51–57 CrossRef CAS .
 C. A. Guido, P. Cortona and C. Adamo, J. Chem. Phys., 2014, 140, 104101 CrossRef PubMed .
 C. A. Guido, P. Cortona, B. Mennucci and C. Adamo, J. Chem. Theory Comput., 2013, 9, 3118–3126 CrossRef CAS PubMed .
 A. Moskvin, O. Yablonskii and L. Bondar, Theor. Exp. Chem., 1966, 2, 469 CrossRef .
 K. Aidas, A. Møgelhøj, E. J. K. Nilsson, M. S. Johnson, K. V. Mikkelsen, O. Christiansen, P. Söderhjelm and J. Kongsted, J. Chem. Phys., 2008, 128, 194503 CrossRef PubMed .
 J. M. Olsen, K. Aidas and J. Kongsted, J. Chem. Theory Comput., 2010, 6, 3721–3734 CrossRef CAS .
 A. H. Steindal, K. Ruud, L. Frediani, K. Aidas and J. Kongsted, J. Phys. Chem. B, 2011, 115, 3027–3037 CrossRef CAS PubMed .
 K. Sneskov, E. Matito, J. Kongsted and O. Christiansen, J. Chem. Theory Comput., 2010, 6, 839–850 CrossRef CAS PubMed .
 S. E. Sheppard and P. T. Newsome, J. Am. Chem. Soc., 1942, 64, 2937–3946 CrossRef CAS .
 P. W. Saddler, J. Org. Chem., 1956, 21, 316–318 CrossRef .
 M. Klessinger and W. Lüttke, Chem. Ber., 1966, 99, 2136–2145 CrossRef CAS .
 E. Wille and W. Lüttke, Angew. Chem., Int. Ed., 1971, 10, 803–804 CrossRef CAS .
 A. R. Monahan and J. E. Kuder, J. Org. Chem., 1972, 37, 4182–4184 CrossRef CAS .
 G. Haucke and G. Graness, Angew. Chem., Int. Ed., 1995, 34, 67–68 CrossRef CAS .
 R. Gerken, L. Fitjer, P. Müller, I. Usón, K. Kowski and P. Rademacher, Tetrahedron, 1999, 55, 14429–14434 CrossRef .
 D. Jacquemin, J. Preat, V. Wathelet and E. A. Perpète, J. Chem. Phys., 2006, 124, 074104 CrossRef PubMed .
 R. Cammi, L. Frediani, B. Mennucci and K. Ruud, J. Chem. Phys., 2003, 119, 5818–5827 CrossRef CAS .
 S. Millefiori, G. Favini, A. Millefiori and D. Grasso, Spectrochim. Acta, Part A, 1977, 33, 21–27 CrossRef .
 I. Hanazaki, J. Phys. Chem., 1972, 76, 1982–1989 CrossRef CAS .
 A. F. Lagalante, R. J. Jacobson and T. J. Bruno, J. Org. Chem., 1996, 61, 6404–6406 CrossRef CAS PubMed .
 S. Budzak, A. D. Laurent, C. Laurence, M. Miroslav and D. Jacquemin, J. Chem. Theory Comput., 2016, 12, 1919–1929 CrossRef CAS PubMed .
 D. Varsano, S. Caprasecca and E. Coccia, J. Phys.: Condens. Matter, 2017, 29, 013002 CrossRef PubMed .
 J. B. Neaton, M. S. Hybertsen and S. G. Louie, Phys. Rev. Lett., 2006, 97, 216405 CrossRef CAS PubMed .
Footnote 
† Electronic supplementary information (ESI) available: Cartesian coordinates of the compounds. See DOI: 10.1039/c8sc00529j 

This journal is © The Royal Society of Chemistry 2018 