Open-ended formulation of self-consistent field response theory with the polarizable continuum model for solvation

The study of high-order absorption properties of molecules is a field of growing importance. Quantumchemical studies can help design chromophores with desirable characteristics. Given that most experiments are performed in solution, it is important to devise a cost-effective strategy to include solvation effects in quantum-chemical studies of these properties. We here present an open-ended formulation of selfconsistent field (SCF) response theory for a molecular solute coupled to a polarizable continuum model (PCM) description of the solvent. Our formulation relies on the open-ended, density matrix-based quasienergy formulation of SCF response theory of Thorvaldsen, et al., [J. Chem. Phys., 2008, 129, 214108] and the variational formulation of the PCM, as presented by Lipparini et al., [J. Chem. Phys., 2010, 133, 014106]. Within the PCM approach to solvation, the mutual solute–solvent polarization is represented by means of an apparent surface charge (ASC) spread over the molecular cavity defining the solute–solvent boundary. In the variational formulation, the ASC is an independent, variational degree of freedom. This allows us to formulate response theory for molecular solutes in the fixed-cavity approximation up to arbitrary order and with arbitrary perturbation operators. For electric dipole perturbations, pole and residue analyses of the response functions naturally lead to the identification of excitation energies and transition moments. We document the implementation of this approach in the Dalton program package using a recently developed open-ended response code and the PCMSolver libraries and present results for one-, two-, three-, fourand five-photon absorption processes of three small molecules in solution.


Introduction
An important challenge in molecular sciences is the study at the quantum molecular mechanical level of systems of growing complexity. The behavior of simple, isolated molecules is in general well understood, but large systems that include multiple constituents pose additional problems, both due to their sheer size and the interaction of the different components. On the experimental side, the best tools to address such systems are provided by spectroscopic techniques, in which photons interact with the system and the response of the system to these perturbations is monitored. Simultaneously, technological improvements both in the intensity of the photon source and the sensitivity of detectors allow new experimental techniques to be developed and applied to complex systems. As a result, methods that used to be proof-of-principle concepts are today routinely employed to study increasingly complex systems.
An example of such a class of experiments is multiphoton absorption (MPA): the simultaneous absorption of several photons. The effect, originally predicted for two-photon absorption (2PA) by Göppert-Mayer in 1931, 1 is too weak to be detected unless laser sources are employed. As a consequence, the first measurement was only made possible in 1961. 2 2PA is still not as widespread as one-photon absorption (1PA), but its use is increasing. In part, this is due to its different symmetry selection rules that allow exploring excited states that are dark in 1PA. In addition, 2PA experiments afford a greater focality than in 1PA, due to their quadratic dependence on the intensity of the incident radiation. More recently, higher-order methods such as three-, 3 four-4 and fivephoton 5 absorption have emerged, although their use is in no way as widespread as 1PA and 2PA.
The growth in available spectroscopic techniques and their application to systems of ever increasing complexity calls for a corresponding effort on the theoretical side to describe correctly the underlying molecular phenomena and to tackle the complexity of the system in a manageable way.
Olsen and Jørgensen have shown how transition moments between ground and excited states can be calculated from residues of response functions of a molecule in its ground state. 6 The residues of the linear response function can thus be used to calculate the strength of transitions in UV/Vis spectroscopy (1PA). 2PA and three-photon absorption (3PA) can in turn be calculated from the residues of the quadratic and cubic response functions, respectively. [6][7][8] Similarly, four-and five-photon absorption cross sections (4PA and 5PA) can be calculated from the corresponding higher-order response functions. The higher the order of the response functions needed, the more complex the working equations become, in particular if attention is given to computational efficiency.
Our group has in the last few years developed an open-ended formulation of response theory 9 and implemented it using recursive programming techniques, 10 enabling the calculation of response properties to any order at the Hartree-Fock (HF) and density-functional theory (DFT) levels and limited only by the degree of generality in connected modules for perturbed one-and two-electron integrals and exchange-correlation (XC) contributions. [9][10][11] The approach has recently been extended to include single residues of response functions. 11 Single residues of these high-order response functions have already been used by our group to calculate 4PA 11 and 5PA 12 absorption cross sections. Arbitrary higher-order processes are also accessible from this open-ended approach. The open-ended response formalism is therefore able to address the challenge of the ever-growing variety of spectroscopic methods available, significantly reducing the development effort and the time required to model new spectroscopic processes for relevant applications.
Several approaches have been developed to tackle large and complex systems, such as solutions and proteins, in the presence of external fields. When the phenomenon studied is localized to a single molecule and its immediate surroundings-as is often the case for MPA where the majority of the response arises from a chromophore in the complex-an efficient strategy is to use focused models that only treat a small portion of the system (e.g. the chromophore molecule) using quantum mechanical (QM) methods, whereas the rest (the environment) is treated classically. A distinction can be drawn between methods keeping atomistic detail in the classical environment, and those which disregard it. The former are commonly known as molecular mechanics (MM) methods, whereas the latter are referred to as dielectric continuum (DC) methods. Both models have strengths and weaknesses: MM methods are better suited to describe specific intermolecular interactions but require configurational sampling.
In contrast, DC methods are more effective at addressing long-range interactions. [13][14][15] and are easier to parametrize, requiring only a handful of macroscopic solvent parameters (dielectric constant e, refractive index n), and the cavity specification (atom-specific radii). Both methods can be augmented with a supermolecular approach (including one or more solvent molecules in the QM system) to address specific quantum effects. A three-layer model which combines a description of the specific effects near the chromophore by an MM part with pre-averaged long-range effects described by the DC part has also been reported. 16 For details we refer to the relevant literature. [16][17][18][19] The above mentioned open-ended approach to response theory has recently been extended to include molecular environment effects for electric dipole properties through a polarizable embedding (PE) QM/MM approach. 20 In this work, we will present an open-ended response formalism for the PCM, 21 in its integral equation formalism (IEF) formulation, 22 which is the most versatile DC method available. For details about the PCM, the reader is referred to two authoritative reviews. 13,23 The model features a molecule-shaped cavity made of interlocking spheres, 24,25 is able to describe a wide variety of environments due to the generality of the IEF formalism, 22,26,27 and can treat dynamical processes thanks to the nonequilibrium formalism. 28,29 All such features are available in the PCMSolver module, an application programming interface (API) for the PCM. 30 Additional features not yet available in PCMSolver are the treatment of nonelectrostatic terms in the solvation energy, 31,32 and the state-specific formalism. [33][34][35][36] Crucial aspects of our work are the variational formulation of the PCM equations 37 and the modular approach employed in the implementation. Both PCMSolver and the open-ended response code 10 are two independent modules which can be interfaced to any quantum chemistry software. This approach has several advantages over a non-modular code: modules can be developed and tested separately, new features can be made available to several programs at once, avoiding lengthy, tedious and error-prone multiple implementations, and the master program can be chosen freely, for instance based on the availability of different functionality.
The rest of the paper is organized as follows: in Section 2 we present the theory for the quasienergy formalism in the context of the PCM. In Section 3 we discuss the details of our modular implementation. After briefly discussing the computational details (Section 4), we will present our results on the MPA processes (up to 5PA) on para-nitroaniline (PNA), para-dinitrobenzene (PDNB) and methylenecyclopropene (MCP) in Section 5. In Section 6 we summarize the main conclusions and future implications of our work.

