A quantum-mechanical perspective on linear response theory within polarizable embedding
Nanna Holmgaard List Patrick Norman Jacob Kongsted Hans Jørgen Aagaard Jensen




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


Author to whom correspondence should be addressed. Electronic mail: nalist@kth.se
Division of Theoretical Chemistry and Biology, School of Biotechnology, KTH Royal Institute of Technology, Roslagstullsbacken 15, SE-106 91 Stockholm, Sweden

panor@kth.se
Division of Theoretical Chemistry and Biology, School of Biotechnology, KTH Royal Institute of Technology, Roslagstullsbacken 15, SE-106 91 Stockholm, Sweden

kongsted@sdu.dk
Department of Physics, Chemistry and Pharmacy, University of Southern Denmark, Campusvej 55, 5230 Odense M, Denmark

hjj@sdu.dk
Department of Physics, Chemistry and Pharmacy, University of Southern Denmark, Campusvej 55, 5230 Odense M, Denmark







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.






33.20.-t, 33.70.Ca, 33.80.-b
                             


                              
[
    Hans Jørgen Aagaard Jensen
    December 30, 2023
==============================





§ 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.<cit.> 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.<cit.> 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 <cit.> and artificial low-lying charge-transfer transitions<cit.> 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,<cit.> 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. . 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 embedding<cit.> (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.<cit.> 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.<cit.> 
Inclusion of the explicit environmental response is particularly important for processes where substantial rearrangements may occur, e.g., upon electronic excitation.<cit.> 
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.<cit.>
In polarizable embedding, this is replaced by a self-consistent determination of the model-specific embedding parameters describing the polarization of the environment.


Both FDE<cit.> and polarizable embedding approaches<cit.> 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,<cit.> but allows in a truncated form to describe a few strongly coupled local excitations as relevant for chromophoric aggregates.<cit.> 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.<cit.> This restriction has, however, been lifted by introducing phenomenological excited-state lifetimes in the response formalism, see Refs. . 

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.<cit.> For the present aim, we will focus on the explicit expressions for the polarizable embedding (PE) model<cit.> 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,<cit.> 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,<cit.> 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,<cit.> 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.  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. . Regardless the debated interpretation of this term as  
a dispersion effect,<cit.> a nonresonant excitonic coupling between subsystems<cit.> or a resonant coupling,<cit.> 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.<cit.> 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,<cit.> 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.<cit.>
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.<cit.> 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.  



§ 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.<cit.>
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.<cit.>

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.




 §.§ 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.<cit.> 
The nonrelativistic Born–Oppenheimer Hamiltonian
for the interacting system decomposes naturally as

    ℋ̂=∑_A=1^Nℋ̂_A+∑_A> B^N𝒱̂_AB ,

where ℋ̂_A is the electronic Hamiltonian of the isolated subsystem A and 
𝒱̂_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

    𝒱̂_AB =   ∑_m∈B^M_BZ_m∑_pq∈Av_pq(𝐑_m)Ê_pq+∑_n∈A^M_AZ_n∑_rs∈Bv_rs(𝐑_n)Ê_rs
    
       +∑_pq∈A∑_rs ∈B v_pq,rsÊ_pqÊ_rs+ ∑_n∈A^M_A∑_m∈B^M_BZ_nZ_m/|𝐑_n-𝐑_m| ,

where Z_n and 𝐑_n denote the charge and the position vector, respectively, of nucleus n of the M_A nuclei in subsystem A, while Z_m, 𝐑_m and M_B are the corresponding quantities belonging to subsystem B.
Ê_pq is the usual second quantization singlet one-electron excitation operator,<cit.> 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

    v_pq(𝐑)    = -∫ϕ_p^*(𝐫)ϕ_q(𝐫)/|𝐑-𝐫|d𝐫 ,
    
    v_pq,rs   = ∬ϕ_p^*(𝐫_1) ϕ_r^*(𝐫_2) ϕ_q(𝐫_1)  ϕ_s(𝐫_2)/|𝐫_1-𝐫_2| d𝐫_1 d𝐫_2 
       = -∫ϕ_p^*(𝐫_1) ϕ_q(𝐫_1) v_rs(𝐫_1) d𝐫_1 .

The first two terms in Eq. (<ref>) 
describe the instantaneous intersubsystem Coulomb interaction between the electrons in 
subsystem A and the M_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

    𝒱̂_AB=∫ρ̂_A(𝐫)𝒱̂_B(𝐫)d𝐫=∫ρ̂_B(𝐫)𝒱̂_A(𝐫)d𝐫=ρ̂^A_𝐫 𝒱̂^B_𝐫 ,

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

    ρ̂_A(𝐫)=   
    ∑_n∈A^M_AZ_nδ(𝐫-𝐑_n)-∑_pq∈Aϕ_p^*(𝐫)ϕ_q(𝐫)Ê_pq ,
    𝒱̂_B(𝐫)=   
    ∑_m∈B^M_BZ_m/|𝐫-𝐑_m|+∑_rs∈Bv_rs(𝐫)Ê_rs .[1]

The last equality of Eq. (<ref>) 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 

    |Ψ_AΨ_B…Ψ_N⟩   = |Ψ_A⟩⊗ |Ψ_B⟩⊗…|Ψ_N⟩ .

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

    ∀ A:ℱ̂_A|0_A⟩=E_A|0_A⟩ ,


where the effective subsystem Hamiltonian is given by

    ℱ̂_A   =ℋ̂_A+∑_B≠A^N⟨ 0_B| 𝒱̂_AB| 0_B⟩
       =
     ℋ̂_A + ∑_pq∈A∑_B≠A[∑_m∈B^M_BZ_mv_pq(𝐑_m)+∑_rs∈Bv_pq,rsD_B,rs]Ê_pq .

Here, D_B,rs=⟨ 0_B | Ê_rs| 0_B⟩ 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 (<ref>) is analogous to Hartree self-consistent field theory, where the orbitals replace the subsystem wave functions. In particular, the parameterization in Eq. (<ref>) 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. 



§ 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. <ref>). This will be used in Sec. <ref>. 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

    V̂^t   = ∑_±ωV̂_α^ωF_α^ωe ^-iω t
       =∑_±ω(V̂_A,α^ω+V̂_B,α^ω)F_α^ωe ^-iω t ,

where F_α^ω are Fourier amplitudes 
associated with the one-electron perturbation operators V̂_α^ω, using Greek subscripts as (possibly composite) Cartesian labels. In order to maintain Hermiticity of V̂^t, we have F_α^ω=[F_α^-ω]^* and 
V̂_α^ω=[V̂_α^-ω]^†.
Note that the additivity of the perturbation operator in the last equality of Eq. (<ref>) 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,<cit.> 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

    {Q_AB}_T    = E_AB + ∑_ω⟨ 0_A0_B|V̂_α^ω|0_A0_B⟩ F^ω_αδ_ω
       +1/2∑_ω_1,ω_2⟨⟨V̂_α^ω_1;V̂_β^ω_2⟩⟩δ_ω_1+ω_2 F_α^ω_1F_β^ω_2
       +1/6∑_ω_1,ω_2⟨⟨V̂_α^ω_1;V̂_β^ω_2,V̂_γ^ω_2⟩⟩δ_ω_1+ω_2+ω_3 F_α^ω_1F_β^ω_2F_γ^ω_2+⋯ ,

where {·}_T indicates that the time-average over one period of oscillation in the external field has been taken. The symbol δ_ω, 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. .  




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

    |0_A 0_B⟩   = e^i[Λ̂_A(t)+Λ̂_B(t)]|0_A0_B⟩
       =e^iΛ̂_A(t)|0_A⟩⊗ e^iΛ̂_B(t)|0_B⟩; [Λ̂_A,Λ̂_B]=0 .

The intersubsystem commutation relation above follows from the zero-overlap assumption.
The time-dependent Hermitian operator Λ̂_A(t) for subsystem A is parameterized 
in terms of a set of time-dependent amplitudes {λ_i_A} and takes the form<cit.>

    Λ̂_A(t) = ∑_i> 0 (λ_i_A(t)q̂_i_A^†+λ^*_i_A(t)q̂_i_A) =𝐐_A^†Λ_A ,
    𝐐_A^†=[ 𝐪_A^†   𝐪_A ];Λ_A=[   λ_A(t) λ_A(t)^* ]^T ,

where an identity operator for subsystem B is implied.
The state-transfer operators
q̂_i_A^†=|i_A⟩⟨ 0_A| and their adjoints 
are built from a set of orthonormalized states {|i_A⟩} 
that spans the orthogonal complement space of the reference state |0_A⟩, 
the latter satisfying the variational condition given by
Eq. (<ref>). Before proceeding, we briefly comment on the employed parameterization: (i) The direct-product parameterization in Eq. (<ref>) is nonlinear: that is, it produces states outside the excitation manifold, defined by the linear action of the operator Λ̂_A(t)+Λ̂_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 Λ̂_AB(t) operator in the exponent of the first equality in Eq. (<ref>) and thus a direct coupling between direct-product states of the type ⟨ i_A0_B| and |i_Aj_B⟩ as well as ⟨ i_A0_B| and |k_Aj_B⟩ 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

    Q_AB(t)= ⟨0_A 0_B | (ℋ̂ +V̂^t-i∂/∂ t) |0_A 0_B⟩ .

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. (<ref>), 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. (<ref>), evaluated at zero field strengths.
The linear response function becomes

    ⟨⟨V̂_α^-ω; V̂_β^ω⟩⟩ =.d^2{Q_AB}_T/dF_α^-ω dF_β^ω|_𝐅=0= -𝐕_α^ω†(𝐄^[2]-ω𝐒^[2])^-1𝐕_β^ω .

Ordering the Fourier components of the configuration amplitudes in the operators in Eq. (<ref>) according to: Λ^ω=(λ^ω_A, λ^ω *_A, λ^ω_B, λ^ω_B^*)^T, leads to the following intra- and intersubsystem blocked forms of the vectors and matrices:

    -i𝐕_β^ω=.∂^2{Q_AB}_T/∂ F_β^ω ∂Λ^ω^*|_𝐅=0=   
    -i
    [ 𝐕^ω_A,β; 𝐕^ω_B,β ] ,


    𝐄^[2]-ω𝐒^[2]=   .∂^2{Q_AB}_T/∂Λ^ω^* ∂Λ^ω|_𝐅=0
    =
    [  𝐄_A^[2] 𝐄_AB^[2]; 𝐄_BA^[2]  𝐄_B^[2] ]
    -ω[ 𝐒_A^[2]       0;       0 𝐒_B^[2] ].

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 (I=A,B), take the same overall form as for the isolated subsystems

    𝐄_I^[2]   = 
    [   𝐀^I   𝐁^I; 𝐁^I^* 𝐀^I^* ]; 𝐒_I^[2] = 
    [  1  0;  0 -1 ];𝐕_I,α^ω=
    [    𝐠_α^I; -𝐠_α^I^* ] ,

where the symbols 0 and 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

    A^A_ij=    ⟨ 0_A 
    0_B |[q̂_i_A,[ℋ̂_A+𝒱̂_AB,q̂_j_A^†]] |0_A0_B⟩
    
     =    ⟨ i_A 0_B | ℋ̂_A+𝒱̂_AB| j_A 0_B⟩-δ_i_Aj_A⟨ 0_A0_B | ℋ̂_A+𝒱̂_AB | 0_A0_B⟩
    
     =    ⟨ i_A | ℋ̂_A+ ⟨ 0_B | 𝒱̂_AB| 0_B⟩ | j_A⟩-δ_i_Aj_A⟨ 0_A | ℋ̂_A+ ⟨ 0_B |𝒱̂_AB | 0_B⟩ | 0_A⟩
    
    B^A_ij=    ⟨ 0_A 0_B |[q̂_i_A,[ℋ̂_A+𝒱̂_AB,q̂_j_A]] |0_A0_B⟩=0  , 
    
    g_α,i^A=    ⟨0_A0_B| [q̂_i_A,V̂_A,α^ω] |0_A0_B⟩
     = ⟨0_A| [q̂_i_A,V̂_A,α^ω] |0_A⟩ .

In addition to the isolated subsystem term, the electronic subsystem Hessian in Eq. (<ref>) 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 𝐁^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 𝐄_AB^[2] and its conjugate transpose 𝐄_BA^[2] in Eq. (<ref>) describe the intersubsystem coupling.
Using the same ordering as in Eq. (<ref>), the structure of this block can be written as

    𝐄_AB^[2]
    =
    [   Γ   Θ; Θ^* Γ^* ] ,

with elements given by

    Γ_ij=    ⟨ 0_A 0_B |[q̂_i_A,[𝒱̂_AB,q̂_j_B^†]] |0_A0_B⟩ 
    = ⟨ i_A 0_B |𝒱̂_AB| 0_A j_B⟩  ,
    Θ_ij=    ⟨ 0_A 0_B |[q̂_i_A,[𝒱̂_AB,q̂_j_B]] |0_A0_B⟩
    = -⟨ i_A j_B |𝒱̂_AB| 0_A 0_B⟩ .

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 𝐒_AB^[2] and 𝐒_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. (<ref>) 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. (<ref>). 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:<cit.>

    𝐄^[2]𝐗   =𝐒^[2]𝐗Ω .

We recall the well-known feature of this generalized eigenvalue problem where the eigenvalues come in pairs ±ω_n and the associated eigenvectors are related,<cit.> and we shall use positive and negative indices (ω_-n=-ω_n) to denote paired solutions. Here, 𝐗 is the matrix of eigenvectors satisfying the orthonormality relation

    𝐗^†𝐒^[2]𝐗   =σ;σ_nm=sgn(n)δ_nm ,
 
and Ω is the diagonal matrix containing the associated eigenvalues ±ω_n.

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

    𝐗=
    [  𝐗^A 𝐗^AB; 𝐗^BA  𝐗^B ],

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

    T^0n_αβ=lim_ω→ω_n(ω-ω_n)⟨⟨V̂_α^-ω;V̂_β^ω⟩⟩ =
     𝐕_α^-ω_n𝐗_n^𝐗_n^†𝐕_β^ω_n ,

and used to compute the dimensionless oscillator strength

    f_0n=2ω_n/3T_αα^0n .


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 <ref>, this alternative approach to the solution of the full system provides a direct link to polarizable embedding.




  §.§.§ 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 (𝐔 and 𝐙) may be written as 

    [[ 𝐔 𝐕; 𝐖 𝐙 ]]^-1=
    [
    [       (𝐔-𝐕𝐙^-1𝐖)^-1 -(𝐔-𝐕𝐙^-1𝐖)^-1𝐕𝐙^-1; -(𝐙-𝐖𝐔^-1𝐕)^-1𝐖𝐔^-1       (𝐙-𝐖𝐔^-1𝐕)^-1 ]] ,

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

    ⟨⟨V̂_α^-ω;V̂_β^ω⟩⟩
    =   
    - (𝐕^ω†_A,α𝐍^ω_A,β
    +𝐕^ω†_B,α𝐍_B,β^ω)
    
    
    =   -𝐕_A,α^ω†(𝐄_A^[2](ω)-ω𝐒^[2]_A)^-1𝐕_A,β^ω-𝐕_B,α^ω†(𝐄_B^[2](ω)-ω𝐒_B^[2])^-1𝐕_B,β^ω .

Since the subsystems are treated on the same footing in this representation,  Eq. (<ref>) will in the following be referred to as the symmetric subsystem decomposition (SD).
In contrast to Eq. (<ref>), the dimensions in the above expression have been reduced to those of the excitation manifolds of the individual subsystems (e.g., dim_A×dim_A) instead of that of the full system (dim_A+dim_B)×(dim_A+dim_B).
As indicated by tildes, the response vectors 𝐍^ω in Eq. (<ref>) are modified quantities that satisfy the following set of effective linear response equations

    (𝐄_A^[2](ω)-ω𝐒_A^[2])𝐍_A,β^ω   = 𝐕^ω_A,β ,
    (𝐄_B^[2](ω)-ω𝐒_B^[2])𝐍_B,β^ω   = 𝐕^ω_B,β .

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., 𝐄_A^[2]), the presence of the other subsystem manifests in explicit contributions to the electronic Hessian and property gradients. In Eqs. (<ref>) and (<ref>), this has been written compactly by introducing effective Hessians and property gradients defined, here for subsystem A, as

    𝐄_A^[2](ω)   = 𝐄^[2]_A-𝐄_AB^[2](𝐄_B^[2]-ω𝐒^[2]_B)^-1𝐄_BA^[2] ,
    𝐕_A,β^ω   = 𝐕_A,β^ω-𝐄_AB^[2](𝐄_B^[2]-ω𝐒^[2]_B)^-1𝐕_B,β^ω .

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 𝐄^[2]_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.<cit.> (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. (<ref>) and (<ref>) 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 ω.<cit.> This can be seen by rewriting the intersubsystem coupling blocks according to

    Γ_ij=   ∑_pq∈A∑_rs∈Bv_pq,rs⟨ 0_A0_B |[q̂_i_A,[Ê_pqÊ_rs,q̂_j_B^†]]|0_A0_B⟩
    
    =   ∑_pq∈A∑_rs∈Bv_pq,rs⟨ 0_A|[q̂_i_A,Ê_pq]| 0_A⟩⟨ 0_B | [Ê_rs,q̂_j_B^†]|0_B⟩
    
    =   -∑_pq∈A∑_rs∈Bv_pq,rs⟨ 0_A|[q̂_i_A,Ê_pq]| 0_A⟩⟨ 0_B | [q̂_j_B^†,Ê_rs]|0_B⟩
    
    =   g_𝒱̂_𝐫,i^Ag_ρ̂_𝐫,j^B* ,

and similarly

    Θ_ij=   -∑_pq∈A∑_rs∈Bv_pq,rs⟨ 0_A|[q̂_i_A,Ê_pq]| 0_A⟩⟨ 0_B | [q̂_j_B,Ê_rs]|0_B⟩
    
    =   -g_𝒱̂_𝐫,i^Ag_ρ̂_𝐫,j^B ,

where the vectors are defined in the same way as the property gradient in Eq. (<ref>) 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



    𝐄_AB^[2]= 𝐕_𝒱̂_𝐫^A𝐕_ρ̂_𝐫^B^† =𝐕_ρ̂_𝐫^A𝐕_𝒱̂_𝐫^B^† .

Substituting this expression into the last term of Eq. (<ref>), we obtain

    𝐄_AB^[2](𝐄_B^[2]-ω𝐒_B^[2])^-1𝐄_BA^[2]
    =   𝐕_𝒱̂_𝐫^A[ 𝐕^†_ρ̂_𝐫^B(𝐄_B^[2]-ω𝐒_B^[2])^-1𝐕_ρ̂^B_𝐫^']𝐕_𝒱̂_𝐫^'^A^†
    
    =   𝐕_𝒱̂_𝐫^A C_𝐫,𝐫'^B(ω)    𝐕_𝒱̂_𝐫^'   ^A† ,

where the factor in square brackets is identified with the frequency-dependent generalized linear polarizability of the polarized ground state of subsystem B:  C_𝐫,𝐫'^B(ω)=⟨⟨ρ̂_B(𝐫);ρ̂_B(𝐫^')⟩⟩_ω, which, upon assuming the dipole approximation for subsystem B, reduces to the electric dipole–dipole polarizability tensor α^B(ω).

The property gradients for subsystem A describe the electrostatic potential generated by a transition density. The last term in Eq. (<ref>) 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. (<ref>) 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 𝐕_A,β^ω and 𝐕_A,β^ω, 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.<cit.> 

As we shall see later, the second terms of Eqs. (<ref>) and (<ref>) 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.<cit.> 
The effective external field effect is of the same origin as the cavity-field effect in polarizable continuum models leading to effective properties,<cit.> but is defined with respect to the external probing field rather than the Maxwell field within the dielectric.<cit.>
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.<cit.>


Finally, we note that neither of the subsystem contributions to the linear response function 
in Eq. (<ref>) 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,<cit.> 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 V̂_A,α^-ω to the actual perturbing field acting on subsystem A in the presence of subsystem B.






  §.§.§ 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. (<ref>) 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. (<ref>), the excitation energies of the combined system may be determined from the effective response equations in Eq. (<ref>) upon zeroing the right-hand sides. Folding the effects of subsystem B into the equation for A yields the following pseudo generalized eigenvalue equation

    𝐄_A^[2](ω_n)𝐗_n=ω_n𝐒_A^[2]𝐗_n ,

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. (<ref>) are identical to those associated with the subsystem A part of the parent system in Eq. (<ref>), 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.<cit.> In other words, this form allows us to compute the poles associated with the A (B) block of Eq. (<ref>) 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. (<ref>) 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 𝐄_A^[2] is included in Eq. (<ref>), 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. (<ref>) couples an excitation in subsystem A with those in B by the interaction of the associated transition densities (see Eq. (<ref>)). 

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,<cit.> each subsystem contribution to the response function (i.e., symmetric decomposition in Eq. (<ref>)) 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. (<ref>)). 


