A quantum-mechanical perspective on linear response theory within polarizable embedding

TL;DR

This paper provides a rigorous quantum-mechanical foundation for linear response theory in polarizable embedding, analyzing subsystem decompositions, response function structures, and implications for molecular property calculations.

Contribution

It introduces a quantum-mechanical derivation of linear response theory in polarizable embedding, clarifies subsystem response functions, and justifies classical approximations from quantum principles.

Findings

Analysis of symmetric and nonsymmetric response function decompositions

Derivation of conservation laws for absorption cross sections

Identification of steps to connect quantum and classical descriptions

Abstract

The derivation of linear response theory within polarizable embedding is carried out from a rigorous quantum-mechanical treatment of a composite system. Two different subsystem decompositions (symmetric and nonsymmetric) of the linear response function are presented, and the pole structures as well as residues of the individual terms are analyzed and discussed. This theoretical analysis clarifies which form of the response function to use in polarizable embedding, and we highlight complications in separating out subsystem contributions to molecular properties. For example, based on the nonsymmetric decomposition of the complex linear response function, we derive conservation laws for integrated absorption cross sections, providing a solid basis for proper calculations of the intersubsystem intensity borrowing inherent to coupled subsystems and how that can lead to negative subsystem intensities. We finally identify steps and approximations required to achieve the transition from a quantum-mechanical description of the composite system to polarizable embedding with a classical treatment of the environment, thus providing a thorough justification for the descriptions used in polarizable embedding models.

Tables

Table 1. Table 1: The polarized subsystem parameters used in the SLM for the p NA–water complex in Fig.  2 . All results are reported in a.u. and have been obtained at the CAM-B3LYP/aug-cc-pVDZ level of theory in the presence of the ground-state-frozen PE embedding potential of the other subsystem.

State	$\Delta E$	$\left|{\mu }_x\right|$	$\left|{\mu }_y\right|$	$\left|{\mu }_z\right|$
$|n{\pi }^{\ast }⟩$	0.14014	0.00102	0.00017	0.00023
$|\pi {\pi }^{\ast }⟩$	0.15428	2.02274	0.00030	0.00539
$|1B_1⟩$	0.26814	0.07302	0.02634	0.54089
$|1A_1⟩$	0.35027	0.42522	0.49612	0.09183

Table 2. Table 2: Absolute (unscaled) electronic couplings used in the SLM for the p NA–water complex in Fig.  2 . Digits in parenthesis refer to negative exponents, i.e., a ​ ( b ) = a × 10 − b 𝑎 𝑏 𝑎 superscript 10 𝑏 a(b)=a\times 10^{-b} . All results are reported in a.u. and have been obtained at the CAM-B3LYP/aug-cc-pVDZ level of theory in the presence of the ground-state-frozen PE embedding potential of the other subsystem.

pNA$\backslash$water	$|1B_1⟩$	$|1A_1⟩$
$|n{\pi }^{\ast }⟩$	2.53(6)	1.03(6)
$|\pi {\pi }^{\ast }⟩$	3.82(6)	1.76(3)

Table 3. Table 3: Excitation energies and one-photon transition strengths for the two lowest-lying (i.e., subsystem A -dominated) electronic transitions in the composite system as computed from the full treatment as well as the nonsymmetric decomposition (NSD), including and excluding various contributions. All results are reported in a.u. and digits in parenthesis refer to negative exponents, i.e., a ​ ( b ) = a × 10 − b 𝑎 𝑏 𝑎 superscript 10 𝑏 a(b)=a\times 10^{-b} . vac: vacuum, FP: frozen ground-state environment polarization, renorm: renormalization, EEF: effective external field, ZF: zero-frequency approximation for the environment.

	$|1_{\text{AB}}⟩$				$|2_{\text{AB}}⟩$
Model	$\Delta E$	${\left|{\mu }_x\right|}^2$	${\left|{\mu }_y\right|}^2$	${\left|{\mu }_z\right|}^2$	$\Delta E$	${\left|{\mu }_x\right|}^2$	${\left|{\mu }_y\right|}^2$	${\left|{\mu }_z\right|}^2$
vac	0.13950	6.0205(7)	3.3384(8)	6.6642(9)	0.15625	4.0014	6.5832(7)	4.3947(5)
FP	0.14014	1.4022(6)	3.0284(8)	5.1435(8)	0.15428	4.0915	8.8700(8)	2.9025(5)
Full system	0.14014	1.9202(6)	4.8179(8)	2.2571(9)	0.15110	4.3996	5.4416(3)	3.5539(4)
NSD (Eq. (45))		1.9202(6)	4.8179(8)	2.2571(9)		4.3996	5.4416(3)	3.5539(4)
NSD($-$renorm)		1.9202(6)	4.8179(8)	2.2571(9)		4.4423	5.4944(3)	3.5884(4)
NSD($-$EEF)		1.8561(6)	3.0263(8)	5.1047(8)		4.1373	8.9645(8)	2.9350(5)
NSD($-$(renorm/EEF))		1.8561(6)	3.0263(8)	5.1047(8)		4.1774	9.0515(8)	2.9635(5)
NSD(+ZF)	0.14014	1.7309(6)	4.4586(8)	9.1521(9)	0.15170	4.3745	3.6180(3)	2.7005(4)

Equations

$$\hat{\mathcal{H}}=\sum _{\text{A}=1}^N{\hat{\mathcal{H}}}_{\text{A}}+\sum _{\text{A}>\text{B}}^N{\hat{\mathcal{V}}}_{\text{AB}}\text{,}$$

$${\hat{\mathcal{V}}}_{\text{AB}}=$$

$$+\sum _{pq\in \text{A}}\sum _{rs\in \text{B}}{v}_{pq,rs}{\hat{E}}_{pq}{\hat{E}}_{rs}+\sum _{n\in \text{A}}^{M_{\text{A}}}\sum _{m\in \text{B}}^{M_{\text{B}}}\frac{Z_n{Z}_m}{|{\mathbf{R}}_n-{\mathbf{R}}_m|}\text{,}$$

$$v_{pq}(\mathbf{R})$$

$$v_{pq,rs}$$

$$\displaystyle=-\int\phi_{p}^{*}(\mathbf{r}_{1})\phi_{q}(\mathbf{r}_{1})v_{rs}(\mathbf{r}_{1})\ \text{d}\mathbf{r}_{1}\ \text{.}\$$

$${\hat{\mathcal{V}}}_{\text{AB}}=\int {\hat{\rho }}_{\text{A}}(\mathbf{r}){\hat{\mathcal{V}}}_{\text{B}}(\mathbf{r})\text{d}\mathbf{r}=\int {\hat{\rho }}_{\text{B}}(\mathbf{r}){\hat{\mathcal{V}}}_{\text{A}}(\mathbf{r})\text{d}\mathbf{r}={\hat{\rho }}_{\mathbf{r}}^{\text{A}}{\hat{\mathcal{V}}}_{\mathbf{r}}^{\text{B}}\text{,}$$

$${\hat{\rho }}_{\text{A}}(\mathbf{r})=$$

$${\hat{\mathcal{V}}}_{\text{B}}(\mathbf{r})=$$

$$|{\Psi }_{\text{A}}{\Psi }_{\text{B}}\dots {\Psi }_{\text{N}}⟩$$

$$\forall \text{A}:{\hat{\mathcal{F}}}_{\text{A}}|0_{\text{A}}⟩=E_{\text{A}}|0_{\text{A}}⟩,$$

$${\hat{\mathcal{F}}}_{\text{A}}$$

$$\displaystyle=\hat{\mathcal{H}}_{{\textit{\tiny A}}}+\sum_{pq\in{{\textit{\tiny A}}}}\sum_{{{\textit{\tiny B}}\neq{\textit{\tiny A}}}}\biggl{[}\sum_{m\in{{\textit{\tiny B}}}}^{M_{{\textit{\tiny B}}}}Z_{m}v_{pq}(\mathbf{R}_{m})+\sum_{rs\in{{\textit{\tiny B}}}}v_{pq,rs}D_{{\textit{\tiny B}},rs}\biggl{]}\hat{E}_{pq}\ \text{.}$$

$${\hat{V}}^t$$

$$=\sum _{\pm \omega }\left({\hat{V}}_{\text{A},\alpha }^{\omega }+{\hat{V}}_{\text{B},\alpha }^{\omega }\right)F_{\alpha }^{\omega }{e}^{-i\omega t}\text{,}$$

$$\{Q_{\text{AB}}{\}}_T$$

$$+\frac{1}{2}\sum _{{\omega }_1,{\omega }_2}⟨⟨{\hat{V}}_{\alpha }^{{\omega }_1};{\hat{V}}_{\beta }^{{\omega }_2}⟩⟩{\delta }_{{\omega }_1+{\omega }_2}{F}_{\alpha }^{{\omega }_1}{F}_{\beta }^{{\omega }_2}$$

$$+\frac{1}{6}\sum _{{\omega }_1,{\omega }_2}⟨⟨{\hat{V}}_{\alpha }^{{\omega }_1};{\hat{V}}_{\beta }^{{\omega }_2},{\hat{V}}_{\gamma }^{{\omega }_2}⟩⟩{\delta }_{{\omega }_1+{\omega }_2+{\omega }_3}{F}_{\alpha }^{{\omega }_1}{F}_{\beta }^{{\omega }_2}{F}_{\gamma }^{{\omega }_2}+\cdots \text{,}$$

$$|\widetilde{0_{\text{A}}{0}_{\text{B}}}⟩$$

$$=e^{i{\hat{\Lambda }}_{\text{A}}(t)}|0_{\text{A}}⟩\otimes e^{i{\hat{\Lambda }}_{\text{B}}(t)}|0_{\text{B}}⟩;[{\hat{\Lambda }}_{\text{A}},{\hat{\Lambda }}_{\text{B}}]=0\text{.}$$

$${\hat{\Lambda }}_{\text{A}}(t)=\sum _{i>0}\left({\lambda }_{i_{\text{A}}}(t){\hat{q}}_{i_{\text{A}}}^{†}+{\lambda }_{i_{\text{A}}}^{\ast }(t){\hat{q}}_{i_{\text{A}}}\right)={\mathbf{Q}}_{\text{A}}^{†}{\mathbf{\Lambda }}_{\text{A}}\text{,}$$

$${\mathbf{Q}}_{\text{A}}^{†}=\left[\begin{array}{cc}{\mathbf{q}}_{\text{A}}^{†}& {\mathbf{q}}_{\text{A}}\end{array}\right];{\mathbf{\Lambda }}_{\text{A}}={\left[\begin{array}{cc}{\mathbf{\lambda }}_{\text{A}}(t)& {\mathbf{\lambda }}_{\text{A}}(t{)}^{\ast }\end{array}\right]}^T\text{,}$$

$$Q_{\text{AB}}(t)=⟨\widetilde{0_{\text{A}}{0}_{\text{B}}}|\left(\hat{\mathcal{H}}+{\hat{V}}^t-i\frac{\partial }{\partial t}\right)|\widetilde{0_{\text{A}}{0}_{\text{B}}}⟩\text{.}$$