Variational formulation of the polarizable continuum model
The variational formulation of the PCM was first presented by Lipparini et al. in ref. 37 and is based on the weak approach to boundary integral equations (BIEs). 38 The weak formulation of partial differential equations (PDEs), boundary value problems (BVPs) associated BIEs is a well-known tool in mathematics. [38][39][40] Given a partition of Euclidean space R 3 into a closed subdomain C, the cavity, with a sufficiently regular boundary G = qC, we want to solve the following transmission problem for a solvent with a homogeneous, isotropic relative permittivity e: where n is the outward-pointing normal vector to the cavity boundary G. The electrostatic potential u(r) in space is sought, given the jump conditions for its traces and conormal derivatives across the boundary, eqn (1c) and (1d), respectively, and the appropriate radiation condition at infinity, eqn (1e). This can be recast in terms of an integral equation: with s(s), the apparent surface charge (ASC), representing the reaction potential arising from solvent polarization and j(s) the molecular electrostatic potential (MEP). The integral operatorŝ R e andR 1 are given in terms of the components of the Calderón projector,Ŝ andD, 38,41 and the identity operatorÎ: such that the operatorŶ ¼R 1 À1R eŜ is self-adjoint and positive definite. TheŜ andD boundary integral operators are mappings between Sobolev spaces of fractional order, which thus are the natural mathematical setting for integral formulations of BVPs. [38][39][40] These are normed spaces, equipped with the scalar product: ð The polarization energy functional: is strictly convex and has a unique minimum, s 0 . This is the unique solution to the IEF-PCM eqn (2): 37,42 @U pol @s ¼Ŷs þ j ¼ 0 This allows us to treat the ASC as an additional, independent, variational density to be optimized. This offers distinct advantages from a theoretical point of view: there is no need to invoke a nonlinear coupling in the Hamiltonian to introduce the classical solute-solvent interaction, 13,43 the functional clearly describes a charge distribution interacting (unfavorably) with itself and (favorably) with its inducing external field and constitutes the polarization energy of the medium, 37 a classical analogue of the Hellmann-Feynman theorem naturally holds for the variational ASC: 44 Simultaneous optimization algorithms can also be successfully employed in practical implementations, 45 but this is not the main topic of this work. Finally, let us note that the use of the term weak formulation of PDEs and BVPs originates from the weaker regularity requirements that can be imposed on the solution, while still handling a well-posed problem (in the sense of Hadamard). The terms ''weak'' and ''variational'' formulation are here used interchangeably, given that the weak formulation of the PCM satisfies the hypotheses of both the Lax-Milgram lemma and its variational corollary. 39