To facilitate the identification of residues, we will explore an alternative partitioning of the matrix resolvent in Eq. (<ref>), which can be obtained by applying the Woodbury matrix identity
(for 𝐔 and 𝐙 nonsingular square matrices)

    (𝐙-𝐖𝐔^-1𝐕)^-1𝐖𝐔^-1=𝐙^-1𝐖(𝐔-𝐕𝐙^-1𝐖)^-1 ,
       (𝐙-𝐖𝐔^-1𝐕)^-1=𝐙^-1+𝐙^-1𝐖(𝐔-𝐕𝐙^-1𝐖)^-1𝐕𝐙^-1 ,

to the third and fourth blocks of Eq. (<ref>) (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

    ⟨⟨V̂_α^-ω;V̂_β^ω⟩⟩=   -𝐕^ω†_A,α(𝐄_A^[2]-ω𝐒_A^[2] )^-1𝐕^ω_A,β	-𝐕_B,α^ω†(𝐄_B^[2]-ω𝐒_B^[2] )^-1𝐕_B,β^ω ,

in which the subsystems, in contrast to Eq. (<ref>), are treated on an unequal footing. Equation (<ref>) 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) theory<cit.> 
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. <ref>, the expressions for polarizable embedding in a linear response framework follows, in a similar spirit as SOPPA and algebraic diagrammatic construction,<cit.> from a perturbation analysis of the individual matrices appearing in Eq. (<ref>).
Note that the second term originates from the first term on the right-hand side of Eq. (<ref>) 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. (<ref>) 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. (<ref>) 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. (<ref>) 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 ± 1 according to

    𝐗^A†_n𝐒_A^ [2]𝐗^A_n+
    𝐗^BA†_n𝐒_B^[2]𝐗^BA_n=σ_n ,

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. (<ref>), 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 ω_n  
 
    Γ_n^A-1   =𝐗_n^A†𝐒_A^ [2]𝐗_n^A+
    𝐗_n^A†𝐄_AB^[2](𝐄_B^[2]-ω_n𝐒_B^[2])^-1𝐒_B^[2](𝐄_B^[2]-ω_n𝐒_B^[2])^-1𝐄_BA^[2]𝐗_n^A .

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

    T_αβ^0n=Γ_n^A𝐕_A,α^-ω_n𝐗^A_n
    𝐗^A†_n𝐕_A,β^ω_n .