$$⟨⟨{\hat{V}}_{\alpha }^{-\omega };{\hat{V}}_{\beta }^{\omega }⟩⟩={\frac{{\text{d}}^2\{Q_{\text{AB}}{\}}_T}{\text{d}{F}_{\alpha }^{-\omega }\text{d}{F}_{\beta }^{\omega }}|}_{\mathbf{F}=0}=-{\mathbf{V}}_{\alpha }^{\omega †}({\mathbf{E}}^{[2]}-\omega {\mathbf{S}}^{[2]}{)}^{-1}{\mathbf{V}}_{\beta }^{\omega }\text{.}$$

$$-i{\mathbf{V}}_{\beta }^{\omega }={\frac{{\partial }^2\{Q_{\text{AB}}{\}}_T}{\partial F_{\beta }^{\omega }\partial {{\mathbf{\Lambda }}^{\omega }}^{\ast }}|}_{\mathbf{F}=0}=$$

$${\mathbf{E}}^{[2]}-\omega {\mathbf{S}}^{[2]}=$$

$${\mathbf{E}}_{\text{I}}^{[2]}$$

$$\displaystyle A^{{\textit{\tiny A}}}_{ij}=\$$

$$\displaystyle=\$$

$$\displaystyle B^{{\textit{\tiny A}}}_{ij}=\$$

$$\displaystyle\mathrm{g}_{\alpha,i}^{\textit{\tiny A}}=\$$

$${\mathbf{E}}_{\text{AB}}^{[2]}=\left[\begin{array}{cc}\mathbf{\Gamma }& \mathbf{\Theta }\\ {\mathbf{\Theta }}^{\ast }& {\mathbf{\Gamma }}^{\ast }\end{array}\right],$$

$$\displaystyle\Gamma_{ij}=\$$

$$\displaystyle\Theta_{ij}=\$$

$${\mathbf{E}}^{[2]}\mathbf{X}$$

$${\mathbf{X}}^{†}{\mathbf{S}}^{[2]}\mathbf{X}$$

$$\mathbf{X}=\left[\begin{array}{cc}{\mathbf{X}}^{\text{A}}& {\mathbf{X}}^{\text{AB}}\\ {\mathbf{X}}^{\text{BA}}& {\mathbf{X}}^{\text{B}}\end{array}\right],$$

Full text

polarizable embedding, response theory, localized electronic transition, subsystem decompositions, effective external field effect

I Introduction

The many types of spectroscopy that have been developed over the years (some of which are still under development) serve as an indispensable tool to gain fundamental understanding of the structure and dynamics of molecular systems.Lindon, Tranter, and Holmes (2010) Equally indispensable are the large number of theoretical methods that nowadays allow for accurate in silico simulations of said systems and spectroscopies, and which can provide complementary information that is inaccessible in the experiment. When designing calculations of molecular properties, it is often necessary to consider the influence of the surrounding environment to obtain adequate accuracy. Given the typical steep computational scaling of quantum-chemical methods with respect to the number of basis functions, the description of large molecular systems is out of question using conventional algorithms.

One way to partially address the challenges of treating large systems is to exploit localization techniques and simplified treatments of long-range interactions to achieve reduced, ultimately linear scaling.Kristensen et al. (2012); Ziółkowski et al. (2010); Ochsenfeld, Kussmann, and Lambrecht (2007); Sherrill (2010) While linear scaling formulations, aided by increasing computer power, significantly extend the applicability range of quantum chemistry, they are generally too costly for routine applications, especially if the system and property of interest put special demands on method, basis set and conformational sampling. For instance, even in cases where (time-dependent) density functional theory ((TD)DFT) is still feasible, it may suffer from electron-transfer overstabilization in the ground state Jakobsen, Kristensen, and Jensen (2013); Olsen et al. (2015a) and artificial low-lying charge-transfer transitionsDreuw, Weisman, and Head-Gordon (2003); Gritsenko and Baerends (2004); Dreuw and Head-Gordon (2004) due to the self-interaction errors associated with approximate standard exchange–correlation functionals. Another complication is that the delocalized picture obtained from such brute-force calculations is not straightforward to interpret in terms of chemical concepts and generally requires additional decomposition analyses of the wave function and properties into local components.

An alternative route is to adopt subsystem approaches in which the system is divided into chemically well-defined constituents which may be treated as separate entities perturbed by the other subsystems, and as such they offer a natural decomposition into environmental effects. Several methods belong to this category and one may generally distinguish two conceptually different approaches: (i) strict subsystem methods that treat the subsystems consistently on an equal footing and which upon recombination yield the properties of the entire system,Gordon et al. (2012); Jacob and Neugebauer (2014) and (ii) the so-called embedding approaches in which an accurate quantum-mechanical description is intended for only a smaller part of the system, while the remaining parts (the environment) and its effects on the former are described with more efficient, approximate methods. A recent review of subsystem and embedding approaches can be found in Ref. 13. Apart from the obvious computational cost argument, embedding approaches are motivated by the fact that the chemical identity of subsystems is largely intact in molecular complexes and that many electronic transitions are localized in nature. Prototypical examples are solute–solvent systems and chromophore–protein complexes. Frozen-density embeddingWesolowski and Warshel (1993) (FDE) provides in principle a full quantum-mechanical description of the system based on a subsystem formulation of DFT and can be applied in both a self-consistent and an embedding mode thereby bridging (i) and (ii). In principle, it is a formally exact framework although its quality with respect to supermolecular DFT in practical calculations is limited by approximations in the employed nonadditive parts of the exchange–correlation and kinetic energy functionals.

More efficient, though more approximate embedding approaches are quantum–classical models, such as hybrid quantum mechanics/molecular mechanics methods, which use a discrete but classical representation of the environment.Warshel and Levitt (1976); Senn and Thiel (2009); van der Kamp and Mulholland (2013) These approaches are classified in a hierarchy according to the extent of the coupling between the quantum and classical subsystems. In electrostatic embedding, the environment is represented as a purely external electrostatic potential perturbing the density of the quantum subsystem, while a more realistic description incorporating the mutual polarization effects between the quantum and classical subsystems is offered by polarizable embedding schemes. A variety of embedding schemes to model environment polarization have been proposed and applied, such as induced dipoles and fluctuating charges.Lopes, Roux, and MacKerell Jr (2009) Inclusion of the explicit environmental response is particularly important for processes where substantial rearrangements may occur, e.g., upon electronic excitation.Sneskov et al. (2011) In the FDE scheme, state-specific polarization of the environment can be included in the self-consistent mode of the formalism by interchanging the roles of the subsystems in so-called freeze-and-thaw iterations.Wesolowski and Weber (1996) In polarizable embedding, this is replaced by a self-consistent determination of the model-specific embedding parameters describing the polarization of the environment.

Both FDECasida and Wesołowski (2004); Höfener, Gomes, and Visscher (2012); Neugebauer (2010, 2009, 2007) and polarizable embedding approachesOlsen, Aidas, and Kongsted (2010); Lipparini, Cappelli, and Barone (2012); Yoo et al. (2008); List, Olsen, and Kongsted (2016); Jensen, van Duijnen, and Snijders (2003); Morton and Jensen (2011); Curutchet et al. (2009) have been generalized to a response formalism to allow for the calculation of response and transition properties of molecules embedded in large molecular complexes. The computational cost associated with the explicit coupling of the subsystem excitation manifolds in a fully coupled FDE scheme generally hinders the inclusion of the dynamical response of the entire environment in large complexes,Neugebauer (2009) but allows in a truncated form to describe a few strongly coupled local excitations as relevant for chromophoric aggregates.Neugebauer (2007) The classical treatment in polarizable embedding schemes, on the other hand, offers an efficient inclusion of the dynamic environmental response, provided the perturbing field is nonresonant with respect to local excitations in the environment.Jacob et al. (2006) This restriction has, however, been lifted by introducing phenomenological excited-state lifetimes in the response formalism, see Refs. 31; 34.

The objective of the present work is to formulate a rigorous derivation of linear response theory within polarizable embedding starting from a quantum-mechanical treatment of the entire system. While various polarizable embedding schemes differ in the specific representation of the environment, the underlying mathematical structure and physical content of their working equations are the same.List et al. (2013) For the present aim, we will focus on the explicit expressions for the polarizable embedding (PE) modelOlsen, Aidas, and Kongsted (2010); Olsen and Kongsted (2011); List, Olsen, and Kongsted (2016) in which the environment is described in terms of a distributed multipole representation, thus providing a well-defined link to the environment charge density. We note that the main conclusions will also be valid for, e.g., polarizable density embedding,Olsen et al. (2015b) which goes beyond pure polarizable embedding. While the formally exact FDE framework has previously been used as a common theoretical platform for discussing embedding models,Gomes and Jacob (2012) our derivation will be based on an explicit parameterization of the wave function of the combined system. The origin and inherent limitations of polarizable embedding become particularly transparent in this framework since approximations are introduced by imposing restrictions directly on the corresponding wave function parameterization, as we shall see.

The implications of the choice of parameterization are pertinent to the extension to response properties: the equivalence between response theory and state-specific approaches is lost when employing a nonlinear parameterization,Olsen and Jørgensen (1995) as is the case for polarizable embedding, such that different physical pictures emerge in the two formalisms. On the basis of a four-state model, it was shown in Ref. 39 that the exact second-order excitation energy for a system of two interacting subsystems contains three contributions in addition to the zeroth-order term describing the excitation within a frozen environment: (i) a term describing the differential polarization in the environment upon excitation in the embedded molecule, i.e., a classical induction effect, (ii) the difference in dispersion interaction between the environment and the embedded molecule in its ground state and excited state as given by the Casimir–Polder formula, and (iii) a term that describes the coupling between the excited state of interest with all lower-lying states of the embedded molecule, i.e., de-excitations, in terms of a screened interaction mediated by the environment through its dynamic polarizability (evaluated at the excitation frequency). For a generalization beyond a few-state model, see Refs. 40; 41. Regardless the debated interpretation of this term as a dispersion effect,Corni et al. (2005); Rösch and Zerner (1994); McRae (1957); Buhmann (2012) a nonresonant excitonic coupling between subsystemsLunkenheimer and Köhn (2012); Schwabe (2016) or a resonant coupling,Daday et al. (2015) it is a manifestation of the mutual correlation of the electrons in the subsystems. Employing a direct-product ansatz in a linear response framework recovers only the third term, although it should be noted that the correspondence is incomplete for all but the lowest excitation since de-excitations to lower-lying states different from the ground state are missing in the response treatment based on a ground-state reference.List (2015) On the other hand, only the state-specific induction term (i) is recovered in a state-specific formalism.

Whereas previous theoretical analyses have primarily focused on excitation energies,Gomes and Jacob (2012); Schwabe (2016); Corni et al. (2005) we shall here be concerned with the general linear response function, and not only its poles. Based on the linear response function of the combined system, we will derive two different subsystem decompositions, referred to as symmetric and nonsymmetric, respectively. These provide a direct link to the standard response formulations of polarizable embedding based on the classical description of the energy contribution of the interaction with the environment.Jensen, Swart, and van Duijnen (2005); List, Olsen, and Kongsted (2016) In addition to justifying the form of the environmental effects in polarizable embedding schemes, the present analysis also sheds light on the basic features of coupled systems and potential complications in separating out subsystem contributions to molecular properties. As discussed by Pavanello in the subsystem DFT framework, transitions in composite systems are not strictly localized in the sense that the pole structure of the combined system is inherited by the individual subsystem contributions to the linear response function.Pavanello (2013) We will extend that analysis by considering the consequences on the evaluation of transition strengths in an embedding framework by contrasting the symmetric and nonsymmetric decompositions. This exposition will further benefit from our consideration of the complex linear response function of the combined system that puts emphasis on conservation laws for integrated intensities and their ramifications for the decomposed forms of the response function used in practical polarizable embedding schemes.