PCM-SCF open-ended response theory
Notation. The PCM equations will be written in the ''complete basis'': we will introduce the usual boundary-element method (BEM) discretization at the very end of the derivation. In other words, we will be working with the exact integral equation and not with its discretized counterpart. As a consequence, the apparent surface charge s and the electrostatic potential j will have a continuous dependence on a ''cavity surface'' index s. Whenever a charge-potential product is present, it is to be interpreted as the surface integral, i.e. the scalar product in the suitable, infinite-dimensional vector space on the cavity boundary G. The following shorthand notation will be adopted: sj = (s, j) G . We use lowercase Latin letters (a, b, c. . .) as a composite index for the perturbation operator and the frequency index (cf. eqn (7)-(16) in ref. 9). The perturbation strength for a given perturbing one-electron operator A associated with a frequency o a will thus be written as e a . Perturbationstrength derivatives will be denoted by lowercase Latin superscripts (a, b, c. . .) to the differentiated quantities. Finally, a tilde will be used for quantities that are considered at general field strengths and thus, in general, are time dependent. As an example, the overlap matrix and its derivative with respect to e a at general perturbation strength will be S and S a , respectively. Equivalently, S and S a denote the overlap matrix and its perturbation-strength derivative at zero field strength, respectively. Tr ¼ will denote that the trace of the expression to follow should be taken. fTrg T ¼ will additionally denote that the tracing is followed by timeaveraging over a period T of the collected perturbations.
Our derivation follows closely the one in ref. 9 and the subsequent developments in ref. 10. The original expressions were developed for a system considered to be in vacuo, and in order to incorporate the effects of the PCM, any energy-like term that appears in these expressions will be augmented by the appropriate solvent term. The solvent term will be derived according to the polarization energy functional given in eqn (5) and the classical Hellmann-Feynman theorem it satisfies, namely eqn (7).
Response functions can be expressed as perturbation-strength derivatives of the perturbation-strength-differentiated timeaveraged quasienergy Lagrangian evaluated at zero perturbation strengths. For example, the linear response function can be written as In an atomic orbital-based density matrix parametrization, the time-averaged quasienergy derivative needed to evaluate response functions is given as L a ðD;s; tÞ ¼ fTrg TG00;a ÀS aW (9) where an element of the overlap matrix S is given by and where the generalized, energy-weighted density matrix W was introduced and is given bỹ This expression for W involves the density matrix D, its timedifferentiated analogue _ D and the generalized Kohn-Sham (KS) matrixF given bỹ F ¼h þṼ t þG g ðDÞ þF xc þsũ À i 2T : The expression forF includes both vacuum-like and PCM contributions. The vacuum-like contributions are expressed in terms of the one-electron matrices h and Ṽ t , and the two-electron matrix G g (D), which are, respectively, defined in the following way: Another part of the vacuum-like contribution is the functional derivative matrix F xc,mn of the XC potential, whose elements are given byF where the integration involves the overlap distributionÕ mn 1 w Ã mwn and the functional derivative of the XC functional in the adiabatic approximation. The x variable refers to both spatial and spin coordinates. The last vacuum-like contribution in eqn (12) is the anti-Hermitian, time-differentiated overlap matrix T whose elements are given bỹ Finally, the PCM contributionsũ involves the electrostatic potential integralsj mn ðsÞ ¼w m À1 r À s j j wn : The first term in eqn (9),G 00;a , involves the generalized KS energyẼ as shown in eqn (97) in ref. 9. The free energy termG including PCM effects is produced by additionally considering solute-solvent interaction terms, so that Due to the implicit time dependence of D and s, higher-order derivatives of the KS generalized energy will require application of the chain rule. The mn, abc. . . superscript describes how and to what extent the chain rule was applied for a given term, i.e. the number of explicit differentiations with respect to the variational densities, so that In this notation, the index m denotes the order of differentiation with respect to the density matrix D, while the index n symbolizes the order of differentiation with respect to the ASC density s. Differentiation with respect to the density matrix will result in a 2m-rank tensor, while differentiation with respect to the ASC density will result in a function of the continuous cavity index s. For higher-order properties, mixed terms involving both density matrix and ASC density differentiation may generally occur. In the fixed-cavity approximation, the cavity is kept frozen at a given molecular geometry. 46 Under this simplifying assumption, only the linear interaction term in the polarization functional eqn (5) will be affected by the movements of the nuclei via the dependence of basis functions on the molecular geometry. Its perturbation-strength derivative will then be where the second term only involves derivatives of the electrostatic potential integrals. We remark that, under the fixed-cavity approximation, both the density matrix -m -and ASC densityn -differentiation indices in eqn (18) can only assume the values 0 or 1 in order for the @ mþnþ3 U pol @ D T ð Þ m @s n @e a @e b @e c term not to be zero. By construction, the density matrix dependence in the polarization functional is at most linear, while, by virtue of the classical Hellmann-Feynman theorem, eqn (7), the ASC variational density will also appear at most linearly inG 00;a . The free energy term perturbation strength derivative is given as G 00;a ðD;s;tÞ ¼Ẽ 0;a þsTrũ aD À Á where the matrixF Response functions can then be obtained by straightforward differentiation with respect to additional perturbations and subsequent evaluation at zero perturbation strength, so that and similarly for higher-order response functions. More detailed expressions for the derivatives of the generalized KS free energy are shown in Appendix A. The expressions (22a)-(22c) adhere to the n + 1 formulation, whereby perturbationstrength derivatives of the variational densities up to order n are required in order to assemble response functions of order n + 1. It is possible to make other formulations of response theory for which truncation rules for perturbed arguments between and including the n + 1 and 2n + 1 rules are possible. 9,10,47 This entails the introduction of Lagrange multipliersk a andf a to take into consideration the idempotency of the density matrix and the time-dependent SCF (TD-SCF) equations, respectively, so that the idempotency condition is expressed with the matrix Ỹ and the TD-SCF condition with the matrix Z, where and where the Lagrange multiplier terms are given bỹ andf a ¼F aDS À where the subscript integers k and n in the various forms shown in this expression denote a given choice of truncation rule. The original expression for systems considered in vacuo contains an energy term E abc... k;n instead of the free energy term G abc k,n but is otherwise unchanged upon solvation, and we will therefore omit further details here about the derivation leading up to eqn (27), referring instead to previous work for more information and for details about the (k, n) truncation rules that can be chosen and applied. 9,10 We note that the task of evaluating eqn (27) and obtaining terms needed for this evaluation can be cast in recursive form, as shown in ref. 10, and we further remark that these routines can be augmented to enable the calculation of single residues of response functions. 11 However, the methodological and algorithmic development needed for residues calculations is not changed by the inclusion of PCM effects, and we will therefore again refer to previous work 9,11 for details.