What we have achieved up to this point is to have recast Eqs. (<ref>) and (<ref>) 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. <ref>, 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.<cit.>

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. (<ref>) 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. (<ref>), 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. 






 §.§ 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 (τ_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 form<cit.>

    ⟨⟨V̂_α^-ω;V̂_β^ω⟩⟩=
    -𝐕^ω†_α
    (𝐄^[2]-(ω+iγ)𝐒^[2])^-1𝐕^ω_β .

In practice, it is customary to adopt a common lifetime, and hence damping parameter γ=(2τ)^-1, for all excited states such that the so-called relaxation matrix becomes γ=γ1.

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

    ⟨⟨V̂_α^-ω;V̂_β^ω⟩⟩=   -𝐕_A,α^ω†(𝐄_A^[2](ω)-(ω+iγ)𝐒_A^[2])^-1𝐕_A,β^ω
       -𝐕_B,α^ω†(𝐄_B^[2](ω)-(ω+iγ)𝐒_B^[2])^-1𝐕_B,β^ω .

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

    𝐄_A^[2](ω)=   𝐄_A^[2]-𝐄_AB^[2](𝐄_B^[2]-(ω+iγ)𝐒_B^[2])^-1𝐄_BA^[2] ,
    𝐕^ω_A,α=   𝐕_A,α^ω-𝐄_AB^[2](𝐄_B^[2]-(ω+iγ)𝐒_B^[2])^-1𝐕_B,α^ω .

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. (<ref>)).
By further rewriting according to Eqs. (<ref>) and (<ref>), we obtain the nonsymmetric decomposition

    ⟨⟨V̂_α^-ω;V̂_β^ω⟩⟩ =   -𝐕_A,α^ω†'(𝐄_A^[2](ω)-(ω+iγ)𝐒_A^[2])^-1𝐕_A,β^ω
        -𝐕_B,α^ω†(𝐄_B^[2]-(ω+iγ)𝐒_B^[2])^-1𝐕_B,β^ω .

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. (<ref>). To be amenable to practical implementation, these may be expressed as 
a coupled set of linear equations for the real and imaginary components<cit.> (indicated by superscripts R and I, respectively). By using Eq. (<ref>) and assuming real wave functions, we obtain the following expression for the effective response equations for subsystem A


    [[ 𝐄_A^[2]-𝐕_𝒱̂^A_𝐫C_𝐫,𝐫'^B,R(ω)𝐕_𝒱̂^A_𝐫'^† -ω𝐒_A^[2]         𝐕_𝒱̂^A_𝐫C_𝐫,𝐫'^B,I(ω)𝐕_𝒱̂ ^A_𝐫'^†+γ𝐒_A^[2];      -𝐕_𝒱̂^A_𝐫 C_𝐫,𝐫'^B,I(ω)𝐕_𝒱̂^A
    _𝐫'^† -γ𝐒_A^[2]  𝐄_A^[2]-𝐕_𝒱̂^A_𝐫C_𝐫,𝐫'^B,R(ω)𝐕_𝒱̂^A_𝐫'^†-ω𝐒_A^[2];                                                    ]]
    [[ 𝐍_A,β^ω,R; 𝐍_A,β^ω,I ]]
    
    =[[ 𝐕_A,β^ω,R-𝐕_𝒱̂^A_𝐫C_𝐫,β^B,R(ω); 𝐕_A,β^ω,I-𝐕_𝒱̂^A_𝐫C_𝐫,β^B,I(ω) ]] ,


where C_𝐫,β^B(ω)=⟨⟨ρ̂^B(𝐫);V̂^ω_B,β⟩⟩_ω. 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 σ(ω).<cit.> For a sample that is isotropic with respect to the light polarization, we have

    σ(ω)=ω/3ϵ_0 c_0Im[0pt12ptα_αα(ω)] ,

where ϵ_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 <ref>)

    I=∫_0^∞σ(ω)  dω=π/2ϵ_0 c_0∑_n>0f_n0 .

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_e. 
Accordingly, for the combined system, we have

    ∑_n>0f_n0^AB=N_e^A+N_e^B ,

and likewise for the subsystems

    ∑_n>0f_n0^vac,I=∑_n>0f_n0^FP,I=N_e^I; I=A,B .

Together, Eqs. (<ref>) and (<ref>) 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. (<ref>)–(<ref>) and recalling that the second term of the complex linear response function in the nonsymmetric subsystem decomposition in Eq. (<ref>) is identical to the linear response functions of subsystem B within the FP approximation, it follows that 

    ∫σ^NSD_1(ω) dω=π/2ϵ_0 c_0N_e^A; ∫σ^NSD_2(ω) dω=π/2ϵ_0 c_0N_e^B ,


where σ^NSD_1 and σ^NSD_2 denote the contribution to the absorption cross section from the first and second term in Eq. (<ref>), respectively.
That is, if an A-dominated transition gains in intensity due to the coupling to excitations in subsystem B, then σ^NSD_1(ω) will take on negative values around poles in subsystem B. Consequently, in these regions of the spectrum, σ^NSD_1(ω) in itself cannot be associated with an absorption spectrum but it rather becomes imperative to consider the total absorption cross section σ(ω). 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 σ^NSD_1(ω) and, as we have seen, caution is called for in the interpretation of the results of such a calculation. 





§ 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.<cit.> 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, …,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)⟩,|j_b^(0)⟩} and {E_0_b^(0),E_j_b^(0)}, for b∈B, are known, where superscripts (n) specify the order with respect to the perturbation (see below).





 §.§ 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. (<ref>) for b∈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

    𝒱̂^int=   ∑_pq∈A∑_b∈B∑_m∈ b ^M_bZ_m[v_pq(𝐑_m)+∑_rs∈bv_pq,rsD_b,rs^(0)]Ê_pq_𝒱̂^es+∑_pq∈A∑_b∈B∑_rs∈ bv_pq,rsD_b,rs^(1)Ê_pq_𝒱̂^ind ,