Finally, the key features of linear response theory of a combined system subjected to a weak external electric field will be illustrated numerically by considering a six-level-model of a hydrogen-bonded dimer complex formed by water and para-nitroaniline. This simple example serves to demonstrate the results otherwise found in the theoretical analysis.

II Theory

Our presentation is intended to be self-contained and thus we provide a rather extensive theory section, starting in the first part with a brief review of the conventional quantum-mechanical direct-product ansatz for the wave function of the combined system.Ángyán (1992); McWeeny and Sutcliffe (1989); Corni et al. (2005) We then provide a rigorous formulation of linear response theory for a composite system within this theoretical framework. Our formulation allows, upon additional approximations, to recover the PE model. In practice, we exploit two different decompositions that allow for solving the response equations in effective subsystem spaces. This leads to the definition of effective quantities that provide insight into how subsystem properties are modified in the presence of another interacting subsystem. The features of the individual terms in the two decompositions are compared with a view to the evaluation of subsystem contributions to response and transition properties. Finally, to overcome some of the challenges related to the subsystem formulation in conventional linear response theory, we extend the treatment to a complex response framework. In addition, such formulation allows to demonstrate the intensity borrowing occurring in a coupled system.

In the second part of the this section, we analyze the problem of recovering the PE model and its extension to quantum-mechanical linear response theory based on the decompositions derived in the first part. Specifically, we assume a perturbation treatment of the environment subsystem as well as a multipole representation of its charge distribution. In this way, we provide a rigorous derivation of the environmental effects emerging in linear response theory within polarizable embedding.List, Olsen, and Kongsted (2016)

Atomic units (a.u.) will be used throughout and, in order to keep equations more compact, we will allow ourselves to leave out physical constants that have a numerical value of 1 a.u.

II.1 Preliminary Considerations

Consider a composite system consisting of $N$ interacting subsystems each with an integer number of electrons and with a fixed relative position. A fundamental assumption made in quantum–classical embedding models is that of nonoverlapping subsystem charge densities (zero-overlap approximation), which implies that exchange–repulsion vanishes. As a consequence, the exact wave function of the combined system can be written in the basis of direct-product states constructed from the complete antisymmetrized subsystem spaces.McWeeny and Sutcliffe (1989); Corni et al. (2005) The nonrelativistic Born–Oppenheimer Hamiltonian for the interacting system decomposes naturally as

[TABLE]

where $\hat{\mathcal{H}}_{{\textit{\tiny A}}}$ is the electronic Hamiltonian of the isolated subsystem A and $\hat{\mathcal{V}}_{{\textit{\tiny AB}}}$ the interaction operator describing the interactions between the nuclei and electrons in subsystem A with those in B. In second quantization, the interaction operator takes the form

[TABLE]

where $Z_{n}$ and $\mathbf{R}_{n}$ denote the charge and the position vector, respectively, of nucleus $n$ of the $M_{{\textit{\tiny A}}}$ nuclei in subsystem A, while $Z_{m}$, $\mathbf{R}_{m}$ and $M_{{\textit{\tiny B}}}$ are the corresponding quantities belonging to subsystem B. $\hat{E}_{pq}$ is the usual second quantization singlet one-electron excitation operator,Helgaker, Jørgensen, and Olsen (2000) and $p,q,r,s$ are used as general spatial molecular orbital indices. The affiliation of the corresponding orbitals to a given subsystem is indicated by the summation. The associated one- and two-electron integrals are defined as

[TABLE]

The first two terms in Eq. (II.1) describe the instantaneous intersubsystem Coulomb interaction between the electrons in subsystem A and the $M_{{\textit{\tiny B}}}$ nuclei in subsystem B and vice versa, whereas the third and fourth terms are the intersubsystem electron–electron and nucleus–nucleus repulsion terms, respectively. Notice that no exchange integrals between the subsystems survive within the zero-overlap approximation, and the two-electron excitation operator therefore factorizes into subsystem contributions. Occasionally, we will use an alternative representation of the interaction operator

[TABLE]

given in terms of the first-order reduced density and electrostatic potential operators:

[TABLE]

The last equality of Eq. (5) introduces a shorthand notation for the spatial integration with respect to repeated space variables.

To maximize the comparability with the effective environment models, we approximate the electronic wave function of the combined system with a single direct product of subsystem wave functions

[TABLE]

The variation of the associated energy functional according to the Rayleigh–Ritz variational principle leads to a set of coupled effective subsystem equations

[TABLE]

where the effective subsystem Hamiltonian is given by

[TABLE]

Here, $D_{{\textit{\tiny B}},rs}=\langle 0_{\textit{\tiny B}}|\hat{E}_{rs}|0_{\textit{\tiny B}}\rangle$ is an element of the first-order reduced density matrix for subsystem B. The two terms collected in the square bracket represent the classical electrostatic potential generated by the ground states of the remaining subsystems in their polarized states, i.e., in the presence of the other subsystems. In practice, any wave function ansatz may be invoked for the individual subsystems and the associated optimization conditions derived upon applying the variational principle. Equation (10) is analogous to Hartree self-consistent field theory, where the orbitals replace the subsystem wave functions. In particular, the parameterization in Eq. (8) discards direct-product states in which the subsystems are simultaneously excited, which implies that no dispersion effects are included in this approximation. However, some intersubsystem electron correlation effects are introduced when the direct-product ansatz is used in a response framework, as mentioned in the Introduction.

III Response theory of composite systems

We will now consider the electronic response of the composite system to an external optical field within the framework introduced in the previous section. In such a case, nuclear motions are neglected and only the pure electronic response will be considered. Special attention will be paid to the physical aspects of the intersubsystem interactions in the presence of the applied field, and how they influence molecular properties. For the sake of notational simplicity, we shall restrict this analysis to two subsystems (leading to an embedded subsystem A and an environment B in polarizable embedding) and work within a configuration interaction (CI) framework for the individual subsystems. The framework can easily be extended to consider the individual subsystems of the environment still assuming nonoverlapping subsystems, by decomposing the environment wave function into a product of subsystem contributions (see Eq. 8). This will be used in Sec. IV. We also note that, even if we are considering explicitly the CI parameterization for wave functions of the individual subsystems, the principles apply more generally also to other variational wave function models.

The field–matter interaction will be described within the semi-classical framework, where the incident field is treated as a classical plane wave that perturbs the molecular system. The perturbation operator describing the action of a monochromatic external field on the composite system can then be expressed in terms of Fourier components as

[TABLE]

where $F_{\alpha}^{\omega}$ are Fourier amplitudes associated with the one-electron perturbation operators $\hat{V}_{\alpha}^{\omega}$, using Greek subscripts as (possibly composite) Cartesian labels. In order to maintain Hermiticity of $\hat{V}^{t}$, we have $F_{\alpha}^{\omega}=[F_{\alpha}^{-\omega}]^{*}$ and $\hat{V}_{\alpha}^{\omega}=[\hat{V}_{\alpha}^{-\omega}]^{\dagger}$. Note that the additivity of the perturbation operator in the last equality of Eq. (11) is a consequence of the zero-overlap assumption. Here and henceforth, the Einstein summation convention is adopted for repeated Greek indices.

Our derivation will follow the quasi-energy formulation of response theory,Christiansen, Jørgensen, and Hättig (1998); Sasagane, Aiga, and Itoh (1993) which can be viewed as a time-dependent generalization of the ordinary energy-differentiation technique from time-independent perturbation theory, which in its time-averaged formulation is restricted to time-periodic perturbations. The response functions are defined as the coefficients of the time-averaged quasi-energy in a Taylor series expansion in terms of the external field strengths

[TABLE]

where $\{\cdot\}_{T}$ indicates that the time-average over one period of oscillation in the external field has been taken. The symbol $\delta_{\omega}$, not to be confused with a Dirac-delta or a discrete Kronecker-delta function, is unity when the continuous frequency variable vanishes and zero otherwise. An instructive exposition of the quasi-energy derivative approach as well as a comparison to the alternative Ehrenfest formulation can be found in Ref. 54.

III.1 The Direct-Product Approximation

The phase-isolated part of the time-dependent direct-product wave function for the composite system may be defined by an exponential unitary parameterization as

[TABLE]

The intersubsystem commutation relation above follows from the zero-overlap assumption. The time-dependent Hermitian operator $\hat{\Lambda}_{{\textit{\tiny A}}}(t)$ for subsystem A is parameterized in terms of a set of time-dependent amplitudes $\{\lambda_{i_{\textit{\tiny A}}}\}$ and takes the formOlsen and Jørgensen (1985)

[TABLE]

where an identity operator for subsystem B is implied. The state-transfer operators $\hat{q}_{i_{{\textit{\tiny A}}}}^{\dagger}{=}|i_{{\textit{\tiny A}}}\rangle\langle 0_{{\textit{\tiny A}}}|$ and their adjoints are built from a set of orthonormalized states $\{|i_{{\textit{\tiny A}}}\rangle\}$ that spans the orthogonal complement space of the reference state $|0_{{\textit{\tiny A}}}\rangle$, the latter satisfying the variational condition given by Eq. (9). Before proceeding, we briefly comment on the employed parameterization: (i) The direct-product parameterization in Eq. (III.1) is nonlinear: that is, it produces states outside the excitation manifold, defined by the linear action of the operator $\hat{\Lambda}_{{\textit{\tiny A}}}(t)+\hat{\Lambda}_{{\textit{\tiny B}}}(t)$. In particular, the exponential parameterization contains states in which both subsystems are excited simultaneously due to products of subsystem excitations, despite the absence of such transitions in the state-transfer operators. This is analogous to the nonlinear exponential parameterization of a Hartree–Fock or Kohn–Sham determinant that is based on a generator of single-electron excitations but yet encompasses multi-electron excited determinants. As mentioned in the Introduction, a direct consequence of the use of a nonlinear parameterization is that the properties derived from a response framework will differ from those based on the state-specific formulation. (ii) The lack of a $\hat{\Lambda}_{{\textit{\tiny AB}}}(t)$ operator in the exponent of the first equality in Eq. (III.1) and thus a direct coupling between direct-product states of the type $\langle i_{\textit{\tiny A}}0_{\textit{\tiny B}}|$ and $|i_{\textit{\tiny A}}j_{\textit{\tiny B}}\rangle$ as well as $\langle i_{\textit{\tiny A}}0_{\textit{\tiny B}}|$ and $|k_{\textit{\tiny A}}j_{\textit{\tiny B}}\rangle$ is the origin of the neglect of state-specific relaxation and London dispersion between the subsystems. (iii) The expansion of the time-dependent wave function in the basis of subsystem states means that only intrasubsystem transitions are included, while intersubsystem transitions are excluded. Therefore, a main approximation in the above parameterization is the exclusion of charge-transfer transitions between subsystems in line with our initial restriction on a fixed number of electrons in a given subsystem.

Having settled on an approximate parameterization of the phase-isolated wave function, we can now construct the time-dependent quasi-energy of the composite system

[TABLE]