Parametrization of the perturbed densities and response equations
In order to compute response properties from eqn (27), the various perturbed D, F and S matrices and the derivatives of the ASC density s that enter into this expression must be obtained. The perturbed overlap matrices can be directly assembled from the relevant one-electron integral derivatives, while the perturbed density and Fock matrices can be obtained from a procedure that involves solving the appropriate response equations. The first step in this procedure is to take perturbationstrength derivatives of the idempotency and TD-SCF conditions of eqn (23) and (24) and evaluating them at zero perturbation strength. 9,10 The evaluation of the perturbation-strengthdifferentiated ASC density introduces an additional response equation, which is constructed by differentiating the equation governing the ASC:Ŷ s þ j ¼ 0 Differentiating eqn (23) and introducing a decomposition of the density matrix into frequency components leads to where b N is the tuple of applied perturbations and o is the sum of the associated frequencies. The RHS matrix K ðnÀ1Þ contains all terms that contain derivatives of the density matrix up to order n À 1.
The perturbed density matrix is partitioned into a particular D b N P and a homogenous D b N P term (H/P partition) as The former may be evaluated in terms of K (nÀ1) o , i.e. lower-order density matrices and differentiated overlap integrals, so that where the projectors P = DS, Q = 1 À P were used. The homogeneous component is parametrized in terms of the nth order response parameters X bN as The governing equations for the perturbed ASC densities are obtained in analogy with the handling of perturbed density matrices outlined above. We introduce a decomposition of the ASC density into frequency components into the perturbationstrength derivative of eqn (28), so that The symbol F (nÀ1) o has been introduced in analogy to the and it contains all terms that depend on lower-order density matrices and differentiated electrostatic potential integrals, for which the latter acts as the metric matrix S in the definition of o always contains at least a first derivative of the electrostatic potential integrals and is thus zero if the basis set is independent of the perturbation tuple being considered. We now introduce the H/P partitioning of s b N o , so that and apply eqn (30), which leads to a separation of the response integral equation into the following system of equations: We note that the particular ASC is nonzero if and only if the basis set depends on the external perturbation.
We finally turn our attention to the TDSCF equation. The perturbation-strength differentiated generalized KS matrix is first separated into its frequency components F b N o . The H/P partition introduced for the variational densities will induce a similar partition into these frequency components: The two-electron and XC contributions depending on the homogeneous perturbed density matrix have been collected A more detailed discussion of these aspects can be found in ref. 9 and 10. The parametrization of the homogeneous part of the perturbed density matrix can be exploited to conveniently reformulate the perturbed TDSCF equation, so that where the generalized Hessian E [2] and metric S [2] matrices were introduced and are defined by their transformations on the response parameters X bN : 48,49 The generalized Hessian matrix E [2] includes two types of solvent contributions: implicit terms included in the zeroth-order Fock matrix, F, and explicit terms, involving the N-th order homogenous ASC variational density. The latter are the last two terms in eqn (39). The theoretical treatment of frequency-dependent properties in solution within the PCM requires adoption of a nonequilibrium response framework. 13,50 The explicit PCM terms in eqn (39) are then evaluated using the optical permittivity, defined as the square of the refractive index of the solvent e N = n 2 , instead of the static permittivity e s , which is employed to compute the implicit contributions. In contrast to E [2] , the generalized metric matrix S [2] is unchanged. The right-hand side (RHS) in the response equation only includes terms that depend on particular contributions up to the desired order or lower-order perturbed density matrices:

PCM-SCF linear response: comparison with previous formulations
Derivations of the linear 50,51 and nonlinear response functions 52,53 for the PCM-SCF model have previously appeared in the literature. All previous derivations exploit the definition of a solute Hamiltonian which is nonlinearly coupled to the classical dielectric continuum. 13,43 In such a framework, the solvent polarization is not treated as an independent, variational degree of freedom. Solvent contributions to the Hamiltonian are partitioned based on their order dependence on the density matrix: zeroth, first (linear) or second (quadratic) order. We remark that one-and two-electron contributions to the energy are also linearly and quadratically dependent, respectively, on the density matrix. Solvent contributions will thus enter into response theory expressions in much the same way as the proper one-and two-electron terms do.
A derivation of open-ended response theory for an SCF solute coupled with a classical description of the solute has already been presented in the context of the PE MM model. 20 There, the abovementioned order dependence on the density matrix of solvent contributions, which arises when a nonlinear Hamiltonian is invoked, was used to facilitate the identification of the polarization terms to be included in the open-ended formulation of electric response properties. That derivation can also be used when the classical solvent model is implicit, such as the PCM considered in the present work, and will in this case lead to a specific implementation strategy, vide infra. However, the converse is also true. As shown by Lipparini et al., 54 a variational formulation can also be used for classical polarizable explicit solvation models.
To the best of our knowledge, the first derivation of the linear response function exploiting the variational formulation for a quantum/classical polarizable Hamiltonian was presented by Lipparini et al. in ref. 19. Our derivation naturally includes general perturbations, if the fixed-cavity approximation is assumed, and avoids the use of nonlinear Hamiltonians, representing a clear theoretical advantage.
An explicit example: first-order, electric response properties. We here report explicit expressions for the first-order response equations. The differentiated TDSCF condition of eqn (24) evaluated at zero perturbation strength is Decomposing into frequency components and introducing the H/P partition for the variational densities yields: where all the contributions not depending on H-type terms are where F b xc;o contains derivative terms of the XC matrix that are independent of the response parameters. We refer to eqn (A26) of the original paper for its explicit expression. 9 Reorganizing eqn (43) to have all terms dependent on X b on the left-hand side (LHS) yields where we recognize the action of the propagator [E [2] À o b S [2] ] on the response vector X b . Finally, the right-hand side M b letting us cast the first-order response equations in the usual form Since perturbation-independent basis sets are usually employed in the calculation of electric response properties, considerable simplifications arise in all expressions. As an example, we will illustrate how the equations look when only electric dipole perturbations are considered, and we will use the symbol f for such perturbations. First of all, the particular perturbed variational densities are identically zero: one also has s f P = 0. Therefore, only terms including F b o will enter the RHS and among these, only those involving V t,b o will be nonzero, so that [E [2] À o b S [2]

Implementation
The algorithmic realization of the open-ended PCM-SCF response theory presented requires the solution of the coupled response equations for the homogeneous density matrix, eqn (38), and ASC, eqn (36a). Hence, one can envision two possible strategies for a computer implementation of the open-ended scheme: Strategy 1 eqn (36a) and (38) are simultaneously solved, in much the same way as described by Lipparini et al. for the SCF equations. 45 A suitable initial guess is provided for both densities and the appropriate iterative linear equation solvers are employed. 55,56 A convergence acceleration scheme might also be exploited. 57,58 Strategy 2 given the initial guess for the response parameters, the homogenous density matrix is formed and the perturbed MEP calculated. Eqn (36a) is solved, eventually allowing the computation of the linear transformation in eqn (39) and the solution of eqn (38).
The two strategies are of course expected to lead to identical results. Strategy 1 could be advantageous when large molecular solutes are considered and might show better convergence properties, at the expense of a nontrivial extension to both the quasienergy derivative Lagrangian framework for the efficient elimination of response parameters 9,47 and to the corresponding recursive implementation. 10,11 The implementation we present in this work follows Strategy 2. This avoids excessive modifications to the recursive core of the open-ended response code and allows a straightforward use of efficient response parameter elimination schemes. 47 The resulting computer code exploits interfaces between the Dalton program package, 59 the PCMSolver library, 30 the implementation of the open-ended, recursive approach to atomic orbital-based density matrix response theory 10 and the subsequent development of the latter for the calculation of single residues of response functions. 11 The interface to PCMSolver provides an alternative implementation of PCM capabilities at the SCF level of theory to the one presented by Cammi et al. 60 To test the linear transformations of the generalized Hessian and metric matrices of eqn (39) and (40), a non-recursive PCM-SCF linear response code was implemented exploiting the PCMSolver library. 30 Further testing for the linear, quadratic and cubic response functions and corresponding single residues for electric dipole perturbations was achieved by comparing to previously published, non-recursive implementations available in Dalton. [50][51][52][53] The geometries were optimized in vacuo using Gaussian 09, 61 the B3LYP XC functional [62][63][64][65] and the cc-pVQZ basis set. 66 No constraints were used on the symmetry during the geometry optimization, yielding point groups C s , D 2h and C 2v for PNA, PDNB and MCP, respectively. The structures of PDNB (identical to the one used in ref. 12) and MCP are planar, whereas the structure of PNA (identical to the one used in ref. 11) is non-planar with a H-N-C-C dihedral angle of 211. All PCM calculations are done on the vacuum geometry so that all differences arise from direct solvent effects rather than from indirect (geometrical) effects. The ESI † contains the molecular structures used in the calculations, in Dalton input format.

MPA calculations
The calculations were performed in a development version of the Dalton program using the open-ended response code of Ringholm et al. 10 One-electron property integrals and their arbitrary-order derivatives were provided by the Gen1Int library. 67 The XC functionals library XCFun 68 and the integrator library XCint were used for the evaluation of the XC terms and their derivatives. Finally, the PCM functionality was provided by the PCMSolver library. 30 The ESI † contains details about software versions and build toolchain used in this work. Vertical excitation energies were calculated for the 7 lowest excited states of the molecules in Fig. 1. MPA transition moments S were calculated for the same excitations from the residues of the response functions using the implementation described in ref. 11. The CAM-B3LYP XC functional 69 and the aug-cc-pVDZ basis set 66 were chosen based on previous results in ref. 11.
Rotationally averaged transition strengths hd MPA i were calculated from the transition moments S and their complex conjugates % S as 70 The dimensionality, in atomic units, of the MPA strengths is systematically given by 11 where a 0 and E h are the atomic units for length and energy, respectively, and M is the number of photons involved. MPA strengths in atomic units are used throughout unless otherwise stated.

PCM details
The No local field effects were included in this study. [71][72][73][74][75][76] The molecular cavities were generated from the van der Waals surface, i.e. from a set of atom-centered, interlocking spheres. For every molecule this was computed from the Bondi-Mantina set of van der Waals radii: 77,78 1.20 Å for hydrogen, 1.70 Å for carbon, 1.55 Å for nitrogen and 1.52 Å for oxygen. All radii were scaled by a factor of 1.2. Cavity generation and discretization was performed according to the one-point centroid collocation GePol scheme, deactivating the addition of non-atom-centered spheres. 24,25 Since vacuum geometries are used in all calculations, the cavity is the same across all solvents.

Results
The excitation energies of the lowest electronic excitations in PNA, PDNB and MCP are shown in Fig. 2 as a function of the Onsager factor e s À 1 e s , where e s is the static relative permittivity. 79 For some states, the excitation energy increases (e.g. 3A 0 in PNA) or decreases (e.g. 1B 2u in PDNB) considerably with solvent polarity, leading to a change in the ordering of the different excited states.  of the non-equilibrium formulation of PCM that we employ. The calculated solvatochromic shifts are smaller than the experimental ones (see Appendix B), which is likely related to the (deliberate) omission of indirect solvent effects. MPA strengths and dominating orbitals for selected excitations in PNA are shown in Fig. 3 and 4. Among the first five excitations of PNA (thus excluding the 4A 00 state), the 2A 0 excitation is the strongest excitation in vacuo for 1PA up to 5PA, as has been observed for 1PA-4PA before. 11 This state is also the brightest in 1PA and 2PA for all solvents examined. The sixth excitation (4A 00 ), however, is one of the excitations that has a higher strength for 3PA-5PA in the more polar solvents. This illustrates that even qualitative results are poorly transferable from vacuum to solvent and from 1PA to MPA.
The discontinuities in the MPA strengths when increasing the solvent polarity can again be related to the use of the nonequilibrium formulation of PCM. Whereas the static permittivity (e s ) increases monotonously from left to right due to the choice of the x-axis, the dynamic permittivity (e N ) does not. The discontinuities thus reflect the variation in e N with nonmonotonous behaviours observed for aniline e s À 1 e s ¼ 0:85 and dimethylsulfoxide e s À 1 e s ¼ 0:98 . This effect is hardly observable for 1PA, but is pronounced for the higher-order MPA strengths, and for 3PA and 5PA in particular (Fig. 3). MPA strengths can be selectively enhanced by an intermediate state at a defined fraction of the energy of the state in question. 12  This resonance enhancement explains the dramatic increase of the 4PA strength of the 4A 00 state in the vacuum calculation. The excitation energy of the 2A 0 excitation happens to be at 0.749 times the excitation energy of the 4A 00 excitation (Fig. 2), giving rise to a resonance condition and enhancement of the 4PA strength. Resonance enhancement also contributes to the high 3PA strength of the 4A 00 in dimethylsulfoxide (e s = 46.7, e s À 1 e s ¼ 0:98), where the excitation energy of the 2A 0 excitation   is 0.678 times the excitation energy of the 4A 00 excitation. It is important to note that our response theory approach breaks down close to such resonance and the results cannot be considered quantitatively accurate. To circumvent this divergence in the MPA strengths, damped response theory should be used, 80,81 as reported for 2PA by Kristensen et al. 82 The variation of the MPA strengths with solvent polarity clearly increases with the number of photons. This has also been observed experimentally, comparing 1PA and 2PA for the same excitation in PNA. 83 The inclusion of solvent effects thus becomes increasingly important for quantitative MPA calculations with increasing number of photons. Comparison with experimental cross sections 83 indicates that the absolute value of the 2PA cross section is underestimated by up to a factor of 2 in the calculations, while the relative strengths across different solvents are not reproduced in the calculations (Appendix B). The underestimation of the 2PA cross sections may be partially caused by the use of TDDFT and can be related to underestimated difference dipole moments. 84,85 Indeed, TDDFT/CAM-B3LYP has been shown to underestimate difference dipole moments in PNA. 84,86,87 We now turn our attention to the centrosymmetric molecule PDNB. MPA strengths and dominating orbitals for selected excitations in PDNB are shown in Fig. 5 and 6. The high point group symmetry determines whether the different excited states will be allowed for the different one-and multiphoton absorption processes. The 1B 3u excitation is the brightest excitation for 1PA, 3PA and 5PA among the first seven excitations, while the B 2g excitation is the brightest for 2PA and 4PA. We note that the intensities of the different odd-or even-order photon excitation processes follow each other both in vacuum and in the different solvents. This is due to the fact that the intensity of oddand even-order multiphoton absorption processes can be shown to be proportional to the one-or two-photon absorption cross section, respectively, 12 with exceptions arising due to near-resonances.
As a final example, we have considered MCP: some of its excited states undergo a sign change in dipole moment across the solvent polarity scale chosen in this study. The MPA strengths and dominating orbitals for selected excitations in MCP are shown in Fig. 7 and 8. The pp* 2A 1 transition has significant orbital overlap between the occupied and the virtual orbitals (Fig. 8) and is the only strong transition in 1PA. The 1B 2 and 2B 2 states, however, are accessible via multi-step processes and are indeed the strongest states in 2PA/4PA and 3PA/5PA, respectively. The maximum in the 4PA strength for the 2A 1 excitation is due to a resonance with the 1B 1 state which is located at    related to the excited-state dipole moment, which changes sign between 1,4-dioxane and toluene. The excited-state dipole moment of the 1B 2 state also has a different direction in vacuum than in solvent. We report the relevant plots of the difference between ground-and excited-state dipole moments in the ESI. † The results presented serve as an illustration of the applicability of our implementation. When interpreting the results, one should bear in mind the limitations of our methodology. First of all, we have calculated vertical excitation energies and neglected vibronic effects, which have been demonstrated to play a role in 2PA. 88,89 Secondly, we have not taken into account the indirect effect of the solvent on the geometry of the chromophore. Indirect solvent effects can be taken into account by using PCM also in the geometry optimization, which has however not been done here, to allow for a comparison of the direct contribution of the various solvents on the MPA strength. Thirdly, the implementation of the solvent model used here does not include explicit local field effects in the molecular cavity 71-76 nor non-electrostatic effects. 31,32,90,91 Explicit solute-solvent interactions are also not included, but can in principle be recovered using a QM cluster model of the system. Finally, DFT is not likely to give MPA strengths (and, in general, excited-state properties) of high accuracy, as shown for 2PA using a coupled-cluster benchmark. 85 There is a clear need for benchmarking DFT MPA strengths against methods of higher accuracy also beyond 2PA. 85,92 All these factors are important if a realistic comparison with experiment is to be attempted. The recent coupling of our open-ended response code to a PE QM/MM framework is able to take into account indirect solvent effects, local field effects and explicit solutesolvent interactions, albeit by necessity introducing configurational sampling in the computational protocol. 20

Conclusion
We have presented the theory and implementation for calculating molecular response properties to arbitrary order in solution within the framework of the polarizable continuum model. The theoretical derivation is based on an energy functional where both the density matrix and the electrostatic polarization in the medium are treated as variational degrees of freedom. Contrary to previous work, the quantum/classical polarizable coupling is not achieved by assuming a nonlinear interaction potential in the Hamiltonian. We have shown that, in the fixed-cavity approximation, molecular response functions to arbitrary order are straightforwardly obtained as higher-order derivatives of the proposed functional. Moreover, differentiation of the stationarity conditions naturally leads to the appropriate response equations determining higher-order perturbed wave function and polarization parameters. Our implementation relies on modular components encapsulating the different tasks required to carry out a response calculation, in line with previous work by some of us. 10,11 In particular, we added the PCM terms to the workflow by means of the PCMSolver library, 30 in the spirit of the PE implementation recently presented by Steindal et al. 20 We have illustrated the implementation by calculating MPA strengths for three small organic molecules. The enhancement of the MPA strength from vacuum to different solvents increases with the number of photons involved in the excitation, clearly emphasizing the importance of including solvent effects in MPA calculations. Relative intensities between features corresponding to different electronic excitations in one-or multiphoton absorption spectra are not necessarily preserved between phenomena involving different numbers of photons absorbed, which is partially related to molecular symmetry. We have also described resonance enhancements in our MPA calculations.

B Comparison with experimental data
Although the paper is essentially methodological, we provide a short comparison of our results with available experimental data. A detailed study of the solvatochromic shift of PNA including 2PA spectra has been presented recently by Wielgus et al. 83 For PDNB, maximum absorption has been observed at 4.75 eV in ethanol. 93 For MCP, two values have been reported in the literature: 4.01 eV in n-pentane and 4.49 eV in methanol. 94 In Table 1, we report our vertical excitation energies (VEEs) together with available experimental results, when the same solvent was available (our calculated values in n-heptane are however compared to experimental data obtained in n-hexane for PNA and in n-pentane for MCP).
For PNA, our results are overall qualitatively well correlated with the experimental absorption maxima. However, the calculations are more compressed (an overall variation of 0.25 eV compared to 0.68 eV for the experiments) and water and 1,4dioxane are slight outliers. A possible explanation for the range compression is the lack of indirect solvent effects (vacuum geometries have been employed throughout, a conscious choice to isolate the direct effects on the MPA cross sections), whereas specific solvent effects might explain why water is an outlier.
For MCP, the computed solvent shift (0.22 eV) is also less pronounced than for the experimental value (0.48 eV) but in good agreement with previous state-specific multiconfigurational calculations (0.18 eV) reported in ref. 60. For PDNB, we cannot draw any conclusions from one isolated number and the observed agreement with experiment is likely to be fortuitous.
2PA strengths for PNA in three solvents were converted to cross sections for comparison with available experimental data 83 using 85 where a is the fine structure constant, a 0 is the bohr radius, c is the velocity of light, o is the photon energy (i.e. half the excitation energy) and G is the half-width half-maximum (HWHM) of the (Lorentzian) broadening function describing homogeneous broadening. The HWHM for 1,4-dioxane, propylene carbonate and dimethyl sulfoxide were taken from reported experimental work (FWHW 4672 cm À1 , 4490 cm À1 , 4014 cm À1 ), 83 giving G = 0.2897 eV, 0.2783 eV and 0.2488 eV, respectively. The dimensionality of the 2PA cross section is cm 4 s photon À1 with 1 Â 10 À50 cm 4 s photon À1 commonly referred to as 1 GM (Table 2).