where an element of the zeroth- and first-order electronic 
densities of subsystem b are defined as 
D_b,rs^(0)=⟨ 0_b^(0)|Ê_rs|0_b^(0)⟩ and D_b,rs^(1)=⟨ 0_b^(1)|Ê_rs|0_b^(0)⟩ +⟨ 0_b^(0)|Ê_rs|0_b^(1)⟩, respectively.
The effective interaction operator acting on subsystem A consists of contributions from the permanent and induced charge distributions, 𝒱̂^es and 𝒱̂^ind, respectively, of the environment subsystems.
In terms of the first-order reduced density and electrostatic potential operators in Eqs. (<ref>) and (<ref>), they read

    𝒱̂^es   =
    ∑_b∈B∫ρ̂_A^e(𝐫)⟨𝒱̂_b(𝐫)⟩_0_b^(0) d𝐫 ,
    𝒱̂^ind   =
    ∑_b∈B∫ρ̂_A^e(𝐫) ⟨𝒱̂^b(𝐫)⟩_0_b^(1) d𝐫 ,

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.,  ⟨𝒱̂_b(𝐫)⟩_0_b^(0)=⟨ 0_b^(0) | 𝒱̂_b(𝐫)| 0_b^(0)⟩.
The 𝒱̂^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<cit.>

    ∀ b∈B:|0_b^(1)⟩ =    -∑_j> 0 |j_b^(0)⟩∫⟨ j_b^(0) | 𝒱̂_b(𝐫)|0_b^(0)⟩/(E_j_b^(0)-E_0_b^(0))
       ×(⟨ρ̂_A(𝐫) ⟩_0_A + ∑_b'∈B\b[⟨ρ̂_b'(𝐫)⟩_0_b'^(0) +⟨ρ̂_b'(𝐫)⟩_0_b'^(1)] ) d𝐫 ,

where the expectation value involving subsystem A is over the fully polarized state within the framework of Eq. (<ref>).
As pointed out by Stone,<cit.> this expression is however inconsistent with the first-order perturbation analysis:  
through the coupling to the other environment subsystems in the third term, Eq. (<ref>) 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. (<ref>).

Substituting the first-order correction to the wave function in Eq. (<ref>) into 𝒱̂^ind
and using the definition of the generalized static polarizability, C_𝐫,𝐫'^b,(0)(ω=0),
of the 
zeroth-order ground states of the environment subsystems, the first-order induction operator can be written as 

    𝒱̂^ind =    - ∑_b∈B∫ρ̂_A^e(𝐫) ∫[∬C_𝐫”,𝐫”'^b,(0)(0)/|𝐫-𝐫”||𝐫'-𝐫”'| d𝐫” d𝐫”' ]
       ×(⟨ρ̂_A(𝐫') ⟩_0_A + ∑_b'∈B\ b[⟨ρ̂_b'(𝐫')⟩_0_b'^(0) +⟨ρ̂_b'(𝐫')⟩_0_b'^(1)] )d𝐫 d𝐫' .


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.<cit.> The norm and factorial 
of a multi-index are defined as |k|=k_x+k_y+k_z and k!=k_x!· k_y!· k_z!, respectively, the sum of two multi-indexes as k± l=(k_x± l_x,k_y± l_y,k_z± l_z), and the multi-index power of a vector as
𝐫^k=x^k_x· y^k_y· z^k_z.
Using this notation, a Taylor series expansion of the electrostatic potential operator can be written as<cit.>

    1/|𝐑_j-𝐑_i|=
    ∑_|k|=0^∞(-1)^|k|/k!(∇_j^k 1/|𝐑_j-𝐑_o|)(𝐑_i-𝐑_o)^k ,

where 𝐑_o is the expansion point and the summation over k runs over the 12(|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 operator<cit.> 

    T_ij^(k) = ∂^k_j1/|𝐫_j-𝐫_i| ;    ∂^k_j = ∂^|k|/∂ x_j^k_x∂ y_j^k_y∂ z_j^k_z ,

where the superscript multi-index notation should not be confused with the perturbation order.
Note that in writing Eq. (<ref>), 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

    𝒱̂_b(𝐑_b')=∑_|k|=0^∞(-1)^|k|/k!T_bb'^(k)M̂_b^(k)(𝐑_b) ,

where the expansion point 𝐑_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

    M̂_b^(k)(𝐑_b)   = ∫ρ̂_b(𝐫)(𝐫-𝐑_b)^k  d𝐫
       =∑_m∈b^M_bZ_m(𝐑_m-𝐑_b)^k+∑_rs∈bm_rs^(k)(𝐑_b)Ê_rs ,

where the associated electronic multipole integral is given by

    m_rs^(k)(𝐑_b)=-∫ϕ_r^*(𝐫)(𝐫-𝐑_b)^kϕ_s(𝐫) d𝐫 .

The expectation value of Eq. (<ref>) corresponds to an electric multipole moment, M^(k)_b(𝐑_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 𝐑_b). For the sake of brevity, 
we will suppress this explicit origin dependence in the following.
Finally, substitution of Eq. (<ref>)
into 𝒱̂^es yields the multipole-expanded
form of the electrostatic interaction operator

    𝒱̂^es=∑_pq∈A∑_b∈B∑_|k|=0^∞(-1)^|k|/k!M_b^(k)t_pq^(k)(𝐑_b)Ê_pq ,

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

    t_pq^(k)(𝐑_b)=-∫ϕ_p^*(𝐫_i)T_bi^(k)ϕ_q(𝐫_i) d𝐫_i  ,


where the index i refers to an electronic coordinate.
Comparing Eq. (<ref>) to the nonexpanded form in Eq. (<ref>),
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. (<ref>) 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

    M̅_b^(k)(𝐑_b)   =∫⟨ρ̂_b(𝐫)⟩_0_b^(1) (𝐫-𝐑_b)^k d𝐫
       = ⟨ 0_b^(1)| M̂_b^(k)|0_b^(0)⟩+ ⟨ 0_b^(0)| M̂_b^(k)|0_b^(1)⟩ ,

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 <ref>,
the induction part of the interaction operator can be written in the multipole-expanded
form as

    𝒱̂^ind =   
    -∑_b∈B∑_|k|=1^∞1/k!F̂_A^e,(k)(𝐑_b)M̅_b^(k)
    
    =   -∑_b∈B∑_|k|=1^∞∑_|l|=1^∞1/k!· l!F̂_A^e,(k)(𝐑_b)P_b^(k,l)
       ×(⟨F̂_A^(l)(𝐑_b)⟩_0_A+∑_b'∈B\ b[⟨F̂_b'^(l)(𝐑_b)⟩^(0)_0_b'+⟨F̂_b'^(l)(𝐑_b)⟩_0_b'^(1)] ) ,

where the two alternative expressions arise by 
expanding Eqs. (<ref>) and (<ref>), 
respectively. 
P_b^(k,l) are static electronic polarizabilities of subsystem b, analogous to that in Eq. (<ref>), defined as

    P_b^(k,l)=   ∑_j_b≠ 0_b[⟨ 0_b^(0) | M̂_b^(k)| j_b^(0)⟩⟨ j_b^(0) | M̂_b^(l)|0_b^(0)⟩/(E_j_b^(0)-E_0_b^(0))+⟨ 0_b^(0) | M̂_b^(l)| j_b^(0)⟩⟨ j_b^(0) | M̂_b^(k)|0_b^(0)⟩/(E_j_b^(0)-E_0_b^(0))] ,

recalling that this definition employs the traced multipole moment operators.
Furthermore, F̂_A^(k)(𝐑_b) is the (k-1)'th-order electric-field derivative,
which probes the field derivative produced by subsystem A at point 𝐑_b

    F̂_A^(k)(𝐑_b)=-∑_n∈A^M_AZ_nT_nb^(k)_F_A^n,(k)(𝐑_b)+(-1)^|k|+1∑_pq∈A t_pq^(k)(𝐑_b)Ê_pq_F̂_A^e,(k)(𝐑_b) ,

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. (<ref>). Accordingly, we can write 
the zeroth- and first-order 
electric-field derivatives of the environment subsystems in Eq. (<ref>), 
in terms of static and induced multipole
moments, respectively, as

    ⟨F̂_b'^(l)(𝐑_b)⟩_0_b'^(0)   =∑_|k|=0^∞(-1)^|k|+1/k!T^(k+l)_b'bM_b'^(k) ,
    
    ⟨F̂_b'^(l)(𝐑_b)⟩_0_b'^(1)   =∑_|k|=1^∞(-1)^|k|+1/k!T^(k+l)_b'bM̅_b'^(k) .

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

    ∀b∈B: 
    M̅_b^(k) =   ∑_|l|=1^∞1/l!P_b^(k,l)F^tot(l)(𝐑_b)
    
    =   ∑_|l|=1^∞1/l!P_b^(k,l)(⟨F̂_A^(l)(𝐑_b)⟩_0_A
    +∑_b'∈B\b[⟨F̂_b'^(l)(𝐑_b)⟩_0_b'^(0)+∑_|m|=1^∞(-1)^|m|+1/m!T_b'b^(m+l)M̅_b'^(m)]) ,

where F^tot(l)(𝐑_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)(𝐑_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. (<ref>) and (<ref>) are terminated at a finite order K_s, and to improve the convergence properties of the multipole representation distributed multipole expansions (using S=∑_b∈BS_b to denote the total number of expansion points) are used instead of one-center expansions. For the expansion over induced moments in Eqs. (<ref>) and (<ref>), the dipole approximation is introduced and only the dipole–dipole polarizability tensor is taken into account. In this case, Eq. (<ref>) gives a set of coupled equations determining the induced dipole moments<cit.>

    μ̅_s(0)= ∑_s'=1^Sℛ_ss'(0)𝐅(𝐑_s') ,

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

    ℛ(ω) =
     ([     α_1(ω)^-1     -𝐓^(2)_12             ⋯     -𝐓^(2)_1S;     -𝐓^(2)_21     α_2(ω)^-1             ⋱             ⋮;             ⋮             ⋱             ⋱ -𝐓^(2)_(S-1)S;     -𝐓^(2)_S1             ⋯ -𝐓^(2)_S(S-1)     α_S(ω)^-1;               ])^-1.
  
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, ℛ(ω) models the dipole–dipole polarizability of the environment.

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

    v̂_PE=   ∑_pq∈A∑_s=1^S∑_|k|=0^K_s(-1)^|k|/k!M_s^(k)t_pq^(k)(𝐑_s)Ê_pq-∑_s=1^Sμ̅_s,α(0)F̂^e_A,α(𝐑_s) .
 
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.



 §.§ 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. (<ref>) and the associated energy functional (see, e.g., Ref. ). For the present analysis, we will instead begin from the response expressions derived in Sec. <ref> and outline the additional assumptions that lead to the expressions within the PE framework. 
As briefly alluded to in Sec. <ref>, 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. (<ref>)), 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 𝐄_B^[2] matrix must be known through zeroth order and the 𝐄_AB^[2] matrix through first order. Accordingly, the electronic Hessian and metric matrices are approximated as

    𝐄^[2]=
    [ 𝐄_A^[2](0,1,2)    𝐄_AB^[2](1);    𝐄_AB^[2](1)     𝐄_B^[2](0) ]; 𝐒^[2]=
    [ 𝐒_A^[2](0,1,2)              0;              0     𝐒_B^[2](0) ].

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 

    𝐄_A^[2] ^(II)=   𝐄_A^[2](0,1,2)_I-𝐄_AB^[2](1)(𝐄_B^[2](0)-ω𝐒_B^[2](0))^-1𝐄_BA^[2](1)_II .

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

    I   =⟨ 0_A | [𝐐_A,[ℋ̂_A+𝒱̂^int,𝐐_A^†]]|0_A⟩ ,
    II   = 𝐕_𝒱̂^A_𝐫C_𝐫,𝐫'^B,(0)(ω)𝐕_𝒱̂^A_𝐫'^† ,

using Eq. (<ref>) for a single environment subsystem and implying that |0_A⟩ has been derived from the effective Hamiltonian including  Eq. (<ref>) rather than the full operator in Eq. (<ref>). The metric matrices retain their structures in Eq. (<ref>).
Based on the chosen truncation, the expression for the effective property gradient for subsystem A becomes

    𝐕_A,α^ω=   𝐕_A,α^ω(2)-𝐄_AB^[2](1)(𝐄_B^[2](0)-ω𝐒_B^[2](0))^-1𝐕_B,α^ω(0) .

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 𝐄^[2](0)_B, is neglected. Taking the simplest two-subsystem environment as an example, the structure of the PE analog of II in Eq. (<ref>) then takes the form illustrated in Fig. <ref>. 

In particular, upon rewrite of the matrix resolvent according to Eq. (<ref>), contraction with the property gradients and repeated use of Eq. (<ref>) 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

    𝐕_μ̂^b_1^(0)†   ((𝐄_b_1^[2](0)-ω𝐒_b_1^[2](0))^-1-𝐄_b_1b_2^[2](1)(𝐄_b_2^[2](0)-ω𝐒_b_2^[2](0))^-1𝐄_b_2b_1^[2](0))^-1𝐕_μ̂	^b_1^(0)
       =(α_b_1(ω)^-1+𝐓^(2)_b_1b_2α_b_2(ω)𝐓^(2)_b_2b_1)^-1 ,

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

    𝐄_A^[2]   PE=⟨ 0_A| [𝐐_A,[ℋ̂_A+v̂_PE,𝐐_A^†]]| 0_A⟩
       -∑_s,s'=1^S⟨ 0_A| [𝐐_A,𝐅̂^e_A(𝐑_s)]| 0_A⟩ℛ_ss'(ω)⟨ 0_A| [𝐐^†,𝐅̂^e_A(𝐑_s')]| 0_A⟩ ,
    𝐕_A,α^ω   PE=⟨ 0_A| [𝐐_A,V̂^ω_A,α]| 0_A⟩ -∑_s,s'=1^S⟨ 0_A| [𝐐_A,𝐅̂^e_A(𝐑_s)]| 0_A⟩ℛ_ss'(ω)𝐞_α ,

taking μ̂_α as perturbation and using 𝐞_α to denote a unit vector in the Cartesian α direction. The above expressions define the environmental effects included in the PE model within a linear response framework. Specifically, Eq. (<ref>) defines the static (through v̂_PE) and dynamic reaction field effects (second term), while the second term of  Eq. (<ref>) defines the effective external field effect.<cit.>
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. <ref>, 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 
𝐒_B^[2]=0. 
This zero-frequency (ZF) approximation offers significant simplifications: (i) the nonlinearity of the effective electronic Hessian is lost such that Eq. (<ref>) reduces to a standard generalized eigenvalue problem, (ii) the renormalization factor, otherwise needed in Eq. (<ref>), 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. (<ref>) corresponding to the FP approximation are removed.
According to previous studies,<cit.> 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.




§ 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. <ref>. 
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π^* state and the intramolecular amino-to-nitro charge-transfer transition referred to as ππ^*—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. <ref>, defining pNA and water as subsystem A and B, respectively. The set of monomer parameters and electronic couplings reported in Tables <ref> and <ref> have been obtained at the TD-DFT level of theory employing CAM-B3LYP<cit.>/aug-cc-pVDZ<cit.> 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. (<ref>) 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 LoProp<cit.> scheme using DALTON and the Loprop-for-Dalton Python script.<cit.>
The calculations were performed using a development version of the DALTON program<cit.> that contains the implementation of electronic couplings between transition densities.<cit.> The nuclear configuration of the complex has been taken from a molecular dynamics simulation.<cit.> 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 ππ^* and 1A_1 states is significant whereas the nπ^* 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 ∼9. 


 

 

Fig. <ref>a shows the isotropic electric dipole–dipole           polarizability, Eq. (<ref>), 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 (ω≈ 0.10-0.20 a.u.), the isotropic polarizability is dominated by the α_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π^*-transition being nearly electric-dipole forbidden and close in energy to the intense ππ^*-transition. 
 
The symmetric decomposition of the polarizability according to Eq. (<ref>) into subsystem A and B contributions is shown in Fig. <ref>b. 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. <ref> by the subsystem A contribution to α_yy (the first term in Eq. (<ref>)), this situation does indeed occur in the present case. As seen in Eq. (<ref>), there are two contributions to the modified property gradient 𝐕^ω_A,β. The contribution from the first term, i.e., 𝐕^ω_A,β, is guaranteed to give a positive contribution to α_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 α_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. (<ref>), Fig. <ref>c), 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. (<ref>). 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. (<ref>) becomes particularly apparent in the frequency region close to the fourth excitation in Fig. <ref>c (see red solid line at ω≈ 0.35 a.u.).






The excitation energies and associated one-photon transition strengths are reported in Table <ref> 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. (<ref>) and (<ref>) and from the consideration of the full system expressed in terms of Eq. (<ref>). 
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. (<ref>) and appearing in Eq. (<ref>) depends on the degree of delocalization of the given transition, and its neglect (denoted by "-renorm" in Table <ref>) 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 <ref>). 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. <ref>, 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. <ref> along with the full frequency-dependent counterpart. Associated excitation energies and transition moments are given in Table <ref>. As anticipated from the weak dispersion of the real polarizability for subsystem B at the resonance frequencies in subsystem A (see Fig. <ref>c), 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.<cit.> Double- and higher-electron excitations and their effects on the Kohn–Sham single excitations are thus neglected in adiabatic TD-DFT,<cit.> as are the excitations in subsystem B in our case.






As discussed in Sec. <ref>, the coupling of subsystem excitations gives rise to intensity borrowing. To illustrate this effect, we report in Fig. <ref>a the linear absorption cross section defined in Eq. (<ref>) for the full system together with that associated with subsystem A, i.e., the first term of Eq. (<ref>). In Fig. <ref>b, 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. <ref>a) across the full frequency range, we obtain a reference value of I_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. (<ref>) 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_A^FP = 0.0606 and I_B^FP = 0.0223 a.u., respectively).

Let us now pass to examine the nonsymmetric decomposition of the interacting system with individual cross sections σ_1^NSD(ω) and σ_2^NSD(ω) for the two terms in Eq. (<ref>), respectively. These two individual cross sections obey their own respective conservation laws as expressed in Eq. (<ref>). Since the nonsymmetric decomposition is made such that the second term corresponds exactly to subsystem B within the FP approximation, the cross section σ_2^NSD(ω) is identical to the red solid line in Fig. <ref>b.

Finally, we include in Fig. <ref> 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 <ref>, the lowest band dominated by the ππ ^* 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 σ_1^NSD(ω). 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, σ_1^NSD(ω) will take on negative values in the region of the fourth transition (ω≈ 0.35 a.u.) as seen in Fig. <ref>a (green solid line). Accordingly, in frequency regions dominated by transitions in subsystem A, the calculation of σ_1^NSD(ω) is to be associated with an absorption spectrum, whereas this cannot be readily done in frequency regions that includes transitions in subsystem B.



§ 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. (<ref>) 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,<cit.> 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. (<ref>) 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.





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



     
  

§ DERIVATION OF EQ. (<REF>)

  In this appendix, we show how to derive Eq. (<ref>). 
  The diagonal representation of the complex linear response function is
  
    - ⟨⟨μ̂_α^-ω;μ̂_β^ω⟩⟩=∑_n>0[T_αβ^0n/ω_n-(ω+iγ)+T_βα^0n/ω_n+(ω+iγ)].

  where we have used the definition of the transition strengths in Eq. (<ref>) as well as the relation between eigenvectors for paired eigenvalues, i.e., 𝐗_n=[ 𝐗_n^1 𝐗_n^2 ] and 𝐗_-n=[ 𝐗_n^2* 𝐗_n^1* ].<cit.> Introducing the dispersion and absorption lineshape functions
  
    𝒟_n(±ω)=   ω_n∓ω/(ω_n∓ω)^2+γ^2 ,
    𝒜_n(±ω)=   γ/(ω_n∓ω)^2+γ^2 ,

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

    -Re[⟨⟨μ̂_α^-ω;μ̂_β^ω⟩⟩]=   ∑_n>0[T_αβ^0n𝒟_n(ω)+T_βα^0n𝒟_n(-ω)],
    
    -Im[⟨⟨μ̂_α^-ω;μ̂_β^ω⟩⟩]=   ∑_n>0[T_αβ^0n𝒜_n(ω)-T_βα^0n𝒜_n(-ω)],


The integrated absorption cross section can therefore be written as

    I=   ∫_0^∞σ(ω) dω = 1/2∫_-∞^∞σ(ω) dω 
    
    =   1/2ϵ_0c_0∑_n>0T_αα^0n/3∫_-∞^∞ω[γ/(ω_n-ω)^2+γ^2-γ/(ω_n+ω)^2+γ^2]
     dω ,

where we have used that σ(ω) 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

    γω/(ω_n∓ω)^2+γ^2=   ∓γ𝒟_n(±ω)±ω_n𝒜_n(±ω) .

Since 𝒟_n(±ω) is an odd function around ±ω_n, its integral over the entire frequency range vanishes. Therefore, only 𝒜_n(±ω), which apart from a factor is a Lorentzian function, contributes

    I=   1/2ϵ_0c_0∑_n>0ω_n/3T_αα^0n∫_-∞^∞[𝒜_n(ω)+𝒜_n(-ω)] dω 
    
    =   π/2ϵ_0c_0∑_n>02ω_n/3T_αα^0n
    
    =   π/2ϵ_0c_0∑_n>0f_0n .

  
  
  
 
  
  
   


§ DERIVATION OF EQS. (<REF>) AND (<REF>)
  

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

    𝒱̂^ind =   ∑_b∈B∫ρ̂^e_A(𝐫)⟨𝒱̂_b(𝐫)⟩_0_b^(1) d𝐫
    
    =    - ∑_b∈B∫ρ̂_A^e(𝐫) ∫[∬C_𝐫”,𝐫”'^b,(0)(0)/|𝐫-𝐫”||𝐫'-𝐫”'| d𝐫” d𝐫”' ]
       ×(⟨ρ̂_A(𝐫') ⟩_0_A + ∑_b'∈B\ b[⟨ρ̂_b'(𝐫')⟩_0_b'^(0) +⟨ρ̂_b'(𝐫')⟩_0_b'^(1)] )d𝐫 d𝐫' .

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. (<ref>), this becomes  

    𝒱̂^ind   = ∑_b∈B∫ρ̂^e_A(𝐫)⟨𝒱̂_b(𝐫)⟩_0_b^(1) d𝐫
       =∑_b∈B∑_|k|=1^∞(-1)^|k|/k!⟨M̂_b^(k)(𝐑_b)⟩_0_b^(1)∫ρ̂^e_A(𝐫)T^(k)_br d𝐫
       =∑_b∈B∑_|k|=1^∞1/k!M̅_b^(k)∫ρ̂^e_A(𝐫)T^(k)_rb d𝐫 ,

where we have used the definition of the induced multipole moments in Eq. (<ref>). Note that the multi-index summation excludes zero, since there is no induced monopole.
We identify 𝒱̂_A^(k)(R_b)=∫ρ̂_A(r)T_rb^(k) dr 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 R_b, generated by the charge density of subsystem A. Since it is common to work in terms of the fields, F̂_A^(k)(R_b)=-𝒱̂_A^(k)(R_b) for |l|=1, and its derivatives, we shall multiply Eq. (<ref>) 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. (<ref>)

    𝒱̂^ind   =-∑_b∈B∑_|k|=1^∞1/k!F̂_A^e,(k)(𝐑_b)M̅_b^(k) ,

where the field and field derivative operators are defined in Eq. (<ref>).

Proceeding to Eq. (<ref>) we start by rewriting 
the double integral in square brackets by introducing a Taylor expansion of the two interaction operators around a point 𝐑_b inside the charge distribution of subsystem b

    ∬C_𝐫”,𝐫”'^b,(0)(0)/|r-r”||r'-r”'| dr” dr”' =   ∑_|k|=1^∞∑_|l|=1^∞(-1)^|k|+|l|/k!· l!T_br^(k)P_b^(k,l)T_br'^(l)
    
    =   ∑_|k|=1^∞∑_|l|=1^∞1/k!· l!T_rb^(k)P_b^(k,l)T_r'b^(l) ,

where the last equality follows from the symmetry of the interaction tensors. P_b^(k,l) is the generalized electronic polarizability defined in Eq. (<ref>). 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. (<ref>) into (<ref>) yields

    𝒱̂^ind =   -∑_b∈B∑_|k|=1^∞∑_|l|=1^∞1/k!· l!∫ρ̂^e_A(r) T_rb^(k) dr P_b^(k,l)
       ×∫ T_r'b^(l)(⟨ρ̂_A(r')⟩_0_A + ∑_b'∈B\ b[⟨ρ̂_b'(r')⟩_0_b'^(0)+⟨ρ̂_b'(r')⟩_0_b'^(1)]) dr' .

As before, we recognize 𝒱̂_b'^(l)(R_b)=∫ T_rb^(l)ρ̂_b'(r) dr 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. (<ref>) by 1=(-1)^2. In this case, we associate a minus with each of the interaction tensors.
Equation (<ref>) can then be rewritten as

    𝒱̂^ind =   -∑_b∈B∑_|k|=1^∞∑_|l|=1^∞1/k!l!F̂_A^e,(k)(R_b)  P_b^(k,l)
       ×(⟨F̂^(l)_A(R_b)⟩_0_A + ∑_b'∈B\ b'
    [⟨F̂_b'^(l)(R_b)⟩_0_b'^(0)+⟨F̂^(l)_b'(R_b)⟩_0_b'^(1)]) .

We shall represent ⟨F̂_b'(R_b)⟩_0_b'^(0) and ⟨F̂_b'(R_b)⟩_0_b'^(1) in terms of the permanent and induced multipole moments. By taking minus the l'th-order derivative of the expectation value of Eq. (<ref>), we obtain 

    ⟨F̂_b'^(l)(R_b)⟩_0_b'^(0)=∑_|k|=0^∞(-1)^|k|+1/k!T_b'b^(k+l)M_b'^(k) .
  
Analogously, the first-order correction can be recast in a multipole-expanded form as

    ⟨F̂^(l)_b'(R_b)⟩_0_b'^(1)=    -∫ T_rb^(l)⟨ρ̂_b'(r)⟩_0_b'^(1) dr
    
    =   -∑_|k|=1^∞(-1)^|k|/k!T_b'b^(k+l)∫(r-R_b)^k⟨ρ̂_b'(r)⟩_0_b'^(1) dr
    
    =   ∑_|k|=1^∞(-1)^|k|+1/k!T_b'b^(k+l)M̅_b'^(k) .

Finally, by substituting Eqs. (<ref>) and (<ref>) into Eq. (<ref>), we arrive at the multipole-expanded representation of the induction operator given by Eq. (<ref>)

    𝒱̂^ind =   -∑_b∈B∑_|k|=1^∞∑_|l|=1^∞1/k!· l!F̂_A^e,(k)(R_b)  P_b^(k,l)
       ×(⟨F̂^(l)_A(R_b)⟩_0_A + ∑_b'∈B\b[⟨F̂_b'^(l)(R_b)⟩_0_b'^(0)+∑_|m|=1^∞(-1)^|m|+1/m!T_b'b^(m+l)M̅_b'^(m)]) .














92
fxundefined [1]
 ifx#1

fnum [1]
 #1firstoftwo
 secondoftwo
 
fx [1]
 #1firstoftwo
 secondoftwo
 





noop [0]secondoftwo

ref[1]@startlink#1@href
href[1]#1@endlink
anitize@url [0]`
12`$12`&12`#12`1̂2`_12`%12
startlink[1]
endlink[0]

rl [1]href #1












@bib@innerbibempty

[Lindon, Tranter, and Holmes(2010)]Lindon2010
  editor J. C. Lindon, editor G. E. Tranter,  and editor J. L. Holmes, eds., @noop
   title Encyclopedia of spectroscopy and
  spectrometry (publisher Academic Press, address San Diego, year 2010)NoStop
[Kristensen et al.(2012)Kristensen, Høyvik, Jansik,
  Jørgensen, Kjærgaard, Reine, and Jakowski]kristensen2012mp2
  author author K. Kristensen, author I.-M. Høyvik, author B. Jansik,
  author P. Jørgensen, author T. Kjærgaard, author
  S. Reine,  and author
  J. Jakowski, title title MP2 Energy and Density for Large Molecular Systems with
  Internal Error Control Using the Divide-Expand-Consolidate Scheme,  @noop  journal journal Phys. Chem. Chem.
  Phys. volume 14, pages 15706–15714
  (year 2012)NoStop
[Ziółkowski et al.(2010)Ziółkowski, Jansik, Kjærgaard,  and Jørgensen]ziolkowski2010linear
  author author M. Ziółkowski, author B. Jansik, author T. Kjærgaard,  and author P. Jørgensen, title title Linear
  Scaling Coupled Cluster Method with Correlation Energy Based Error
  Control, @noop  journal journal J.
  Chem. Phys. volume 133, pages 014107
  (year 2010)NoStop
[Ochsenfeld, Kussmann, and  Lambrecht(2007)]ochsenfeld2007linear
  author author C. Ochsenfeld, author J. Kussmann,  and author D. S. Lambrecht, title title
  Linear-Scaling Methods in Quantum Chemistry, @noop  journal journal Rev. Comp. Ch. volume
  23, pages 1 (year 2007)NoStop
[Sherrill(2010)]sherrill2010frontiers
  author author C. D. Sherrill, title title Frontiers in
  Electronic Structure Theory, @noop  journal
  journal J. Chem. Phys. volume 132,  pages 110902 (year 2010)NoStop
[Jakobsen, Kristensen, and  Jensen(2013)]jakobsen2013electrostatic
  author author S. Jakobsen, author K. Kristensen,  and author F. Jensen, title title Electrostatic
  Potential of Insulin: Exploring the Limitations of Density Functional Theory
  and Force Field Methods, @noop  journal journal J. Chem. Theory Comput. volume 9,  pages 3978–3985 (year 2013)NoStop
[Olsen et al.(2015a)Olsen, List,
  Kristensen, and Kongsted]olsen2015accuracy
  author author J. M. H.  Olsen, author N. H. List, author K. Kristensen,   and author J. Kongsted,  title title Accuracy of Protein
  Embedding Potentials: An Analysis in Terms of Electrostatic Potentials,  @noop  journal journal J. Chem. Theory
  Comput. volume 11, pages 1832–1842
  (year 2015a)NoStop
[Dreuw, Weisman, and Head-Gordon(2003)]dreuw2003long
  author author A. Dreuw, author J. L. Weisman,
   and author M. Head-Gordon,  title title Long-Range Charge-Transfer
  Excited States in Time-Dependent Density Functional Theory Require Non-Local
  Exchange, @noop  journal journal J.
  Chem. Phys. volume 119, pages
  2943–2946 (year 2003)NoStop
[Gritsenko and Baerends(2004)]gritsenko2004asymptotic
  author author O. Gritsenko and author E. J. Baerends, title title Asymptotic
  Correction of the Exchange–Correlation Kernel of Time-Dependent Density
  Functional Theory for Long-Range Charge-Transfer Excitations, @noop
   journal journal J. Chem. Phys. volume 121, pages 655–660 (year
  2004)NoStop
[Dreuw and Head-Gordon(2004)]dreuw2004failure
  author author A. Dreuw and author M. Head-Gordon, title title Failure of
  Time-Dependent Density Functional Theory for Long-Range Charge-Transfer
  Excited States: The Zincbacteriochlorin-Bacteriochlorin and
  Bacteriochlorophyll-Spheroidene Complexes, @noop  journal journal J. Am. Chem. Soc. volume 126, pages 4007–4016 (year
  2004)NoStop
[Gordon et al.(2012)Gordon,
  Fedorov, Pruitt, and Slipchenko]gordon2012fragmentation
  author author M. S. Gordon, author D. G. Fedorov,
  author S. R. Pruitt,  and  author L. V. Slipchenko,  title title Fragmentation Methods: a
  Route to Accurate Calculations on Large Systems, @noop 
  journal journal Chem. Rev volume 112, pages 632–672 (year
  2012)NoStop
[Jacob and Neugebauer(2014)]jacob2014subsystem
  author author C. R. Jacob and author J. Neugebauer, title title Subsystem
  Density-Functional Theory, @noop  journal journal WIREs Comput. Mol. Sci. volume 4,  pages 325–362 (year 2014)NoStop
[Gomes and Jacob(2012)]gomes2012quantum
  author author A. S. P.  Gomes and author C. R.  Jacob, title title
  Quantum-Chemical Embedding Methods for Treating Local Electronic
  Excitations in Complex Chemical Systems, @noop  journal journal Annu. Rep. Prog. Chem., Sect. C: Phys.
  Chem. volume 108, pages 222–277
  (year 2012)NoStop
[Wesolowski and Warshel(1993)]wesolowski1993frozen
  author author T. A. Wesolowski and author A. Warshel, title title Frozen Density
  Functional Approach for Ab Initio Calculations of Solvated Molecules,  @noop  journal journal J. Phys. Chem.  volume 97, pages 8050–8053 (year 1993)NoStop
[Warshel and Levitt(1976)]warshel1976theoretical
  author author A. Warshel and author M. Levitt, title title Theoretical
  Studies of Enzymatic Reactions: Dielectric, Electrostatic and Steric
  Stabilization of the Carbonium Ion in the Reaction of Lysozyme,  @noop  journal journal J. Mol. Biol.  volume 103, pages 227–249 (year 1976)NoStop
[Senn and Thiel(2009)]senn2009qm
  author author H. M. Senn and author W. Thiel,  title title QM/MM Methods for
  Biomolecular Systems, @noop  journal journal Angew. Chem. Int. Ed. volume 48,  pages 1198–1229 (year 2009)NoStop
[van der Kamp and Mulholland(2013)]van2013combined
  author author M. W. van der Kamp and author A. J. Mulholland, title title Combined
  Quantum Mechanics/Molecular Mechanics (QM/MM) Methods in Computational
  Enzymology, @noop  journal journal
  Biochemistry volume 52, pages
  2708–2728 (year 2013)NoStop
[Lopes, Roux, and MacKerell Jr(2009)]lopes2009molecular
  author author P. E. M.  Lopes, author B. Roux,  and author A. D. MacKerell Jr, title title Molecular
  Modeling and Dynamics Studies with Explicit Inclusion of Electronic
  Polarizability: Theory and Applications, @noop  journal journal Theor. Chem. Acc. volume 124, pages 11–28 (year
  2009)NoStop
[Sneskov et al.(2011)Sneskov, Schwabe, Christiansen, and  Kongsted]sneskov2011scrutinizing
  author author K. Sneskov, author T. Schwabe,
  author O. Christiansen,  and  author J. Kongsted, title title Scrutinizing the Effects of
  Polarization in QM/MM Excited State Calculations, @noop 
  journal journal Phys. Chem. Chem. Phys.  volume 13, pages 18551–18560 (year 2011)NoStop
[Wesolowski and Weber(1996)]wesolowski1996kohn
  author author T. A. Wesolowski and author J. Weber, title title Kohn-Sham
  Equations with Constrained Electron Density: an Iterative Evaluation of the
  Ground-State Electron Density of Interacting Molecules, @noop 
  journal journal Chem. Phys. Lett. volume 248, pages 71–76 (year
  1996)NoStop
[Casida and Wesołowski(2004)]casida2004generalization
  author author M. E. Casida and author T. A. Wesołowski, title title
  Generalization of the Kohn–Sham Equations with Constrained Electron
  Density Formalism and its Time-Dependent Response Theory Formulation,  @noop  journal journal Int. J. Quant.
  Chem. volume 96, pages 577–588
  (year 2004)NoStop
[Höfener, Gomes, and  Visscher(2012)]hofener2012molecular
  author author S. Höfener, author A. S. P. Gomes,  and author L. Visscher, title title Molecular
  Properties via a Subsystem Density Functional Theory Formulation: A Common
  Framework for Electronic Embedding, @noop  journal
  journal J. Chem. Phys. volume 136,  pages 044104 (year 2012)NoStop
[Neugebauer(2010)]neugebauer2010chromophore
  author author J. Neugebauer, title title
  Chromophore-Specific Theoretical Spectroscopy: From Subsystem Density
  Functional Theory to Mode-Specific Vibrational Spectroscopy, @noop
   journal journal Phys, Rep. volume 489, pages 1–87 (year
  2010)NoStop
[Neugebauer(2009)]neugebauer2009calculation
  author author J. Neugebauer, title title On the
  Calculation of General Response Properties in Subsystem Density Functional
  Theory, @noop  journal journal J.
  Chem. Phys. volume 131, pages 084104
  (year 2009)NoStop
[Neugebauer(2007)]neugebauer2007couplings
  author author J. Neugebauer, title title Couplings
  between Electronic Transitions in a Subsystem Formulation of Time-Dependent
  Density Functional Theory, @noop  journal journal J. Chem. Phys. volume 126, pages 134116 (year 2007)NoStop
[Olsen, Aidas, and Kongsted(2010)]olsen2010
  author author J. M. H.  Olsen, author K. Aidas,  and author J. Kongsted, title title Excited
  States in Solution through Polarizable Embedding, @noop 
  journal journal J. Chem. Theory Comput.  volume 6, pages 3721–3734 (year 2010)NoStop
[Lipparini, Cappelli, and  Barone(2012)]lipparini2012linear
  author author F. Lipparini, author C. Cappelli,
   and author V. Barone,  title title Linear Response Theory and
  Electronic Transition Energies for a Fully Polarizable QM/Classical
  Hamiltonian, @noop  journal journal
  J. Chem. Theory Comput. volume 8, pages 4153–4165 (year 2012)NoStop
[Yoo et al.(2008)Yoo,
  Zahariev, Sok, and Gordon]EFP-excited-states
  author author S. Yoo, author F. Zahariev,
  author A. Sok,  and author M. S. Gordon, title
  title Solvent effects on optical properties of
  molecules: A combined time-dependent density functional theory/effective
  fragment potential approach, @noop  journal
  journal J. Chem. Phys. volume 129,  pages 144112 (year 2008)NoStop
[List, Olsen, and Kongsted(2016)]list2016excited
  author author N. H. List, author J. M. H. Olsen,
   and author J. Kongsted,  title title Excited States in Large
  Molecular Systems through Polarizable Embedding, @noop 
  journal journal Phys. Chem. Chem. Phys.  volume 18, pages 20234–20250 (year 2016)NoStop
[Jensen, van Duijnen, and  Snijders(2003)]jensen2003discrete2
  author author L. Jensen, author P. T. van
  Duijnen,  and author J. G. Snijders, title title A Discrete
  Solvent Reaction Field Model for Calculating Molecular Linear Response
  Properties in Solution, @noop  journal journal J. Chem. Phys. volume 119, pages 3800–3809 (year 2003)NoStop
[Morton and Jensen(2011)]morton2011discrete
  author author S. M. Morton and author L. Jensen, title title A Discrete
  Interaction Model/Quantum Mechanical Method to Describe the Interaction of
  Metal Nanoparticles and Molecular Absorption, @noop  journal journal J. Chem. Phys. volume
  135, pages 134103 (year 2011)NoStop
[Curutchet et al.(2009)Curutchet, Muñoz-Losa, Monti,
  Kongsted, Scholes, and Mennucci]curutchet2009electronic
  author author C. Curutchet, author A. Muñoz-Losa, author S. Monti, author J. Kongsted,
  author G. D. Scholes,  and  author B. Mennucci, title title Electronic Energy Transfer in Condensed
  Phase Studied by a Polarizable QM/MM Model, @noop  journal journal J. Chem. Theory Comput. volume 5, pages 1838–1848 (year
  2009)NoStop
[Jacob et al.(2006)Jacob,
  Neugebauer, Jensen, and Visscher]jacob2006comparison
  author author C. R. Jacob, author J. Neugebauer,
  author L. Jensen,  and author L. Visscher, title
  title Comparison of Frozen-Density Embedding and
  Discrete Reaction Field Solvent Models for Molecular Properties,  @noop  journal journal Phys. Chem. Chem.
  Phys. volume 8, pages 2349–2359
  (year 2006)NoStop
[Payton et al.(2012)Payton,
  Morton, Moore, and Jensen]payton2012discrete
  author author J. L. Payton, author S. M. Morton,
  author J. E. Moore,  and  author L. Jensen, title title A Discrete Interaction Model/Quantum
  Mechanical Method for Simulating Surface-Enhanced Raman Spectroscopy,  @noop  journal journal J. Chem. Phys.  volume 136, pages 214103 (year 2012)NoStop
[List et al.(2013)List,
  Jensen, Kongsted, and Hedegård]list2013unified
  author author N. H. List, author H. J. Aa. Jensen, author J. Kongsted,   and author E. D. Hedegård,  title title A Unified Framework for the
  Polarizable Embedding and Continuum Methods Within Multiconfigurational
  Self-Consistent Field Theory, @noop  journal
  journal Adv. Quantum Chem. volume
  66, pages 195 (year 2013)NoStop
[Olsen and Kongsted(2011)]olsen2011
  author author J. M. H.  Olsen and author J. Kongsted, title title Molecular Properties through Polarizable Embedding, @noop
   journal journal Adv. Quantum Chem. volume 61, pages 107–143 (year
  2011)NoStop
[Olsen et al.(2015b)Olsen, Steinmann,
  Ruud, and Kongsted]olsen2015polarizable
  author author J. M. H.  Olsen, author C. Steinmann, author K. Ruud,   and author J. Kongsted,  title title Polarizable Density
  Embedding: A New QM/QM/MM-based Computational Strategy, @noop 
  journal journal J. Phys. Chem. A volume 119, pages 5344–5355 (year
  2015b)NoStop
[Olsen and Jørgensen(1995)]olsen1995
  author author J. Olsen and author P. Jørgensen, @noop  title Modern Electronic
  Structure Theory, edited by editor D. R.  Yarkony, Vol. volume 2 (publisher
  World Scientific, year 1995) pp. pages
  857–990NoStop
[Corni et al.(2005)Corni,
  Cammi, Mennucci, and Tomasi]corni2005electronic1
  author author S. Corni, author R. Cammi,
  author B. Mennucci,  and  author J. Tomasi, title title Electronic Excitation Energies of
  Molecules in Solution within Continuum Solvation Models: Investigating the
  Discrepancy between State-Specific and Linear-Response Methods,  @noop  journal journal J. Chem. Phys.  volume 123, pages 134512 (year 2005)NoStop
[List(2015)]list2015thesis
  author author N. H. List, title Theoretical Description of Electronic
  Transitions in Large Molecular Systems in the Optical and X-Ray Regions,  @noop  Ph.D. thesis, school University of Southern
  Denmark, address Odense, Denmark (year 2015),  note
  <http://www.diva-portal.org/smash/record.jsf?pid=diva2:1072871&dswid=7352>NoStop
[Schwabe(2016)]schwabe2016general
  author author T. Schwabe, title title General Theory
  for Environmental Effects on (Vertical) Electronic Excitation Energies,  @noop  journal journal J. Chem. Phys.  volume 145, pages 154105 (year 2016)NoStop
[Rösch and Zerner(1994)]roesch1994calculation
  author author N. Rösch and author M. C. Zerner, title title Calculation of
  Dispersion Energy Shifts in Molecular Electronic Spectra, @noop 
  journal journal J. Phys. Chem. volume 98, pages 5817–5823 (year
  1994)NoStop
[McRae(1957)]mcrae1957theory
  author author E. G. McRae, title title Theory of
  Solvent Effects on Molecular Electronic Spectra. Frequency Shifts,  @noop  journal journal J. Phys. Chem.  volume 61, pages 562–572 (year 1957)NoStop
[Buhmann(2012)]buhmann2012dispersion
  author author S. Y. Buhmann, @noop  title Dispersion Forces
  II (publisher Springer Berlin Heidelberg, year
  2012)NoStop
[Lunkenheimer and Köhn(2012)]lunkenheimer2012solvent
  author author B. Lunkenheimer and author A. Köhn, title title Solvent
  Effects on Electronically Excited States Using the Conductor-Like Screening
  Model and the Second-Order Correlated Method ADC (2), @noop 
  journal journal J. Chem. Theory Comput.  volume 9, pages 977–994 (year 2012)NoStop
[Daday et al.(2015)Daday,
  Curutchet, Sinicropi, Mennucci, and Filippi]daday2015chromophore
  author author C. Daday, author C. Curutchet,
  author A. Sinicropi, author B. Mennucci,  and author C. Filippi, title
  title Chromophore–Protein Coupling Beyond
  Nonpolarizable Models: Understanding Absorption in Green Fluorescent
  Protein, @noop  journal journal
  J. Chem. Theory Comput. volume 11, pages 4825–4839 (year 2015)NoStop
[Jensen, Swart, and van Duijnen(2005)]jensen2005microscopic
  author author L. Jensen, author M. Swart,   and author P. T. van Duijnen,  title title Microscopic and Macroscopic
  Polarization within a Combined Quantum Mechanics and Molecular Mechanics
  Model, @noop  journal journal J.
  Chem. Phys. volume 122, pages 034103
  (year 2005)NoStop
[Pavanello(2013)]pavanello2013subsystem
  author author M. Pavanello, title title On the
  Subsystem Formulation of Linear-Response Time-Dependent DFT, @noop
   journal journal J. Chem. Phys. volume 138, pages 204118 (year
  2013)NoStop
[Ángyán(1992)]angyan1992common
  author author J. G. Ángyán, title title Common
  Theoretical Framework for Quantum Chemical Solvent Effect Theories,  @noop  journal journal J. Math. Chem.  volume 10, pages 93–137 (year 1992)NoStop
[McWeeny and Sutcliffe(1989)]mcweeny1969methods
  author author R. McWeeny and author B. T. Sutcliffe, @noop  title Methods of Molecular
  Quantum Mechanics, edition 2nd ed., Vol. volume 2 (publisher Academic press London, year 1989)NoStop
[Helgaker, Jørgensen, and  Olsen(2000)]helgaker2000molecular
  author author T. Helgaker, author P. Jørgensen,  and author J. Olsen, @noop  title Molecular
  Electronic-Structure Theory (publisher Wiley, year 2000)NoStop
[Christiansen, Jørgensen, and  Hättig(1998)]christiansen1998response
  author author O. Christiansen, author P. Jørgensen,  and author C. Hättig, title title Response
  Functions from Fourier Component Variational Perturbation Theory applied to a
  Time-Averaged Quasienergy, @noop  journal journal Int. J. Quant. Chem. volume 68,  pages 1–52 (year 1998)NoStop
[Sasagane, Aiga, and Itoh(1993)]sasagane1993
  author author K. Sasagane, author F. Aiga,   and author R. Itoh, title title Higher–Order Response Theory Based on
  the Quasienergy Derivatives: The Derivation of the Frequency–Dependent
  Polarizabilities and Hyperpolarizabilities, http://dx.doi.org/10.1063/1.466123 journal journal J. Chem. Phys. volume 99, pages 3738–3778 (year 1993)NoStop
[Norman(2011)]norman2011perspective
  author author P. Norman, title title A Perspective
  on Nonresonant and Resonant Electronic Response Theory for Time-Dependent
  Molecular Properties, @noop  journal journal Phys. Chem. Chem. Phys. volume 13,  pages 20519–20535 (year 2011)NoStop
[Olsen and Jørgensen(1985)]olsen1985linear
  author author J. Olsen and author P. Jørgensen, title title Linear
  and Nonlinear Response Functions for an Exact State and for an MCSCF
  State, @noop  journal journal J.
  Chem. Phys. volume 82, pages
  3235–3264 (year 1985)NoStop
[List et al.(2014)List,
  Coriani, Christiansen, and Kongsted]list2014identifying
  author author N. H. List, author S. Coriani,
  author O. Christiansen,  and  author J. Kongsted, title title Identifying the Hamiltonian Structure
  in Linear Response Theory, @noop  journal journal J. Chem. Phys. volume 140, pages 224103 (year 2014)NoStop
[Löwdin(1962a)]lowdin1962bstudies
  author author P.-O. Löwdin, title title Studies in
  Perturbation Theory. V. Some Aspects on the Exact Self-Consistent Field
  Theory, @noop  journal journal J.
  Math. Phys. volume 3, pages
  1171–1184 (year 1962a)NoStop
[Löwdin(1963)]lowdin1963studies
  author author P.-O. Löwdin, title title Studies in
  Perturbation Theory: Part I. An Elementary Iteration-Variation Procedure for
  Solving the Schrödinger Equation by Partitioning Technique,  @noop  journal journal J. Mol.
  Spectrosc. volume 10, pages 12–33
  (year 1963)NoStop
[Löwdin(1962b)]lowdin1962astudies
  author author P.-O. Löwdin, title title Studies in
  Perturbation Theory. IV. Solution of Eigenvalue Problem by Projection
  Operator Formalism, @noop  journal journal J. Math. Phys. volume 3, pages 969–982 (year 1962b)NoStop
[Hsu et al.(2001)Hsu,
  Fleming, Head-Gordon, and Head-Gordon]hsu2001excitation
  author author C.-P. Hsu, author G. R. Fleming,
  author M. Head-Gordon,  and  author T. Head-Gordon,  title title Excitation Energy Transfer
  in Condensed Media, @noop  journal journal J. Chem. Phys. volume 114, pages 3065–3072 (year 2001)NoStop
[Rinaldi, Morton, and Jensen(2013)]rinaldi2013discrete
  author author J. M. Rinaldi, author S. M. Morton,
   and author L. Jensen,  title title A Discrete Interaction
  Model/Quantum Mechanical Method for Simulating Nonlinear Optical Properties
  of Molecules Near Metal Surfaces, @noop  journal
  journal Mol. Phys. volume 111,  pages 1322–1331 (year 2013)NoStop
[Böttcher and Bordewijk(1973)]botcher1973theory
  author author C. J. F.  Böttcher and author P. Bordewijk, @noop  title
  Theory of Electric Polarization, Vol. 1: Dielectrics in Static Fields  (publisher Elsevier Amsterdam, year
  1973)NoStop
[Wortmann and Bishop(1998)]wortmann1998effective
  author author R. Wortmann and author D. M. Bishop, title title Effective
  Polarizabilities and Local Field Corrections for Nonlinear Optical
  Experiments in Condensed Media, @noop  journal
  journal J. Chem. Phys. volume 108,  pages 1001–1007 (year 1998)NoStop
[Cammi, Mennucci, and Tomasi(1998)]cammi1998calculation
  author author R. Cammi, author B. Mennucci,   and author J. Tomasi,  title title On the Calculation of Local
  Field Factors for Microscopic Static Hyperpolarizabilities of Molecules in
  Solution with the Aid of Quantum-Mechanical Methods, @noop 
  journal journal J. Phys. Chem. A volume 102, pages 870–875 (year
  1998)NoStop
[Pipolo, Corni, and Cammi(2014)]pipolo2014cavity2
  author author S. Pipolo, author S. Corni,   and author R. Cammi, title title The Cavity Electromagnetic Field
  within the Polarizable Continuum Model of Solvation, @noop 
  journal journal J. Chem. Phys. volume 140, pages 164114 (year
  2014)NoStop
[List, Jensen, and Kongsted(2016)]list2016local
  author author N. H. List, author H. J. Aa. Jensen,  and author J. Kongsted, title title Local
  Electric Fields and Molecular Properties in Heterogeneous Environments
  through Polarizable Embedding, DOI: 10.1039/C6CP00669H
  journal journal Phys. Chem. Chem. Phys.  volume 18, pages 10070–10080 (year 2016)NoStop
[Nielsen, Jørgensen, and  Oddershede(1980)]nielsen1980transition
  author author E. S. Nielsen, author P. Jørgensen,  and author J. Oddershede, title title Transition
  Moments and Dynamic Polarizabilities in a Second Order Polarization
  Propagator Approach, @noop  journal journal J. Chem. Phys. volume 73, pages 6238–6246 (year 1980)NoStop
[Dreuw and Wormit(2015)]dreuw2015algebraic
  author author A. Dreuw and author M. Wormit,  title title The Algebraic Diagrammatic
  Construction Scheme for the Polarization Propagator for the Calculation of
  Excited States, @noop  journal journal
  WIREs Comput. Mol. Sci. volume 5, pages 82–95 (year 2015)NoStop
[Norman et al.(2005)Norman,
  Bishop, Jensen, and Oddershede]norman2005nonlinear
  author author P. Norman, author D. M. Bishop,
  author H. J. Aa. Jensen,  and author J. Oddershede, title title Nonlinear
  Response Theory with Relaxation: The First-Order Hyperpolarizability,  @noop  journal journal J. Chem. Phys.  volume 123, pages 194103 (year 2005)NoStop
[Norman et al.(2001)Norman,
  Bishop, Jensen, and Oddershede]norman2001near
  author author P. Norman, author D. M. Bishop,
  author H. J. Aa. Jensen,  and author J. Oddershede, title title
  Near-Resonant Absorption in the Time-Dependent Self-Consistent Field and
  Multiconfigurational Self-Consistent Field Approximations, @noop 
  journal journal J. Chem. Phys. volume 115, pages 10323–10334 (year
  2001)NoStop
[Kauczor, Jørgensen, and  Norman(2011)]kauczor2011onthe
  author author J. Kauczor, author P. Jørgensen,  and author P. Norman, title title On the
  Efficiency of Algorithms for Solving Hartree-Fock and Kohn-Sham Response
  Equations, 10.1021/ct100729t journal
  journal J. Chem. Theory Comput. volume
  7, pages 1610–1630 (year 2011)NoStop
[Boyd(2003)]boyd2003nonlinear
  author author R. W. Boyd, @noop  title Nonlinear Optics  (publisher Academic press, year 2003)NoStop
[List et al.(2015)List,
  Kauczor, Saue, Jensen, and  Norman]list2015beyond
  author author N. H. List, author J. Kauczor,
  author T. Saue, author
  H. J. Aa. Jensen,  and author P. Norman, title
  title Beyond the Electric-Dipole Approximation: A
  Formulation and Implementation of Molecular Response Theory for the
  Description of Absorption of Electromagnetic Field Radiation, @noop
   journal journal J. Chem. Phys. volume 142, pages 244111 (year
  2015)NoStop
[Stone(1989)]stone1989induction
  author author A. J. Stone, title title The Induction
  Energy of an Assembly of Polarizable Molecules, @noop  journal journal Chem. Phys. Lett. volume 155, pages 102–110 (year
  1989)NoStop
[Saint Raymond(1991)]saint1991elementary
  author author X. Saint Raymond, @noop  title Elementary
  Introduction to the Theory of Pseudodifferential Operators,  Vol. volume 3 (publisher CRC Press, year 1991)NoStop
[Olsen(2012)]olsen2012thesis
  author author J. M. H.  Olsen, title Development of Quantum
  Chemical Methods towards Rationalization and Optimal Design of Photoactive
  Proteins, @noop  Ph.D. thesis, school University
  of Southern Denmark, address Odense, Denmark (year
  2012), note DOI: 10.6084/m9.figshare.156852NoStop
[Stone(2002)]stone2013theory
  author author A. Stone, @noop  title The Theory of
  Intermolecular Forces (publisher Oxford University Press,  year 2002)NoStop
[Buckingham(1967)]buckingham1967
  author author A. D. Buckingham, title title Permanent
  and Induced Molecular Moments and Long-Range Intermolecular Forces,  @noop  journal journal Adv. Chem. Phys.  volume 12, pages 107–142 (year 1967)NoStop
[Applequist, Carl, and  Fung(1972)]applequist1972atom
  author author J. Applequist, author J. R. Carl,  and author K.-K. Fung, title title Atom Dipole
  Interaction Model for Molecular Polarizability. Application to Polyatomic
  Molecules and Determination of Atom Polarizabilities, @noop 
  journal journal J. Am. Chem. Soc. volume 94, pages 2952–2960 (year
  1972)NoStop
[Harczuk, Vahtras, and  Ågren(2015)]harczuk2015frequency
  author author I. Harczuk, author O. Vahtras,   and author H. Ågren,  title title Frequency-Dependent Force
  Fields for QMMM Calculations, @noop  journal
  journal Phys. Chem. Chem. Phys. volume
  17, pages 7800–7812 (year 2015)NoStop
[Nørby et al.(2016)Nørby, Vahtras, Norman, and  Kongsted]norby2016assessing
  author author M. S. Nørby, author O. Vahtras,
  author P. Norman,  and author J. Kongsted, title
  title Assessing Frequency-Dependent Site
  Polarisabilities in Linear Response Polarisable Embedding, @noop 
  journal journal Mol. Phys. , pages
  1–9 (year 2016)NoStop
[Yanai, Tew, and Handy(2004)]yanai2004new
  author author T. Yanai, author D. P. Tew,   and author N. C. Handy,  title title A New Hybrid
  Exchange–Correlation Functional using the Coulomb-Attenuating Method
  (CAM-B3LYP), @noop  journal journal
  Chem. Phys. Lett. volume 393, pages
  51–57 (year 2004)NoStop
[Dunning(1989)]dunning
  author author T. H. Dunning, title title Gaussian Basis
  Sets for Use in Correlated Molecular Calculations. I. The Atoms Boron through
  Neon and Hydrogen, @noop  journal journal J. Chem. Phys. volume 90, pages 1007–1023 (year 1989)NoStop
[Gagliardi, Lindh, and  Karlström(2004)]loprop
  author author L. Gagliardi, author R. Lindh,   and author G. Karlström,  title title Local Properties of Quantum
  Chemical Systems: The LoProp Approach. @noop  journal journal J. Chem. Phys. volume
  121, pages 4494–4500 (year 2004)NoStop
[Vahtras(2014)]Vahtras:13276
  author author O. Vahtras, @noop  title LoProp for
  Dalton, howpublished see
  <http://dx.doi.org/10.5281/zenodo.13276> (year
  2014)NoStop
[Aidas et al.(2014)Aidas,
  Angeli, Bak, Bakken,
  Bast, Boman, Christiansen,
  Cimiraglia, Coriani, Dahle,
  Dalskov, Ekström, Enevoldsen, Eriksen, Ettenhuber,
  Fernández, Ferrighi, Fliegl, Frediani, Hald, Halkier, Hättig, Heiberg, Helgaker, Hennum, Hettema, Hjertenæs, Høst, Høyvik,
  Iozzi, Jansik, Jensen,
  Jonsson, Jørgensen, Kauczor, Kirpekar, Kjærgaard,
  Klopper, Knecht, Kobayashi,
  Koch, Kongsted, Krapp,
  Kristensen, Ligabue, Lutnæs, Melo, Mikkelsen, Myhre, Neiss, Nielsen, Norman, Olsen, Olsen, Osted, Packer, Pawlowski, Pedersen, Provasi, Reine, Rinkevicius, Ruden, Ruud, Rybkin, Sałek, Samson, de Merás, Saue, Sauer, Schimmelpfennig, Sneskov, Steindal,
  Sylvester-Hvid, Taylor, Teale, Tellgren, Tew, Thorvaldsen, Thøgersen, Vahtras,
  Watson, Wilson, Ziolkowski,  and Ågren]dalton
  author author K. Aidas, author C. Angeli,
  author K. L. Bak, author V. Bakken, author
  R. Bast, author L. Boman, author O. Christiansen, author R. Cimiraglia, author S. Coriani, author P. Dahle, author E. K.  Dalskov, author U. Ekström, author T. Enevoldsen, author J. J. Eriksen, author P. Ettenhuber,
  author B. Fernández, author L. Ferrighi, author
  H. Fliegl, author L. Frediani, author K. Hald, author A. Halkier, author C. Hättig,
  author H. Heiberg, author T. Helgaker, author
  A. C. Hennum, author
  H. Hettema, author E. Hjertenæs, author S. Høst, author I.-M.  Høyvik, author M. F.  Iozzi, author B. Jansik, author H. J. Aa.  Jensen, author D. Jonsson, author P. Jørgensen, author J. Kauczor, author S. Kirpekar,
  author T. Kjærgaard, author W. Klopper, author
  S. Knecht, author R. Kobayashi, author H. Koch, author J. Kongsted, author A. Krapp,
  author K. Kristensen, author A. Ligabue, author
  O. B. Lutnæs, author
  J. I. Melo, author
  K. V. Mikkelsen, author
  R. H. Myhre, author
  C. Neiss, author C. B. Nielsen, author P. Norman, author J. Olsen, author J. M. H.  Olsen, author A. Osted, author M. J. Packer,
  author F. Pawlowski, author T. B. Pedersen, author
  P. F. Provasi, author
  S. Reine, author Z. Rinkevicius, author T. A.  Ruden, author K. Ruud, author V. Rybkin,
  author P. Sałek, author C. C. M. Samson, author
  A. S. de Merás, author
  T. Saue, author S. P. A. Sauer, author B. Schimmelpfennig, author K. Sneskov, author A. H.  Steindal, author K. O.  Sylvester-Hvid, author P. R.  Taylor, author A. M.  Teale, author E. I. Tellgren, author D. P. Tew,
  author A. J. Thorvaldsen,
  author L. Thøgersen, author O. Vahtras, author
  M. A. Watson, author
  D. J. D. Wilson, author
  M. Ziolkowski,  and author
  H. Ågren, title title The Dalton Quantum Chemistry Program System,  @noop  journal journal WIREs
  Comput. Mol. Sci. volume 4, pages
  269–284 (year 2014)NoStop
[dal()]dalton2
  @noop  title Dalton, A Molecular Electronic
  Structure Program, release dalton2016.0 (2016), howpublished
  see <http://daltonprogram.org/>NoStop
[Steinmann and Kongsted(2015)]steinmann2015electronic
  author author C. Steinmann and author J. Kongsted, title title Electronic
  Energy Transfer in Polarizable Heterogeneous Environments: A Systematic
  Investigation of Different Quantum Chemical Approaches, @noop 
  journal journal J. Chem. Theory Comput.  volume 11, pages 4283–4293 (year 2015)NoStop
[Maitra et al.(2004)Maitra,
  Zhang, Cave, and Burke]maitra2004double
  author author N. T. Maitra, author F. Zhang,
  author R. J. Cave,  and  author K. Burke, title title Double Excitations Within
  Time-Dependent Density Functional Theory Linear Response, @noop 
  journal journal J. Chem. Phys. volume 120, pages 5932–5937 (year
  2004)NoStop
[Cave et al.(2004)Cave,
  Zhang, Maitra, and Burke]cave2004dressed
  author author R. J. Cave, author F. Zhang,
  author N. T. Maitra,  and  author K. Burke, title title A Dressed TDDFT Treatment of the 2 1 A
  g States of Butadiene and Hexatriene, @noop  journal journal Chem. Phys. Lett. volume 389, pages 39–42 (year
  2004)NoStop
[Casida and Huix-Rotllant(2015)]casida2015many
  author author M. E. Casida and author M. Huix-Rotllant, title title
  Many-Body Perturbation Theory (MBPT) and Time-Dependent Density-Functional
  Theory (TD-DFT): MBPT Insights About What Is Missing In, and Corrections To,
  the TD-DFT Adiabatic Approximation, in @noop  booktitle Density-Functional Methods for Excited States (publisher Springer, year 2015) pp. pages
  1–60NoStop
[Caballero, Moreira, and  Bofill(2013)]caballero2013comparison
  author author M. Caballero, author I. d. P. Moreira,  and author J. M. Bofill, title title A Comparison
  Model between Density Functional and Wave Function Theories by Means of the
  Löwdin Partitioning Technique, @noop  journal
  journal J. Chem. Phys. volume 138,  pages 174107 (year 2013)NoStop