The time-dependent amplitudes are expanded in orders of the perturbation and determined by imposing the variational principle for the time-averaged quasi-energy to be satisfied at the various orders. Following Eq. (III), explicit expressions for response functions of the combined system can then be obtained as perturbation-strength derivatives of the time-averaged quasi-energy in Eq. (15), evaluated at zero field strengths. The linear response function becomes

[TABLE]

Ordering the Fourier components of the configuration amplitudes in the operators in Eq. (III.1) according to: $\mathbf{\Lambda}^{\omega}=(\bm{\lambda}^{\omega}_{{\textit{\tiny A}}},\,\bm{\lambda}^{\omega\,*}_{{\textit{\tiny A}}},\,\bm{\lambda}^{\omega}_{{\textit{\tiny B}}},\,{\bm{\lambda}^{\omega}_{{\textit{\tiny B}}}}^{*})^{T}$, leads to the following intra- and intersubsystem blocked forms of the vectors and matrices:

[TABLE]

[TABLE]

The diagonal blocks of the electronic Hessian and overlap matrices, arising upon differentiation of the quasi-energy with respect to the wave function parameters belonging to the same subsystem (${\textit{I}}={\textit{A}},{\textit{B}}$), take the same overall form as for the isolated subsystems

[TABLE]

where the symbols $\mathbf{0}$ and $\mathbf{1}$ are used to denote appropriately sized null and identity matrices. As exemplified for subsystem A, the elements of the subsystem blocks take the form

[TABLE]

In addition to the isolated subsystem term, the electronic subsystem Hessian in Eq. (19) incorporates a contribution from the intersubsystem coupling. In particular, it describes the coupling between intrasubsystem excitations under the influence of the electrostatic potential produced by the electronic ground state of the other subsystem. Note that the vanishing $\mathbf{B}^{{\textit{\tiny I}}}$ blocks result from our choice of a CI parameterization within the subsystems. However, we choose to keep these blocks in our above illustration of the structure of the electronic Hessian of the subsystems as they could be nonzero for other wave function models. The generally rectangular off-diagonal block $\mathbf{E}_{{\textit{\tiny AB}}}^{[2]}$ and its conjugate transpose ${\mathbf{E}_{{\textit{\tiny BA}}}^{[2]}}$ in Eq. (18) describe the intersubsystem coupling. Using the same ordering as in Eq. (19), the structure of this block can be written as

[TABLE]

with elements given by

[TABLE]

The off-diagonal blocks couple excitations in one subsystem with those in the other, and as follows from these expressions, the coupling is described through a Coulombic interaction of transition densities in the two subsystems. Note that the off-diagonal blocks $\mathbf{S}_{{\textit{\tiny AB}}}^{[2]}$ and $\mathbf{S}_{{\textit{\tiny BA}}}^{[2]}$ vanish due to the zero-overlap assumption.

By analogy to exact-state theory, a pole and residue analysis of the linear response function in Eq. (16) determines expressions for excitation energies and transition strengths for ground- to excited-state transitions in the combined system for the specific choice of parameterization in Eq. (III.1). Accordingly, excitation energies are found as eigenvalues of the generalized eigenvalue equation involving the electronic Hessian and the metric in terms of the overlap matrix:Olsen and Jørgensen (1985)

[TABLE]

We recall the well-known feature of this generalized eigenvalue problem where the eigenvalues come in pairs $\pm\omega_{n}$ and the associated eigenvectors are related,Olsen and Jørgensen (1985); List et al. (2014) and we shall use positive and negative indices ($\omega_{-n}=-\omega_{n}$) to denote paired solutions. Here, $\mathbf{X}$ is the matrix of eigenvectors satisfying the orthonormality relation

[TABLE]

and $\bm{\Omega}$ is the diagonal matrix containing the associated eigenvalues $\pm\omega_{n}$.

For later reference, we introduce the partitioned form of the eigenvector matrix according to the blocked structure in Eq. (18)

[TABLE]

where the off-diagonal blocks describe the degree of delocalization of the excitations in the combined system. The transition strength associated with excitation $n$ in the composite system can be obtained from the residue of the linear response function as

[TABLE]

and used to compute the dimensionless oscillator strength

[TABLE]

Having outlined the response formalism for the combined system, the main objective in the two following subsections is to obtain expressions for subsystem contributions to properties of the combined system by solving effective equations within the subsystem spaces. As will be shown in Section IV, this alternative approach to the solution of the full system provides a direct link to polarizable embedding.

III.1.1 Subsystem decomposition: electronic response properties

To decompose the linear response function of the combined system into subsystem contributions, we use that the inverse of a blocked matrix with nonsingular square diagonal blocks ($\mathbf{U}$ and $\mathbf{Z}$) may be written as

[TABLE]

which corresponds to using the Löwdin partitioning technique.Löwdin (1962a, 1963, b) Using this identity, we can rewrite the matrix resolvent in Eq. (16) in a partitioned form and obtain the following subsystem decomposition of the linear response function:

[TABLE]

Since the subsystems are treated on the same footing in this representation, Eq. (III.1.1) will in the following be referred to as the symmetric subsystem decomposition (SD). In contrast to Eq. (16), the dimensions in the above expression have been reduced to those of the excitation manifolds of the individual subsystems (e.g., $\text{dim}_{\textit{\tiny A}}{\times}\text{dim}_{\textit{\tiny A}}$) instead of that of the full system $(\text{dim}_{\textit{\tiny A}}{+}\text{dim}_{\textit{\tiny B}}){\times}(\text{dim}_{\textit{\tiny A}}{+}\text{dim}_{\textit{\tiny B}})$. As indicated by tildes, the response vectors $\widetilde{\mathbf{N}}^{\omega}$ in Eq. (III.1.1) are modified quantities that satisfy the following set of effective linear response equations

[TABLE]

In addition to the implicit modifications through the polarization of the reference state vectors (changes induced by interactions within the subsystem through the interaction part of the diagonal blocks, e.g., $\mathbf{E}_{{\textit{\tiny A}}}^{[2]}$), the presence of the other subsystem manifests in explicit contributions to the electronic Hessian and property gradients. In Eqs. (III.1.1) and (32), this has been written compactly by introducing effective Hessians and property gradients defined, here for subsystem A, as

[TABLE]

These two expressions are key of the present work, since they show how the properties of subsystem A are affected by the presence of another interacting subsystem. As follows, the coupling interactions between subsystems are governed by three different kinds of mechanisms. The response vector for subsystem A describes changes induced by: (i) the indirect coupling through the modification of the pure intrasubsystem term $\mathbf{E}^{[2]}_{\textit{\tiny A}}$ to include the ground-state electrostatic potential of subsystem B. (ii) The explicit coupling to subsystem B through the frequency-dependent term in the electronic Hessian. A similar result was obtained in the work by Neugebauer in which a subsystem partitioning of the eigenvalue problem within the linear response generalization of subsystem DFT is discussed.Neugebauer (2007) (iii) The direct interaction between the applied external field and subsystem B through the modified property gradient. In particular, we recognize the matrix resolvent of the second term in Eqs. (33) and (34) as a matrix representation of the frequency-dependent linear polarizability of the polarized reference state of subsystem B (i.e., self-consistently polarized by subsystem A), evaluated at the optical frequency $\omega$.Hsu et al. (2001) This can be seen by rewriting the intersubsystem coupling blocks according to

[TABLE]

and similarly

[TABLE]

where the vectors are defined in the same way as the property gradient in Eq. (19) but with the perturbation operator being replaced by either the electrostatic potential operator or the density operator. The off-diagonal block of the electronic Hessian can then be written as

[TABLE]

Substituting this expression into the last term of Eq. (33), we obtain

[TABLE]

where the factor in square brackets is identified with the frequency-dependent generalized linear polarizability of the polarized ground state of subsystem B: $C_{\mathbf{r},\mathbf{r}^{\prime}}^{{\textit{\tiny B}}}(\omega)=\langle\langle\hat{\rho}_{\textit{\tiny B}}(\mathbf{r});\hat{\rho}_{\textit{\tiny B}}(\mathbf{r}^{\prime})\rangle\rangle_{\omega}$, which, upon assuming the dipole approximation for subsystem B, reduces to the electric dipole–dipole polarizability tensor $\bm{\alpha}^{{\textit{\tiny B}}}(\omega)$.

The property gradients for subsystem A describe the electrostatic potential generated by a transition density. The last term in Eq. (33) can thus be interpreted as the linear response of subsystem B induced by the electrostatic potential due to a transition density of subsystem A, which in turn produces an electrostatic potential acting on A. In the same way, the last term of the effective property gradient in Eq. (34) can be interpreted as the electrostatic potential acting on subsystem A due to the linear polarization induced in subsystem B by its direct interactions with the external field. In other words, subsystem B acts as a source of a field that gives rise to an effective field strength at the location (but devoid) of subsystem A different from the external field, as represented by $\widetilde{\mathbf{V}}_{{\textit{\tiny A}},\beta}^{\omega}$ and ${\mathbf{V}}_{{\textit{\tiny A}},\beta}^{\omega}$, respectively. Note that because of the tensorial nature of the environment polarizability, the external field can be screened differently in different directions. It is also worth underlining the relative distance dependency of these two interaction mechanisms. Assuming a dipole approximation for the interaction operator, it follows that the former effect decays quickly with increasing intermolecular distance $r$, displaying a $r^{-6}$ behavior, whereas the latter effect is comparably long range, dropping of as $r^{-3}$. This difference in interaction range is important to consider in practical calculations in order to get a balanced description of the various contributing terms when dealing with finite-sized systems.Rinaldi, Morton, and Jensen (2013)

As we shall see later, the second terms of Eqs. (33) and (34) reduce exactly to the form of the so-called dynamic reaction field and the effective external field (EEF) effects, the latter also referred to as the local field effect, appearing in polarizable embedding.Jensen, Swart, and van Duijnen (2005); List, Olsen, and Kongsted (2016) The effective external field effect is of the same origin as the cavity-field effect in polarizable continuum models leading to effective properties,Böttcher and Bordewijk (1973); Wortmann and Bishop (1998); Cammi, Mennucci, and Tomasi (1998); Pipolo, Corni, and Cammi (2014) but is defined with respect to the external probing field rather than the Maxwell field within the dielectric.List, Jensen, and Kongsted (2016) It should also be noted that the same terms describe the so-called image- and near-field effects in the context of treating molecules adsorbed to metal nanoparticles.Morton and Jensen (2011); Payton et al. (2012)

Finally, we note that neither of the subsystem contributions to the linear response function in Eq. (III.1.1) are symmetric with respect to the left and right property gradients. For that reason, there is no guarantee that the effective subsystem polarizability tensors are symmetric, as discussed before,Neugebauer (2009) or that the diagonal elements are positive in the static limit. This is in contrast to the polarizability tensor of the combined system, which is symmetric, as it should be. Physically, the nonsymmetric form of the individual terms can be viewed as describing the linear response of a property associated with the operator $\hat{V}_{{\textit{\tiny A}},\alpha}^{-\omega}$ to the actual perturbing field acting on subsystem A in the presence of subsystem B.

III.1.2 Subsystem decomposition: electronic transition properties

In this section, we perform an analysis of the transition properties of the composite system as relevant for resonant external fields. As pointed out in the Introduction, the lowest electronic transitions are often localized in nature and can often be attributed predominantly to a particular subsystem, i.e., the excitation vectors have the dominant contribution in the excitation manifold of one of the subsystems. The condition under which such excitations occur follows from perturbation theory: the magnitude of the coupling strengths of the excited state of interest with electronic states in the other subsystem should be small compared with their energetic separation. This ratio dictates the magnitude difference between the diagonal and off-diagonal blocks in Eq. (27) and thus the extent of the delocalization of a given electronic excitation in the combined system.

Rather than solving the full generalized eigenvalue equation in Eq. (25), the excitation energies of the combined system may be determined from the effective response equations in Eq. (32) upon zeroing the right-hand sides. Folding the effects of subsystem B into the equation for A yields the following pseudo generalized eigenvalue equation

[TABLE]

which has the dimension of the excitation manifold of subsystem A. Provided there are no degenerate states located in the B part, the excitation energies derived from the reduced system in Eq. (39) are identical to those associated with the subsystem A part of the parent system in Eq. (25), i.e., the full system. This reformulation thus turns the generalized eigenvalue problem for the full system into a dressed subsystem problem, as has been shown before in the subsystem TDDFT framework.Neugebauer (2007) In other words, this form allows us to compute the poles associated with the A (B) block of Eq. (25) without having to consider the full problem with the dimension of both subsystems. This is particularly attractive for determining localized transitions as relevant from the perspective of embedding calculations. The nonlinearity of Eq. (39) introduced by the frequency-dependent effective Hessian, however, requires knowledge of the solutions beforehand or that the problem is solved iteratively one excitation at a time, including the construction of a new Hessian for each eigenvalue in question. As a consequence, the subsystem A components will not be orthogonal to each other (only the full eigenvectors will). If only the subsystem term $\mathbf{E}_{{\textit{\tiny A}}}^{[2]}$ is included in Eq. (39), the result will be denoted by ”FP” for frozen ground-state polarization—implying that subsystem B is not allowed to respond to the density changes in A upon excitation (that is, neglecting intersubsystem couplings). As discussed in the previous section, the second term in Eq. (33) couples an excitation in subsystem A with those in B by the interaction of the associated transition densities (see Eq. (37)).

Let us now return to the evaluation of the transition strengths associated with the excitations in the composite system, but now taking a decomposed form of the linear response function as a starting point. As has been shown by Pavanello,Pavanello (2013) each subsystem contribution to the response function (i.e., symmetric decomposition in Eq. (III.1.1)) displays poles at all excitations in the combined system. This complication reflects the delocalization necessarily present in coupled quantum systems, which is not compliant with our heuristic view on local excitations. As a consequence, transition strengths of electronic transitions ”localized” in subsystem A cannot be identified as the residues of the subsystem A contribution to the response function (the first term in Eq. (III.1.1)).

To facilitate the identification of residues, we will explore an alternative partitioning of the matrix resolvent in Eq. (16), which can be obtained by applying the Woodbury matrix identity (for $\mathbf{U}$ and $\mathbf{Z}$ nonsingular square matrices)

[TABLE]

to the third and fourth blocks of Eq. (30) (alternatively, the first and second blocks to get the residues related to the B part). Hereby, the linear response function of the combined system can be written as

[TABLE]

in which the subsystems, in contrast to Eq. (III.1.1), are treated on an unequal footing. Equation (42) will therefore be referred to as the nonsymmetric decomposition (NSD). The same approach has been used in the original formulations of second-order approximate polarization propagator (SOPPA) theoryNielsen, Jørgensen, and Oddershede (1980) in which the effective quantities describe the doubles correction to the particle-hole spectrum, while the pure double excitations are given by the second term. As will be discussed in Sec. IV.2, the expressions for polarizable embedding in a linear response framework follows, in a similar spirit as SOPPA and algebraic diagrammatic construction,Dreuw and Wormit (2015) from a perturbation analysis of the individual matrices appearing in Eq. (42). Note that the second term originates from the first term on the right-hand side of Eq. (41) and corresponds to the linear response function for the ground-state polarized subsystem B. An important outcome of this alternative representation of the response function of the combined system is that the poles of the A-dominated excitations (all excitations of the composite system) are contained entirely in the first term. However, as a consequence of the appearance of the second term, the partitioning in Eq. (42) introduces additional poles in the individual terms at the transition frequencies of the ground-state polarized subsystem B. Hence, the correct pole structure for the composite system is only recovered upon taking the sum of the two terms. Nevertheless, the first term in Eq. (42) can be used to identify the correct decomposed expression for the residues of the excitations mainly in subsystem A, since the second term does not affect excitation energies and transition moments but is needed only in the calculation of the response function.

Until now the normalization of the effective eigenvectors obtained from Eq. (39) has been of no concern since it does not affect excitation energies. However, it is necessary to consider a renormalization of the effective eigenvectors before the transition strengths can be evaluated in the decomposed formulation. In particular, since the eigenvector of the full equation, corresponding to a transition mainly localized in subsystem A, is normalized to $\pm 1$ according to

[TABLE]

the subsystem A component (first term) must have a norm less than unity. The A and B components of the eigenvector are related through Eq. (25), which can be used to rewrite the normalization condition in terms of the A component. This leads to the following renormalization factor for the effective eigenvector for the given pole $\omega_{n}$

[TABLE]

The transition strength between the ground state and an excited state mainly located in subsystem A can then be written, in the partitioned form, as

[TABLE]

What we have achieved up to this point is to have recast Eqs. (16) and (25) into dressed subsystem expressions that separate out subsystem contributions to molecular and transition properties of the combined system. As will become clear from the final steps taken in Sec. IV, this provides a justification of the various environmental effects appearing in polarizable embedding. The theoretical analysis further allowed to identify what decomposed form of the response function to use for a specific purpose. We underline that such clarification does not follow from a derivation anticipating the classical description of the environment in the first place, where the Ehrenfest and quasi-energy derivative formulations of response theory lead to different response functions.List, Jensen, and Kongsted (2016)

Although the subsystem decompositions based on standard response theory are illustrative and provide insight into the mechanisms governing the interaction between subsystems, their practical application for transition properties is limited for several reasons: (i) the solution of the pseudo generalized eigenvalue problem in Eq. (39) requires an iterative scheme, (ii) a new effective subsystem A Hessian has to be constructed for each ”eigenvalue” of interest, and (iii) the problem is ill-defined if poles in subsystem B are too close to the one being solved for. Another important point is that the normalization factor in Eq. (44), necessary for the evaluation of absorption intensities, cannot be straightforwardly converted into an effective environment analog as needed when turning to embedding models. As will be discussed next, these complications can be avoided by considering the combined system within a complex response theory framework in which absorption properties can be computed without having to resolve the individual excitations, which is also advantageous for systems in which the density of electronic states is high. In addition, such framework clearly illustrates the intensity borrowing that occurs between interacting subsystems.

III.2 Subsystem decomposition in a complex response framework

In complex response theory, effects of radiative and nonradiative relaxation mechanisms for the decay of the excited states are modeled in a phenomenological manner by assigning finite lifetimes ($\tau_{n}$) to the excited states. This leads to complex-valued response functions that are well-defined across the entire frequency range and thus provides resonant-convergent properties. The complex linear response function takes the following formNorman (2011); Norman et al. (2005, 2001)

[TABLE]

In practice, it is customary to adopt a common lifetime, and hence damping parameter $\gamma=(2\tau)^{-1}$, for all excited states such that the so-called relaxation matrix becomes $\bm{\gamma}=\gamma\mathbf{1}$.

Similar to the conventional response framework, the complex linear response function may be expressed in alternative subsystem decomposed forms. Applying Eq. (30) yields its symmetric subsystem decomposition

[TABLE]

The effective subsystem vector and matrix quantities are now complex, here given for subsystem A:

[TABLE]

Note that even in the usual case of a real (or purely imaginary) external perturbation, the property gradient will be complex as a result of the damped response of the other subsystem (the second term of Eq. (49)). By further rewriting according to Eqs. (40) and (41), we obtain the nonsymmetric decomposition

[TABLE]

It should be noted that the conjugate transpose of the effective quantities are here assumed, as indicated by the prime, to act only on vectors and matrices but without changing the sign in front of the damping parameter.

In actual calculations, the value of the response function is determined by solving the complex linear response equation. For the subsystem decompositions, this implies solving the complex analogs of Eq. (32). To be amenable to practical implementation, these may be expressed as a coupled set of linear equations for the real and imaginary componentsKauczor, Jørgensen, and Norman (2011) (indicated by superscripts $R$ and $I$, respectively). By using Eq. (III.1.1) and assuming real wave functions, we obtain the following expression for the effective response equations for subsystem A

[TABLE]

where $C_{\mathbf{r},\beta}^{{\textit{\tiny B}}}(\omega)=\langle\langle\hat{\rho}^{{\textit{\tiny B}}}(\mathbf{r});\hat{V}^{\omega}_{{\textit{\tiny B}},\beta}\rangle\rangle_{\omega}$. As seen, the coupling between the real and imaginary components of the effective response vector is, in addition to the vacuum contribution from the damping within subsystem A, mediated by the imaginary part of the generalized linear polarizability for subsystem B.

The real part of the complex electric dipole–dipole polarizability is related to the refractive index of the system and, as follows from energy-loss considerations of the perturbing electromagnetic field, the imaginary part is directly proportional to the linear absorption cross section $\sigma(\omega)$.Boyd (2003); List et al. (2015) For a sample that is isotropic with respect to the light polarization, we have

[TABLE]

where $\epsilon_{0}$ is the vacuum permittivity and $c_{0}$ is the speed of light in vacuum.

The physical significance of the imaginary component of the polarizability can alternatively be recognized by the relation between the integrated absorption cross section and the sum of oscillator strengths (see Appendix A)

[TABLE]

In the framework of exact state theory or variational approximate state theory in the complete basis set limit, the Thomas–Reiche–Kuhn sum rule further implies that the sum of oscillator strengths is equal to the number of the electrons in the system $N_{\mathrm{e}}$. Accordingly, for the combined system, we have

[TABLE]

and likewise for the subsystems

[TABLE]

Together, Eqs. (53) and (54) reflect that excitations in one subsystem can borrow intensity from transitions in the other subsystem, while the total integrated cross section is preserved. By combining Eqs. (53)–(55) and recalling that the second term of the complex linear response function in the nonsymmetric subsystem decomposition in Eq. (50) is identical to the linear response functions of subsystem B within the FP approximation, it follows that

[TABLE]

where $\sigma^{\text{NSD}}_{1}$ and $\sigma^{\text{NSD}}_{2}$ denote the contribution to the absorption cross section from the first and second term in Eq. (50), respectively. That is, if an A-dominated transition gains in intensity due to the coupling to excitations in subsystem B, then $\sigma^{\text{NSD}}_{1}(\omega)$ will take on negative values around poles in subsystem B. Consequently, in these regions of the spectrum, $\sigma^{\text{NSD}}_{1}(\omega)$ in itself cannot be associated with an absorption spectrum but it rather becomes imperative to consider the total absorption cross section $\sigma(\omega)$. This discussion is important due to its implications in the context of polarizable embedding (see next section) in which one focuses on the calculation of $\sigma^{\text{NSD}}_{1}(\omega)$ and, as we have seen, caution is called for in the interpretation of the results of such a calculation.

IV Polarizable Embedding

Having derived the expressions for the direct-product ansatz for the combined system, we will in this section detail the additional steps that lead to the definition of the PE model, following the work of Ángyán for the derivation of the effective embedding operator.Ángyán (1992) In particular, the subsystems will now be treated at different levels, where a classical description will be adopted for subsystem B. To reflect this distinction, subsystem A will in this section be referred to as the quantum region and B as the environment. Furthermore, instead of considering the environment as a whole the individual subsystems constituting the environment ${B}=\{b_{1},b_{2},\dots,b_{N-1}\}$ will be treated separately by decomposing the environment wave function into a product of subsystem contributions, still assuming nonoverlapping subsystem charge densities. We will assume that the unperturbed (i.e., vacuum) environment subsystem eigenfunctions and -energies $\{|0_{{b}}^{(0)}\rangle,|j_{{b}}^{(0)}\rangle\}$ and $\{E_{0_{b}}^{(0)},E_{j_{b}}^{(0)}\}$, for $b\in{\textit{B}}$, are known, where superscripts $(n)$ specify the order with respect to the perturbation (see below).

IV.1 Working Equations

In the PE model, we invoke a perturbation treatment of all but subsystem A and assume that the environment is only linearly responsive. This corresponds to requiring that Eq. (9) for $b\in{\textit{B}}$ is fulfilled only through first order in terms of the electrostatic potentials from the ground states of the other subsystems. Within this approximation, the interaction operator acting on subsystem A takes the form

[TABLE]

where an element of the zeroth- and first-order electronic densities of subsystem $b$ are defined as $D_{b,rs}^{(0)}=\langle 0_{b}^{(0)}|\hat{E}_{rs}|0_{b}^{(0)}\rangle$ and $D_{b,rs}^{(1)}=\langle 0_{b}^{(1)}|\hat{E}_{rs}|0_{b}^{(0)}\rangle+\langle 0_{b}^{(0)}|\hat{E}_{rs}|0_{b}^{(1)}\rangle$, respectively. The effective interaction operator acting on subsystem A consists of contributions from the permanent and induced charge distributions, $\hat{\mathcal{V}}^{\text{es}}$ and $\hat{\mathcal{V}}^{\text{ind}}$, respectively, of the environment subsystems. In terms of the first-order reduced density and electrostatic potential operators in Eqs. (6) and (7), they read

[TABLE]

where superscript ”e” signifies that only the electronic part of the operator is included. We have further introduced a shorthand notation for expectation values, e.g.,  $\langle\hat{\mathcal{V}}_{b}(\mathbf{r})\rangle_{0_{b}}^{(0)}=\langle 0_{b}^{(0)}|\hat{\mathcal{V}}_{b}(\mathbf{r})|0_{b}^{(0)}\rangle$. The $\hat{\mathcal{V}}^{\text{es}}$ operator is straightforwardly constructed from the charge densities of the unperturbed environment subsystems and contains the contributions from both the nuclei and electrons in the environment. The induction operator requires the first-order densities to be known. Using standard Rayleigh–Schrödinger perturbation theory yields the following first-order perturbation expression for the wave function of the environment subsystemsÁngyán (1992)

[TABLE]

where the expectation value involving subsystem A is over the fully polarized state within the framework of Eq. (9). As pointed out by Stone,Stone (1989) this expression is however inconsistent with the first-order perturbation analysis: through the coupling to the other environment subsystems in the third term, Eq. (IV.1) contains contributions to infinite order in the electrostatic potential generated by subsystem A. A strict first-order expression can be obtained by neglecting the many-body polarization among the environment subsystems, i.e., removing the first-order term of Eq. (IV.1).

Substituting the first-order correction to the wave function in Eq. (IV.1) into $\hat{\mathcal{V}}^{\text{ind}}$ and using the definition of the generalized static polarizability, $C_{\mathbf{r},\mathbf{r}^{\prime}}^{b,(0)}(\omega=0)$, of the zeroth-order ground states of the environment subsystems, the first-order induction operator can be written as

[TABLE]

In the PE model, the charge distributions of the environment subsystems are represented by multipole expansions rather than densities. To this end, it is expedient to use a 3-dimensional multi-index notation in which a multi-index $k=(k_{x},k_{y},k_{z})$ is an ordered list of positive integers.Saint Raymond (1991); Olsen (2012) The norm and factorial of a multi-index are defined as $|k|=k_{x}+k_{y}+k_{z}$ and $k!=k_{x}!\cdot k_{y}!\cdot k_{z}!$, respectively, the sum of two multi-indexes as $k\pm l=(k_{x}\pm l_{x},k_{y}\pm l_{y},k_{z}\pm l_{z})$, and the multi-index power of a vector as $\mathbf{r}^{k}=x^{k_{x}}\cdot y^{k_{y}}\cdot z^{k_{z}}$. Using this notation, a Taylor series expansion of the electrostatic potential operator can be written asStone (2002)

[TABLE]

where $\mathbf{R}_{\text{o}}$ is the expansion point and the summation over $k$ runs over the $\tfrac{1}{2}(|k|+1)(|k|+2)$ Cartesian components. The multipole form of the potential operator involves components of the Cartesian interaction tensors defined as derivatives of the potential operatorBuckingham (1967); Stone (2002)

[TABLE]

where the superscript multi-index notation should not be confused with the perturbation order. Note that in writing Eq. (62), we have used the symmetry properties of the interaction tensors, that is, $T_{ij}^{(k)}{=}({-}1)^{|k|}T_{ji}^{(k)}$. Applying this to the electrostatic potential operators for the environment subsystems yields

[TABLE]

where the expansion point $\mathbf{R}_{b}$ is chosen to reside within the charge density of subsystem $b$. A Cartesian component of the $|k|$’th-order multipole moment operator acting on subsystem $b$ is given by

[TABLE]

where the associated electronic multipole integral is given by

[TABLE]

The expectation value of Eq. (IV.1) corresponds to an electric multipole moment, $M^{(k)}_{b}(\mathbf{R}_{b})$, of subsystem $b$. For instance, letting $|k|=0$ gives the charge, while $|k|=1$ covers the three Cartesian components, i.e., $(1,0,0)$, $(0,1,0)$ and $(0,0,1)$, of a dipole moment. Only the lowest-order nonvanishing multipole moment of the permanent charge distribution is independent of the choice of origin (here $\mathbf{R}_{b}$). For the sake of brevity, we will suppress this explicit origin dependence in the following. Finally, substitution of Eq. (64) into $\hat{\mathcal{V}}^{\text{es}}$ yields the multipole-expanded form of the electrostatic interaction operator

[TABLE]

where the electrostatic potential integral is defined over the interaction tensors in Eq. (63) as

[TABLE]

where the index $i$ refers to an electronic coordinate. Comparing Eq. (67) to the nonexpanded form in Eq. (9), it is clear that the two-electron integrals have been replaced by one-electron integrals for subsystem A over the interaction tensors multiplied by the permanent multipole moments of the environment subsystems.

We proceed in the same way for the induction operator in Eq. (IV.1) by replacing all occurrences of the potential operator by its Taylor series expansion. We further define a component of the $|k|$’th-order induced multipole moment belonging to an environment subsystem $b$ as

[TABLE]

where we have introduced the bar notation to distinguish induced multipole moments from their permanent counterparts. By contrast to the permanent moments, there are no induced monopoles and the nuclear contributions vanish irrespective of the multipole order. This is a result of the intermediate normalization of the corrections to the wave functions of the environment subsystems. As detailed in Appendix B, the induction part of the interaction operator can be written in the multipole-expanded form as

[TABLE]

where the two alternative expressions arise by expanding Eqs. (59) and (IV.1), respectively. $P_{{b}}^{(k,l)}$ are static electronic polarizabilities of subsystem $b$, analogous to that in Eq. (III.1.1), defined as

[TABLE]

recalling that this definition employs the traced multipole moment operators. Furthermore, $\hat{F}_{{\textit{\tiny A}}}^{(k)}(\mathbf{R}_{b})$ is the $(k{-}1)$’th-order electric-field derivative, which probes the field derivative produced by subsystem A at point $\mathbf{R}_{b}$

[TABLE]

as customary defining the field ($|k|=1$) as minus the gradient of the electrostatic potential. The operator has been partitioned into nuclear and electronic contributions. For the environment subsystems, we invoke a multipole expansion of the field operator by analogy to Eq. (64). Accordingly, we can write the zeroth- and first-order electric-field derivatives of the environment subsystems in Eq. (71), in terms of static and induced multipole moments, respectively, as

[TABLE]

Note the difference between the lower limits of the two summations. Finally, by equating the right-hand sides of Eqs. (70) and (71), we obtain the equation determining the induced multipole moments

[TABLE]

where $F^{\text{tot}(l)}(\mathbf{R}_{b})$ is the total $(l{-}1)$’th-order electric-field derivative acting on subsystem $b$. As follows from the second equality, it consists of the physical electric-field contributions from the nuclei and electrons in subsystem A and the permanent multipoles of the other environment subsystems, collectively denoted $F^{(l)}(\mathbf{R}_{b})$, as well as the contribution from the remaining first-order induced multipole moments. Hence, the first term describes the mutual coupling between subsystem A with all environment subsystems, whereas the second and third terms account for the mutual polarization between the environment subsystems.

In practice, the multipole expansions in Eqs. (67) and (70) are terminated at a finite order $K_{s}$, and to improve the convergence properties of the multipole representation distributed multipole expansions (using $S=\sum_{b\in{\textit{\tiny B}}}S_{b}$ to denote the total number of expansion points) are used instead of one-center expansions. For the expansion over induced moments in Eqs. (70) and (IV.1), the dipole approximation is introduced and only the dipole–dipole polarizability tensor is taken into account. In this case, Eq. (IV.1) gives a set of coupled equations determining the induced dipole momentsApplequist, Carl, and Fung (1972)

[TABLE]

where the polarizability tensors for the individual sites have been replaced by a $(3S{\times}3S)$-dimensional classical linear response matrix (also known as the relay matrix) given by

[TABLE]

This matrix holds the inverse of the distributed electronic dipole–dipole polarizability tensors on the diagonal and second-order interaction tensors in the off-diagonal blocks. Upon contraction with unit vectors, $\bm{\mathcal{R}}(\omega)$ models the dipole–dipole polarizability of the environment.

By combining Eqs. (67) and (71) in truncated and distributed forms, we finally obtain the embedding operator defining the PE model:

[TABLE]

The induced dipoles, and in turn the embedding operator, depend on the wave function of subsystem A through the electric fields. In other words, upon averaging over the environment wave functions, the Hamiltonian turns into a nonlinear effective operator.

IV.2 Response theory framework

The extension of the PE model within a quantum-mechanical response framework usually proceeds by assuming the classical description of the environment from the outset, i.e., starting from Eq. (82) and the associated energy functional (see, e.g., Ref. 66). For the present analysis, we will instead begin from the response expressions derived in Sec. III.1 and outline the additional assumptions that lead to the expressions within the PE framework. As briefly alluded to in Sec. III.1, the differentiated treatment in the PE model is achieved by a perturbation analysis of the quantities in the linear response function, although, as we will presently discuss, the choice of truncation is not fully coherent from a perturbation theory point of view. Since special attention is given to subsystem A, and the sole purpose of subsystem B in this context is to obtain a realistic description of the properties of A, the order counting will be performed on the effective subsystem A quantities.

As a first step toward the PE model, one includes terms in the pure subsystem A blocks and vectors through second order (in the meaning of Eq. (IV.1)), using a first-order corrected wave function for the environment that is normalized through second order. The pure B blocks as well as the coupling blocks are evaluated only through lowest nonvanishing order. That is, only the zeroth-order contribution to the wave function of the environment subsystem is included. For the electronic Hessian, this implies that the $\mathbf{E}_{{\textit{\tiny B}}}^{[2]}$ matrix must be known through zeroth order and the $\mathbf{E}_{\textit{\tiny AB}}^{[2]}$ matrix through first order. Accordingly, the electronic Hessian and metric matrices are approximated as

[TABLE]

Although the pure B block is treated only through zeroth order, its effect on excitations in A is correct through second order as can be seen from the resulting effective electronic Hessian for subsystem A

[TABLE]

This order truncation thus provides excitation energies of A-dominated transitions that are consistent through second order. The explicit expressions for the terms are given by

[TABLE]

using Eq. (57) for a single environment subsystem and implying that $|0_{\textit{\tiny A}}\rangle$ has been derived from the effective Hamiltonian including Eq. (57) rather than the full operator in Eq. (10). The metric matrices retain their structures in Eq. (20). Based on the chosen truncation, the expression for the effective property gradient for subsystem A becomes

[TABLE]

The analogy to SOPPA is thus imperfect, since keeping only the zeroth-order correction to the B part of the property gradients means that the effective transition moments and thereby the linear response function are not consistent through second order.

To arrive at the PE model, we further need to decompose the environment into individual subsystems and invoke a truncated multipole representation of the interaction operator with respect to the environment subsystems. For practical feasibility but without theoretical justification, the lowest-order approximation invoked for the combined environment is also employed for all individual subsystems constituting the environment, meaning that the ground-state polarization among the environment subsystems, otherwise implied in $\mathbf{E}^{[2](0)}_{\textit{\tiny B}}$, is neglected. Taking the simplest two-subsystem environment as an example, the structure of the PE analog of II in Eq. (86) then takes the form illustrated in Fig. 1.

In particular, upon rewrite of the matrix resolvent according to Eq. (30), contraction with the property gradients and repeated use of Eq. (41) on the resulting subsystem blocks, we recognize the series expansions of the corresponding blocks in the relay matrix. For instance, the first block can be rewritten as

[TABLE]

that is, in terms of the relay matrix in Eq. (81). Rewriting the second term of Eq. (87) in a similar manner, we finally obtain the PE analogs of the effective electronic Hessian and effective property gradient defined in Ref. 66:

[TABLE]

taking $\hat{\mu}_{\alpha}$ as perturbation and using $\mathbf{e}_{\alpha}$ to denote a unit vector in the Cartesian $\alpha$ direction. The above expressions define the environmental effects included in the PE model within a linear response framework. Specifically, Eq. (89) defines the static (through $\hat{v}_{\text{PE}}$) and dynamic reaction field effects (second term), while the second term of Eq. (90) defines the effective external field effect.List, Jensen, and Kongsted (2016) In this way, we have shown the transition from a full quantum-mechanical description of the linear response of the combined system to the differentiated subsystem treatment in the PE model. Based on the insight from the theoretical analysis in Sec. III.1, these effective quantities can then be used in the first term of the SD (NSD) to obtain the subsystem A contributions to molecular (transition) properties.

A commonly adopted possibility to further reduce the computational complexity of the calculations is to assume frequency-independent environment subsystems, corresponding to imposing $\mathbf{S}_{{\textit{\tiny B}}}^{[2]}=\mathbf{0}$. This zero-frequency (ZF) approximation offers significant simplifications: (i) the nonlinearity of the effective electronic Hessian is lost such that Eq. (39) reduces to a standard generalized eigenvalue problem, (ii) the renormalization factor, otherwise needed in Eq. (45), becomes unity because the excitation is restricted to subsystem A in this approximation, and (iii) the additional zeroth-order poles in the first term of Eq. (50) corresponding to the FP approximation are removed. According to previous studies,Harczuk, Vahtras, and Ågren (2015); Nørby et al. (2016) the zero-frequency limit is a good approximation at off-resonant and optical frequencies because of the larger excitation energy where the frequency dispersion in the environment subsystems is typically small. In such cases, the dynamical response of the environment on the excitations in subsystem A is captured reasonably well by the static limit.

V Numerical illustration

To illustrate the basic features of the response of a combined system to a perturbing external field and the importance of the various intersubsystem interactions, we will in this section perform a numerical inspection of the working expressions presented in Sec. II. For this purpose, we consider a simplified six-level-model (SLM) for a para-nitroaniline(pNA)–water complex and its linear response to a uniform electric-field perturbation. In addition to the respective ground states, the SLM includes also the first and second singlet excited states for pNA—these are the $n\pi^{*}$ state and the intramolecular amino-to-nitro charge-transfer transition referred to as $\pi\pi^{*}$—and the first and third singlet excited states for water—these are states $1B_{1}$ and $1A_{1}$, respectively, using the symmetry labels referring to the irreducible representations of the $C_{2v}$ point group of the parent molecule. The manifold of states in the SLM is depicted in Fig. 2, defining pNA and water as subsystem A and B, respectively. The set of monomer parameters and electronic couplings reported in Tables 1 and 2 have been obtained at the TD-DFT level of theory employing CAM-B3LYPYanai, Tew, and Handy (2004)/aug-cc-pVDZDunning (1989) in the presence of the ground-state-frozen embedding potential of the other system (i.e., corresponding to the FP approximation), as obtained from a PE calculation. In other words, the second term in Eq. (89) is excluded in the response calculation. The embedding potentials consisted of atom-centered permanent electric multipoles up to quadrupoles and anisotropic electric dipole–dipole polarizabilities and were computed according to the LoPropGagliardi, Lindh, and Karlström (2004) scheme using DALTON and the Loprop-for-Dalton Python script.Vahtras (2014) The calculations were performed using a development version of the DALTON programAidas et al. (2014); dal that contains the implementation of electronic couplings between transition densities.Steinmann and Kongsted (2015) The nuclear configuration of the complex has been taken from a molecular dynamics simulation.Sneskov et al. (2011) We use this polarized basis as approximate representation for the eigenvectors of the subsystem electronic Hessians, such that the diagonal blocks of the full Hessian are diagonal. As expected from the strengths of the leading-order transition dipoles and their relative orientations (the charge-transfer transition is directed along the $x$-axis), only the electronic coupling between the $\pi\pi^{*}$ and $1A_{1}$ states is significant whereas the $n\pi^{*}$ and $1B_{1}$ states will essentially be unaffected. However, the state mixing still remains small because their energy difference is significantly larger than their electronic coupling. Therefore, to better illustrate the different aspects of the response of a combined system, the coupling blocks of the electronic Hessian have been scaled by a factor of 12. The largest absolute intersubsystem coupling element is then equal to 0.0212 a.u., which is smaller than the relevant difference between the excitation energies in the subsystems by a factor of $\sim$9.

Fig. 3a shows the isotropic electric dipole–dipole polarizability, Eq. (16), of the pNA–water complex within the SLM. The excitation energies of the full system (vertical dotted lines) can be identified as the poles of the linear response function. In the frequency region around the two lowest poles ($\omega\approx 0.10-0.20$ a.u.), the isotropic polarizability is dominated by the $\alpha_{xx}$ component of the tensor with a dispersion that is in turn dictated by the charge-transfer transition (second pole). This leads to a seemingly absence of a pole at the first excitation in subsystem A, but it is thus merely a consequence of the $n\pi^{*}$-transition being nearly electric-dipole forbidden and close in energy to the intense $\pi\pi^{*}$-transition.

The symmetric decomposition of the polarizability according to Eq. (III.1.1) into subsystem A and B contributions is shown in Fig. 3b. Note that while the diagonal elements of the polarizability of the combined system in the static limit are guaranteed to be positive, the same does not hold for the individual subsystem contributions. As exemplified in Fig. 4 by the subsystem A contribution to $\alpha_{yy}$ (the first term in Eq. (III.1.1)), this situation does indeed occur in the present case. As seen in Eq. (34), there are two contributions to the modified property gradient $\widetilde{\mathbf{V}}^{\omega}_{{\textit{\tiny A}},\beta}$. The contribution from the first term, i.e., $\mathbf{V}^{\omega}_{{\textit{\tiny A}},\beta}$, is guaranteed to give a positive contribution to $\alpha_{yy}$ for subsystem A but, due the large response in subsystem B, the second term in the effective property gradient becomes dominant and leads to an overall negative value for $\alpha_{yy}(0)$ of subsystem A.

Furthermore, we note that both the subsystem A and B contributions in the symmetric decomposition contain poles at all the excitation energies of the combined system, also referred to as the physical excitations. As we discussed previously, this implies that transition moments, in contrast to excitation energies, cannot be determined from any single one of the terms. With the nonsymmetric decomposition (Eq. (42), Fig. 3c), on the other hand, all the physical poles are collected in the first term, and it is thus the proper choice when determining residues for the transitions mainly located in subsystem A. In addition to the physical poles, however, there are unphysical zeroth-order poles in both the first and second terms of Eq. (42). Specifically, they contain poles at the excitation energies of subsystem B within the FP approximation, i.e., where the environment polarization is fixed during the response calculation. This inclusion of both physical and unphysical poles in the first term of Eq. (42) becomes particularly apparent in the frequency region close to the fourth excitation in Fig. 3c (see red solid line at $\omega\approx 0.35$ a.u.).

The excitation energies and associated one-photon transition strengths are reported in Table 3 for the two lowest transitions in the model system, i.e., those predominantly localized in subsystem A.

First, the results provide a clear evidence for the equivalence between the properties obtained from the decomposed subsystem expressions in Eqs. (39) and (45) and from the consideration of the full system expressed in terms of Eq. (25). We further consider various approximate models that are defined according to what terms are retained in the response expressions. The impact of the renormalization factor defined in Eq. (44) and appearing in Eq. (45) depends on the degree of delocalization of the given transition, and its neglect (denoted by ”$-$renorm” in Table 3) would lead to an overestimation of transition strengths. The effective field strengths experienced by subsystem A in the presence of subsystem B can be different for different directions as a consequence of the anisotropy of the polarizability of B and the relative orientation of the interacting subsystems, as seen by comparing to the assumption that only subsystem A interacts with the external field (denoted by ”$-$EEF” in Table 3). For excitations in subsystem A far from resonances in B (the environment), as in the present case, the impact of renormalization is smaller than that of the EEF effect, as expected from the relative distance dependency of the two contributions.

As discussed in Sec. IV.2, the various polarizable embedding models typically assume the effective electronic Hessian of subsystem A to be frequency-independent, thereby losing the nonlinearity of the linear response equations and the poles associated with excitations predominantly in subsystem B. The zero-frequency limit of the inverse of the isotropic polarizability of the combined system is illustrated in Fig. 5 along with the full frequency-dependent counterpart. Associated excitation energies and transition moments are given in Table 3. As anticipated from the weak dispersion of the real polarizability for subsystem B at the resonance frequencies in subsystem A (see Fig. 3c), the neglect of transitions in subsystem B only leads to a small changes in the transition properties of the A-dominated excitations. It is noted that the ZF approximation is analogous to the familiar adiabatic approximation in TD-DFT where the exchange–correlation kernel is assumed to be frequency-independent and hence leads to a linear eigenvalue problem with solutions only at Kohn–Sham one-electron excitations.Maitra et al. (2004); Cave et al. (2004); Casida and Huix-Rotllant (2015) Double- and higher-electron excitations and their effects on the Kohn–Sham single excitations are thus neglected in adiabatic TD-DFT,Caballero, Moreira, and Bofill (2013) as are the excitations in subsystem B in our case.

As discussed in Sec. III.2, the coupling of subsystem excitations gives rise to intensity borrowing. To illustrate this effect, we report in Fig. 6a the linear absorption cross section defined in Eq. (52) for the full system together with that associated with subsystem A, i.e., the first term of Eq. (50). In Fig. 6b, the corresponding cross sections for the subsystems within the FP approximation are reported. First, we note that upon integrating the absorption cross section for the full system (black solid line in Fig. 6a) across the full frequency range, we obtain a reference value of $I_{{\textit{\tiny AB}}}=0.0829$ a.u. Due to the very limited description of the electronic structure of the pNA–water complex by means of the SLM, this value is far from the exact value of 11.81 a.u. as obtained from the conservation law in Eq. (53) for a system with 82 electrons. But this discrepancy is of no concern here since the key point is that the value for the integrated cross section is identical to the corresponding summed result for the two subsystems within the FP approximation ($I_{\textit{\tiny A}}^{\text{FP}}=0.0606$ and $I_{\textit{\tiny B}}^{\text{FP}}=0.0223$ a.u., respectively).

Let us now pass to examine the nonsymmetric decomposition of the interacting system with individual cross sections $\sigma_{1}^{\mathrm{NSD}}(\omega)$ and $\sigma_{2}^{\mathrm{NSD}}(\omega)$ for the two terms in Eq. (50), respectively. These two individual cross sections obey their own respective conservation laws as expressed in Eq. (56). Since the nonsymmetric decomposition is made such that the second term corresponds exactly to subsystem B within the FP approximation, the cross section $\sigma_{2}^{\mathrm{NSD}}(\omega)$ is identical to the red solid line in Fig. 6b.

Finally, we include in Fig. 6 also the partially integrated absorption cross section for two separate finite energy intervals, shown as gray-shaded areas in the figure. In agreement with the increase in transition strength found in Table 3, the lowest band dominated by the $\pi\pi^{*}$ state is intensified upon coupling the excitations in the two subsystems. This gain in intensity in one part of the spectrum is counteracted by a reduction of the intensity in another part and, in the nonsymmetric decomposition, this will be seen only in $\sigma_{1}^{\mathrm{NSD}}(\omega)$. Since the coupling of the lowest subsystem A-dominated band (second transition) is only effective to the subsystem B-dominated $1A_{1}$-band (fourth transition), there will be an intensity borrowing from the latter to the former. As a consequence, $\sigma_{1}^{\mathrm{NSD}}(\omega)$ will take on negative values in the region of the fourth transition ($\omega\approx 0.35$ a.u.) as seen in Fig. 6a (green solid line). Accordingly, in frequency regions dominated by transitions in subsystem A, the calculation of $\sigma_{1}^{\mathrm{NSD}}(\omega)$ is to be associated with an absorption spectrum, whereas this cannot be readily done in frequency regions that includes transitions in subsystem B.

VI Summary and Conclusions

In this work, we have provided a rigorous derivation of the various contributing terms appearing in linear response theory of a quantum molecular system embedded in a polarizable environment. The origin of the three distinct types of mechanisms for the environmental effects within polarizable embedding—the static and dynamic environment responses entering the intra- and intersubsystem blocks of the electronic Hessian, respectively, as well as the effective external field effect modifying the property gradient—follow in a straightforward manner from a subsystem partitioning of a quantum-mechanical direct-product treatment of the response of the entire system. In particular, the effective external field effect is a consequence of the direct interaction between the environment and the probing external field and leads to the definition of effective subsystem properties.

A crucial point is what decomposed form of the response function that ought to be used in a given context. In the present theoretical analysis, we have demonstrated the basic features of two alternative subsystem decompositions, which clarify such discussions and highlight potential issues in defining subsystem contributions to response and transition properties.

Our first decomposition provided in Eq. (III.1.1) treats subsystems on an equal footing and results in two linear response function terms that express the permutation of subsystem indices. For that reason, we have denoted this decomposition scheme as symmetric. This is the natural choice for computing subsystem contributions to molecular properties, but, as already shown in the framework of subsystem DFT,Pavanello (2013) the coupling of the subsystems manifests itself in that each subsystem contribution contains the poles of the entire system—a feature that emphasizes the approximate nature of our picture of ”localized” transitions. As a result that is made clear in the present work, transition strengths for a subsystem-dominated excitation cannot be found from a residue analysis of the symmetric subsystem linear response function.

Our second decomposition provided in Eq. (42) treats, on the other hand, subsystems on an unequal footing and results in two linear response function terms that both are symmetric with respect to left and right property gradients. This is an unnatural choice for the development of polarizable embedding models due to the fact that the first term of the response function (describing the response properties of the chromophore of interest) contains poles not only from all transitions in the fully interacting system but also poles from the ground-state polarized, but otherwise uncoupled, environment. It is demonstrated that the linear response function in the nonsymmetric decomposition lends itself to a determination of transition strengths but only after taking proper account of renormalization of the subsystem excitation vectors. However, the renormalization factor is not accessible within the framework of polarizable embedding because it cannot be rewritten in terms of the response kernel of the environment. This complication can be avoided in practice by assuming the static limit for the environment response or by turning to the framework of complex linear response theory. We have shown that the integrated absorption cross sections of the two terms in the symmetric decomposition of the linear response function are preserved independently of each another and, as a consequence, it is demonstrated (also in the numerical example of the water–pNA dimer) that subsystem absorption cross sections (and likewise oscillator strengths) will take on negative values if intersubsystem intensity borrowing takes place.

Acknowledgements.

N. H. L. thanks J. Oddershede (University of Southern Denmark, Odense, Denmark) for helpful discussions, and the Carlsberg Foundation for a postdoctoral fellowship (Grant No. CF15-0792). P. N. acknowledges financial support from the Swedish Research Council (Grant No. 621-2014-4646). J. K. thanks the Danish Council for Independent Research (the Sapere Aude program) and the Villum Foundation for financial support. Computation/simulation for the work described in this paper has been supported by the DeIC National HPC Center, University of Southern Denmark (SDU).

Appendix A Derivation of Eq. (53)

In this appendix, we show how to derive Eq. (53). The diagonal representation of the complex linear response function is

[TABLE]

where we have used the definition of the transition strengths in Eq. (28) as well as the relation between eigenvectors for paired eigenvalues, i.e., $\mathbf{X}_{n}=\begin{bmatrix}\mathbf{X}_{n}^{1}&\mathbf{X}_{n}^{2}\end{bmatrix}$ and $\mathbf{X}_{-n}=\begin{bmatrix}\mathbf{X}_{n}^{2*}&\mathbf{X}_{n}^{1*}\end{bmatrix}$.Olsen and Jørgensen (1985) Introducing the dispersion and absorption lineshape functions

[TABLE]

the real and imaginary components of the complex linear response function can be written as

[TABLE]

The integrated absorption cross section can therefore be written as

[TABLE]

where we have used that $\sigma(\omega)$ is an even function. Each term in the square bracket can be written as a weighted sum of the dispersion and absorption lineshape functions according to

[TABLE]

Since $\mathcal{D}_{n}(\pm\omega)$ is an odd function around $\pm\omega_{n}$, its integral over the entire frequency range vanishes. Therefore, only $\mathcal{A}_{n}(\pm\omega)$, which apart from a factor is a Lorentzian function, contributes

[TABLE]

Appendix B Derivation of Eqs. (70) and (71)

In this appendix, we shall derive the two multipole-expanded representations of the induction operator given by the first and second equality in Eq. (71). The corresponding nonexpanded representations of the operator were given in Eqs. (59) and (IV.1), respectively, but to ease the discussion, we reiterate both expressions here

[TABLE]

We begin by considering the expression provided by the first equality. By introducing the multipole-expanded form of the electrostatic potential operator of subsystem $b$ given in Eq. (64), this becomes

[TABLE]

where we have used the definition of the induced multipole moments in Eq. (IV.1). Note that the multi-index summation excludes zero, since there is no induced monopole. We identify $\hat{\mathcal{V}}_{{\textit{\tiny A}}}^{(k)}(\textbf{R}_{b}){=}\int\hat{\rho}_{{\textit{\tiny A}}}(\textbf{r})T_{{\textit{r}}{b}}^{(k)}\,\text{d}\textbf{r}$ as a component of the $k$’th-order derivative of the potential operator, whose expectation value gives the $k$’th-order derivative of the electrostatic potential at $\textbf{R}_{b}$, generated by the charge density of subsystem A. Since it is common to work in terms of the fields, $\hat{F}_{{\textit{\tiny A}}}^{(k)}(\textbf{R}_{b}){=}{-}\hat{\mathcal{V}}_{{\textit{\tiny A}}}^{(k)}(\textbf{R}_{b})$ for $|l|{=}1$, and its derivatives, we shall multiply Eq. (A.101) by $1{=}({-}1)^{2}$, associating a minus with the interaction operator. Hereby, we arrive at the first of the two alternative multipole-expanded representations of the induction operator given by Eq. (70)

[TABLE]

where the field and field derivative operators are defined in Eq. (73).

Proceeding to Eq. (A.100) we start by rewriting the double integral in square brackets by introducing a Taylor expansion of the two interaction operators around a point $\mathbf{R}_{b}$ inside the charge distribution of subsystem $b$

[TABLE]

where the last equality follows from the symmetry of the interaction tensors. $P_{b}^{(k,l)}$ is the generalized electronic polarizability defined in Eq. (72). The summations exclude zero, since $P_{b}^{(k,l)}$ involves transition matrix elements, which vanish for the constant monopole operator as a consequence of the orthogonality of the unperturbed wave functions of subsystem $b$. For the same reason, $P_{b}^{(k,l)}$ contains no nuclear contribution. Substituting Eq. (B) into (A.100) yields

[TABLE]

As before, we recognize $\hat{\mathcal{V}}_{b^{\prime}}^{(l)}(\textbf{R}_{b}){=}\int T_{{\textit{r}}{b}}^{(l)}\hat{\rho}_{b^{\prime}}(\textbf{r})\,\text{d}\textbf{r}$ as a component of the $l$’th-order derivative of the potential operator, which can be translated to the corresponding field and field derivative operators upon multiplication of Eq. (A.104) by $1{=}({-}1)^{2}$. In this case, we associate a minus with each of the interaction tensors. Equation (A.104) can then be rewritten as

[TABLE]

We shall represent $\langle\hat{F}_{{b^{\prime}}}(\textbf{R}_{b})\rangle_{0_{b^{\prime}}}^{(0)}$ and $\langle\hat{F}_{{b^{\prime}}}(\textbf{R}_{b})\rangle_{0_{b^{\prime}}}^{(1)}$ in terms of the permanent and induced multipole moments. By taking minus the $l$’th-order derivative of the expectation value of Eq. (64), we obtain

[TABLE]

Analogously, the first-order correction can be recast in a multipole-expanded form as

[TABLE]

Finally, by substituting Eqs. (A.106) and (A.107) into Eq. (A.105), we arrive at the multipole-expanded representation of the induction operator given by Eq. (71)

[TABLE]