Nonequilibrium Solvent Polarization Effects in Real-Time Electronic Dynamics of Solute Molecules Subject to Time-Dependent Electric Fields: A New Feature of the Polarizable Continuum Model.

We develop an extension of the time-dependent equation-of-motion formulation of the polarizable continuum model (EOM-TDPCM) to introduce nonequilibrium cavity field effects in quantum mechanical calculations of solvated molecules subject to time-dependent electric fields. This method has been implemented in Octopus, a state-of-the-art code for real-space, real-time time-dependent density functional theory (RT-TDDFT) calculations. To show the potential of our methodology, we perform EOM-TDPCM/RT-TDDFT calculations of trans-azobenzene in water and in other model solvents with shorter relaxation times. Our results for the optical absorption spectrum of trans-azobenzene show (i) that cavity field effects have a clear impact in the overall spectral shape and (ii) that an accurate description of the solute shape (as the one provided within PCM) is key to correctly account for cavity field effects.

Introduction

The electronic and optical properties of a molecule embedded in a dielectric environment, such as a solvent, can be substantially different from those of its gas-phase counterpart. The presence of the molecule induces a polarization in the surrounding that in turn interacts with the molecule itself via the so-called reaction field. Additionally, if the molecule is subject to an external electric field, the effective field acting locally on the molecule (namely, the cavity field) contains another polarization contribution of the environment due to the applied field.[1] Hence, when simulating a solvated molecule interacting with any electromagnetic perturbation (e.g., light), both the reaction and cavity field effects must be considered.

A suitable framework to account for polarization effects in quantum-mechanical calculations of molecules surrounded by a dielectric environment is the polarizable continuum model (PCM).[2−5] Within PCM, (i) the environment is regarded as a continuum and infinite dielectric medium featuring a frequency-dependent dielectric function; (ii) the molecule is treated explicitly, that is, by considering its atomic structure, and generally within quantum mechanics; (iii) a void cavity hosts the molecule, conforming with its shape and separating it from the dielectric environment by a sharp interface; and (iv) the reaction- and cavity-field polarization of the medium are described in terms of apparent surface charge (ASC) densities sitting on the cavity surface that are induced by the molecule and the applied electric fields, respectively.

The PCM formulation can accommodate nonequilibrium effects in the molecule-environment dynamics when shifting from the electrostatic[2,5] to the time-dependent quasi-electrostatic[3,6−9] regime. However, until recently, nonequilibrium features were only accessible through the sudden excitation approximation, where the molecular density changes abruptly from the ground to the excited state. In response to this sudden excitation, the solvent degrees of freedom are crudely divided in two groups: the “slow” ones that are fully unable to adapt to the molecular density rearrangements, and the “fast” ones that follow instantaneously the molecular density evolution. Although, this slow/fast splitting of the environment polarization contributions represents a qualitatively correct scheme to treat nonequilibrium polarization effects, it neglects the fact that even the slow environment polarization (albeit dominated by its inertial character) has a finite relaxation time and can as well adjust to the molecular excitation in the long term.[2]

Currently, a computationally convenient time-dependent extension of PCM (TDPCM) has made possible to study the real-time evolution of the molecular density subject to time-dependent quasi-electrostatic perturbations, together with the polarization of the external medium.[10−13] Such a TDPCM scheme introduces polarization fields that are history-dependent and not instantaneously in equilibrium with the molecule or the applied fields. These nonequilibrium polarization effects arise ultimately from the frequency-dependent dielectric function describing the environment, and they are physically linked to the fact that the environment takes a finite amount of time to react in response to the electrostatic perturbations, that is, the molecular density fluctuations and the applied fields. Mathematically, this formulation can be cast in terms of a set of equations of motion (EOM) for the ASCs dictating the evolution of the medium polarization in real-time. In the case of solvated molecules, EOM-TDPCM has been formulated for the reaction field;[13] however, a similar account for the cavity field is still missing.

On a different note, even though treating the reaction field within the quasi-electrostatic regime is generally pertinent, an equivalent description of the cavity field does not apply, in principle, when a molecule surrounded by a dielectric medium is irradiated with light. To overcome this limitation, in ref (14), cavity field effects have been rigorously tackled within classical electrodynamics, by using the full set of Maxwell equations to describe the cavity modification of electromagnetic fields traveling in a medium when impinging a vacuum cavity. Notwithstanding, under convenient gauges and the relevant long-wavelength limit,[14,29] the electrodynamics result has been brought to match the standard PCM account of the cavity field effects,[3−5] where an effective dipole of the embedded molecule interacts with the time-dependent electric field of the light wave.

The aim of this paper is to present an extension of EOM-TDPCM to consider nonequilibrium cavity field polarization effects for molecules embedded in an homogeneous dielectric medium, such as a solvent. This contribution represents a step further in considering more realistic solvent polarization effects. Nonequilibrium cavity field effects may impact, strongly and nontrivially, the electronic and optical properties obtained from real-time quantum mechanical simulations of solvated molecules, especially when at least some of the solvent degrees of freedom relax in a short time scale (0.1–1 of ps). Cavity field effects treated within the PCM framework enable us to consider in a more realistic manner complex molecules, where the solute shape departs strongly from symmetric shapes (such as a sphere or an ellipsoid) and furnishes quantitative an qualitative improvement over some standard descriptions for the cavity field, such as Onsager’s model.

To show the potential of our methodology, we perform real-space real-time time-dependent density functional theory (TDDFT) simulations of solvated trans-azobenzene to assess the influence of the cavity field in the optical spectrum and time-dependent properties. The molecule of choice does not have spherical symmetry and is largely anisotropic. As a consequence, some of its relevant excited states feature transition dipoles oriented along different directions. trans-Azobenzene is, therefore, a good candidate to emphasize the nontrivial role of cavity field effects in changing the optical spectrum shape. We present results for the molecule immersed in water. Furthermore, to highlight the role of nonequilibrium cavity field effects, we provide comparison between the results obtained for the molecule solvated in water and in different fictitious solvents with faster orientational relaxation.

The present account is organized as follows. In section , we recall the PCM theory of cavity field effects for a molecule embedded in a dielectric environment. Section is devoted to the electrostatics formulation of cavity field effects. In section , we extend the latter to tackle the case of time-dependent electric fields in the long-wavelength limit and using the length gauge for the case of a molecule in solution. Section introduces time-dependent quantum mechanics for electrons in a solvated molecule subject to an external electric field. Section focuses on a specific time-dependent quantum-mechanics treatment, namely, real-time TDDFT (RT-TDDFT). In section , we obtain the expression for the full photoabsorption cross-section spectrum for solvated molecules, emphasizing cavity-field contributions and the ways to account for them through real-time propagation of an effective molecular dipole. Sections and 2.5 are specific for the case of solvent, while the rest are rather general for any medium considered. In section , we present the computational details (section ) and the results for the optical absorption spectra (section ) of trans-azobenzene in water. Section focuses on the cavity field and nonequilibrium cavity field effects in the optical absorption spectrum of trans-azobenzene. Finally, in section , we summarize the relevant points of the Article, and we stress the potential of our methodology for future and more accurate studies of molecules in solution.

PCM Theory of Cavity Field Effects

Electrostatics

Let us consider a static electric field applied to the solution. Far from the vacuum cavity hosting the solute, the dielectric representing the solvent is homogeneous and the electric field there is the Maxwell field, FM.[16] Close to and within the solute cavity, the electric field is locally modified by the presence of the cavity itself. Any distribution of charges placed inside the cavity (e.g., a molecule) is affected not only by the Maxwell field as if it were applied in an homogeneous medium but also by the medium polarization because of the cavity boundary. The effects of such electric field contribution supplementing the Maxwell field are called cavity field effects. To treat them within a PCM framework, one defines an effective electrostatic potential.[14]

Such effective applied potential (that includes cavity field effects) acting on any distribution of charges, vC(r) to keep the nomenclature of ref (14), is the solution of a generalized Laplace problem, where the dielectric constant is position-dependent, with boundary conditions at the cavity surface and far from it (see refs (1, 2, and 5)). The latter generalized Laplace problem can be equivalently written as a generalized Poisson problem, considering a fictitious distribution of charges ρM(r) placed outside the vacuum cavity (however far from it) as the source of the Maxwell field. The PCM framework,[2] in particular, in its integral equation formalism (IEF) approach[17] allow us to solve this kind of Poisson problem. Ref (17) shows a detailed IEF-PCM formulation in the reaction field case that we followed thoroughly to obtain an analogous account for the cavity field case.

As it is customary within PCM framework, we split the solution of the aforementioned Poisson problem into two contributions, namely,[18]However, the specific partition of the potential we employ differs from the usual choice in the case of the reaction field.[2,17] Here, the first contribution, vM(r), is the potential associated with the Maxwell field, whereas the remaining term, vcf(r), represents a polarization potential that solves the Laplace equation inside and outside the cavity and fulfills the boundary conditions from the original Poisson problem. Thence, vcf(r) is a so-called layer potential that can be written in terms of ASC density σcf(r) located at the interface Γ, that iswhere dA(s) is the surface area element corresponding to s ∈ Γ.[2,17]

Within the numerical strategy known as boundary element method (BEM), we consider a discrete (but fine) tessellation of the cavity surface and approximate the latter as[2]where Ntess is the number of elements (tesserae) in Γ’s tessellation, whereas s, A, and qcf are the position, the surface area and the polarization point charge associated with the ith tessera.

IEF-PCM equations for the polarization charges {qcf} in terms of the potential that induces them can be written as[2,17] or in matrix form, . , with j, k = 1, ..., Ntess, is the PCM response matrix depending on the cavity shape, the tessellation of its surface, and the dielectric constant of the medium.[2,17]

To get the cavity field PCM response matrix starting from our partition , we follow an analogous derivation to that laid down in ref (17) to obtain the reaction-field PCM response matrix. It yieldswherewith j, k = 1, ..., Ntess, and n(s) is a unitary normal vector to the k-th tessera, pointing outward from the cavity. Matrices and are the BEM representation of Calderon operators in the IEF.[17] Heretofore, we have introduced very general definitions and, in particular, we have not assumed that the dielectric constant for the external medium is position-independent, as it is for a typical solvent. Indeed any in-homogeneous external dielectric medium, with a position-dependent dielectric constant of its own, can be tackled in the same way provided that the and matrices are known, which is to say that the Green’s function of the external generalized Poisson problem were found.[2,17] Particularizing to the solvent case, the PCM response matrix takes the form

The cavity field PCM response matrices and 7 differ from their reaction-field counterparts by two sign changes: one to the full expression and another one affecting the term inside the second (noninverted) bracket; both “–” for reaction-field and “+” for the cavity-field case (cf., eq in ref (13)). PCM response matrices for the reaction- and cavity-field cases are different because the source charges inducing the medium polarization are located in different regions of the space in each case, that is, inside the vacuum cavity for the reaction field and outside (in the external medium) for the cavity field.

Partition is not the only starting point to solve the original generalized Laplace problem for the applied electrostatic potential within the PCM framework. In Appendix , we develop, for a different partition, an alternative derivation to cavity field effects within PCM.

Electrodynamics in the Long-Wavelength Limit

In the previous section, we considered the static polarization response of a dielectric medium with a vacuum cavity due to an externally applied electrostatic field. Here, we study the dynamic polarization response of the medium because of a time-dependent electromagnetic field. This can be, for example, the electromagnetic field associated with a monochromatic light pulse probing the molecule in a specific spectroscopy.

Far from the cavity, in the solvent bulk, the electric field is again that proper for an homogeneous medium and is defined as the Maxwell field. Close to the molecule, we work in the long-wavelength limit and under the length-gauge,[14] so that locally the scalar Maxwell potential can be written asin the frequency domain. Within these assumptions, we may solve the generalized Poisson problem in the frequency domain, considering the complex dielectric function dependence on the frequency, that is, ϵ(r;ω) . We solve the latter by means of the PCM eq and the polarization potential (eq ).

Let us focus on an external medium with the characteristics of a solvent (i.e., position-independent dielectric constant). In such a case, the PCM response matrix is known in terms of the frequency-dependent dielectric function of the external medium (eq with ϵ → ϵ(ω)). Hence, the PCM eq holds in the frequency domain as .[2,13] To obtain the time-dependence of the polarization charges, we have to Fourier transform the latter PCM equation, that is,[13,19]The former expression is central in this work and reflects the physical meaning of the nonequilibrium cavity field effects. It shows that cavity field polarization charges—and therefore, the cavity field—depend upon the full previous history of the Maxwell potential and not just on its values at the current time.

Debye EOM for the Cavity–Field Polarization Charges

In the case of a solvent described by Debye’s dielectric model, the evolution of cavity field polarization charges—governed by the latter expressions—can be reformulated in terms of an EOM. Such EOM was derived by Corni, Pipolo, and Cammi in ref (13) for the reaction-field polarization charges. It is also valid for the cavity field polarization charges once the molecular potential is replaced by the Maxwell electric field and the proper propagation matrices are obtained from the new form of the PCM response matrix, see eq .

We recall that Debye’s model represents the dielectric medium by a set of freely rotating point-like dipoles whose orientation is determined by random collisions and by the applied electric field.[20] Therefore, Debye’s model captures the polarization response of the medium to an electric perturbation due to the orientational relaxation of the solvent molecules. Within Debye’s model, the external medium is characterized by the complex dielectric functionwhere τD is the so-called Debye relaxation time, which is a physical measure of how fast the dielectric environment can adapt to the changes of the potential it is subject to. Mathematically, it interpolates between the static and dynamic dielectric constants in the limit of small and large relaxation time, ϵ0 and ϵd, respectively.

The resulting EOM for the cavity field polarization charges in the Debye case iswhere Λ and τ are diagonal matrices. The only changes in the EOM propagation matrices for the cavity field with respect to reaction-field case (cf., ref (13)) are located in the definition of in eq ; more explicitly, they are sign changes (i) to the full expression (“+” for “–”), and (ii) to the matrix Λ in the noninverted bracket (“–” for “+”). For a propagation with a finite Maxwell field switched-on at the initial time t = 0, q(0) = 0. If we consider a medium with an infinite relaxation time, we getThe case of a Dirac delta time-dependence is special and it deserves a separated treatment (see Appendix ).

Finally, the time-dependent equivalent of eq is

Time-Dependent Quantum Mechanics

The time-dependent effective electronic Hamiltonian of a molecule embedded in a dielectric environment and interacting with a probing electric field in the long-wavelength limit, using the length gauge and within the PCM framework is[14] is the Hamiltonian of the embedded molecule in absence of external fields, collecting the isolated molecule (in vacuo) Hamiltonian and the interaction with the polarization of the surrounding medium due to the molecule itself (reaction field). On the other hand, V̂C(t) = V̂M(t) + V̂cf(t) describes the interaction with the total electrostatic field acting on the molecule, including the interaction with the bare external field in the dielectric medium (Maxwell field) and with the polarization of the medium due to the applied electrostatic field (cavity field term). As long as the nuclei are fixed, the isolated molecule Hamiltonian Ĥmol is time-independent. The time-dependent interaction Hamiltonian can be written aswhere Ne and {r} are the number of electrons of the molecule and their position operators, respectively. Notice that the time-dependent reaction-field potential, vrf(r,t), depends on the time-dependent molecular potentialat any earlier time t′ < t.[2,13] These time-dependencies can be taken into account altogether by means of an EOM for polarization charges equivalent to that in eq that only depends on the molecular potential at current times.[13] In eq , ρmol(r,t) is the time-dependent molecular charge density, collecting electronic and nuclei contributions.[2]

RT-TDDFT

In this Article, the quantum-mechanical description of the molecule embedded in a dielectric medium and interacting with an external electric field (including reaction and cavity field effects) is tackled at the level of DFT[21] and TDDFT.[22] Particularly, we use the Kohn–Sham (KS) approach within DFT,[23] justified by the Hohenberg–Kohn theorems,[24] and its time-dependent extension within TDDFT, justified by the Runge-Gross[25] and the van Leeuwen[26] theorems.

We start from the time-independent KS equations[21]where{ψ(r)} and {ε} are the KS Hamiltonian, molecular orbitals and energies, respectively, whereas ρ0(r) is the ground-state electronic charge density. ĥmol[ρ0] is the KS Hamiltonian for the isolated molecule (in vacuo and in the absence of external fields) including the kinetic energy operator, nuclear attraction, Hartree, exchange and correlation potentials.

The time-dependent KS equations,[22]whereand {ψ(r,t)} are the time-dependent KS Hamiltonian and molecular orbitals, respectively, whereas ρ(r,t) is the time-dependent electronic charge density. The adiabatic approximation for the exchange and correlation term of KS Hamiltonian is usually assumed, that is to say that ĥmol(t) is a functional of the density at the current time, that is, ĥmol[ρ(r,t)].[22] However, the reaction-field term, vrf[ρ] (r,t), is instead intrinsically diabatic and it depends on the density at earlier times.[2]

Provided that the solution for the starting time t = 0 is known, the propagation scheme[27] solves eq . Within RT-TDDFT schemes using the adiabatic approximation, the initial KS orbitals are obtained from the time-independent KS equation (see eq ), that is, ψ(r, 0) = ψ(r).

Optical Absorption Spectrum of a Solvated Molecule

In this section, we find the photoabsorption cross-section of a solvated molecule, aiming at including and singling-out cavity field effects. In addition, we assume an electromagnetic plane wave impinging on the system as the source of photoabsorption processes. Furthermore, we assume to be in the long-wavelength limit (within the length-gauge) and in the linear regime.[28]

The intensity of a monochromatic electromagnetic radiation with frequency ω traveling in a nonmagnetic medium is[29]where FM(ω) is the amplitude of the Maxwell field and n(ω) is the refraction index of the medium. We recall that . Within our framework of assumptions, the rate of energy dissipation of a monochromatic electromagnetic radiation applied to the system can be written as[29]where ρmol(r,ω) is the charge density of the molecule induced by the total potential vC(r,ω), both in the frequency domain (i.e., Fourier transformed from ρmol(r,t) and vC(r,t), respectively), collecting the Maxwell and cavity field terms. Finally, the photoabsorption cross-section is the ratio between eqs and 22, yieldingUsing partition and vM(r,ω) = −FM(ω)·r in eq , we getwhere dmolcf(ω) is the dipole of the molecule induced by the total potential in the frequency domain. By incorporating a superscript “cf”, we emphasize that the latter includes implicitly cavity field effects, but it is not the only source of them in the cross-section spectrum. We now focus on the term Σcf(ω), which is the contribution to the cross-section coming from the polarization of the medium because of the Maxwell field and, henceforth, accounts explicitly for further cavity field effects. Using the expression of the potential vcf(r,ω) in terms of its sources σcf(r,ω), we can writewhere vmol(s,ω) is the Fourier transform of the quantity defined in eq . The polarization charge density due to the Maxwell field, σcf(s,ω), is a linear function of the Maxwell field, FM; therefore, by applying Euler theorem for linear functions, we get[14]where

Using eq , we may write the complete cross-section aswhere d̅mol(ω) = dmolcf(ω) + d̃(ω) is the effective induced dipole of the solvated molecule (external dipole in Onsager’s nomenclature[1]). We recall that the cross-section of a molecule without cavity field effects can be obtained from eq replacing vC(r,ω) by vM(r,ω), yieldingwhere dmol(ω) is the molecular dipole induced by the Maxwell field.

Notice that we can write the induced dipoles aswhere is the frequency-dependent first-order polarizability tensor. The first superscript labels the observable under examination and the second superscript keeps track of the magnitude coupled with the perturbator.[22,30] In eqs and 31b, we monitor different observables, namely, the induced molecular and effective dipole (second superscript d and d̅, respectively), considering different interaction Hamiltonians as a perturbation, that is, V̂M(t) = −d̂mol·FM(t) and V̂C(t) = −d̂®mol·FM(t) (first superscript d and d̅), respectively. In eq , we fully consider cavity field effects, whereas in eq , we do not.

From eqs and 31b, we get for the cavity field casePerforming an orientational average over all the possible field orientations, we getAnalogously, for the case without cavity field effects, we have

The key quantities in the optical spectra of eqs and 34 are the polarizability tensors and with and without cavity field effects, respectively. The polarizability tensor components may be written asParticularly, we are interested only in the diagonal terms, since such are the components relevant for a photoabsorption cross section spectrum obtained from unpolarized light (see in eqs and 34). We stress that d(ω) and d̅(ω) are induced dipole moments, not permanent ones. We can write them using the Fourier transform from the time-domain asTo obtain the optical spectrum, an impulsive electric field independent of the position and with a Dirac delta time-dependence, that is, FM(r, t) = FM0δ(t) is often employed as a perturbator of the system. Being frequency-independent, the latter field is suitable to perturb equally all the dipole-allowed excitations of the system. Using this kind of external (Maxwell) electric field, we get the full expression for the polarizability tensor components as

The time-dependent effective and bare molecular dipole moments are written aswhere dmolcf(t) is also computed from eq in the case of a cavity-field time-dependent simulation, and

The time-dependent convolution in eq is unsuitable for an implementation due its high computational cost. Therefore, to compute this part of the cavity field contribution to the optical absorption spectrum we follow a different approach stemming from an alternative partition for the total electrostatic potential (see Appendix , cf., eq ), and described in the Appendix .

Results

Computational Details

To illustrate the operation of our methodology accounting for nonequilibrium cavity field effects within TDPCM, we perform RT-TDDFT simulations and compute the optical absorption spectrum of a case-study solvated system. We choose azobenzene in its trans configuration as a solute molecule. trans-Azobenzene geometry is planar, and it was optimized in vacuo at the DFT level. We work under the adiabatic local density approximation (LDA), and we use standard pseudopotentials consistent with the LDA exchange-correlation functional. The solvent is described using a Debye dielectric model with parameters τD = 3.37 ps, ϵ0 = 78.39, and ϵd = 1.776, mimicking water.[31,32] We perform as well linear response TDDFT (LR-TDDFT) calculations of trans-azobenzene in vacuo to obtain the first-excited state transition dipole moment (TDM). To that aim, all the single-particle transitions in an energy window of ∼10 eV about the Fermi energy are considered. Molecule-solvent interactions via the reaction field are always taken into account, regardless of the inclusion of cavity field effects or not. To carry out all computations, we use a modified version of Octopus (version 8.0)[33] (see Appendix for details on the implementation).

In our calculation, the real-space grid is built from interlocking spheres of 5 Å of radius centered at the atomic positions, with a spacing between grid points of 0.18 Å. The real-time simulations were performed for 15 fs with a time step of 1.5 as. We start the propagation of the KS orbitals from the ground state after applying a weak and impulsive external electric field with an amplitude of ∼5 × 10–3 a.u.

Optical Absorption Spectrum with Cavity Field Effects

In this section, we present the results for the optical absorption spectrum of trans-azobenzene in vacuo and in aqueous solution, focusing on cavity field effects.

As anticipated, we obtain the optical absorption spectrum of the solvated molecule by applying an impulsive electric field oriented along each of the three Cartesian directions. The averaged photoabsorption cross-section is calculated using eqs and 34, respectively, if we consider cavity field effects or not. and are computed from the Fourier transform of the time-dependent effective or bare dipole moments (eq ). Finally, we compute the time-dependent dipoles via the time-dependent density and reaction-field polarization charges (eqs , 53, and 58) obtained from RT-TDDFT calculations coupled with EOM-PCM, as described in ref (13) and section for the reaction and cavity field solvent polarization, respectively.

In Figure , we show the results for the optical absorption spectra (strength function ∝ cross-section) for the molecule in vacuo and in solution, with and without cavity field effects. The first two panels contain the absorption spectra using an electric field in the plane of the molecule (see Figure ). In the last panel, we plot the mean value of the absorption spectra in the three Cartesian directions.

Figure 1

Absorption spectra of trans-azobenzene in vacuum (black lines) and water. The case with (without) cavity field effects is represented with red (blue) lines. For the first two panels, the impulsive electric field is parallel to the plane of the molecule (see Figure ). The bottom panel shows the absorption spectrum averaged over the three Cartesian directions of the field (eqs and 34).

Figure 2

Optimized structure of trans-azobenzene in vacuo, along with TDMs associated with A and B bands from Table .

Table 1

TDMs Associated with A and B Bands in Figure a

	d0A (D)	d0B1 (D)	d0B2 (D)
x	5.02	5.42	–2.46
y	–6.07	–2.46	–5.42
z	0.00	0.00	0.00

Since the B band is degenerate, results are shown for two excited states, B1 and B2, contributing the most to B absorption.

In addition to the RT-TDDFT computations, we perform a LR-TDDFT calculation in vacuo to obtain the TDMs corresponding to the selected transitions A an B indicated in the lower panel of Figure . The TDM associated with the first optically active excited state (namely, transition A), can be directly computed from LR-TDDFT results, using the appropriate transition density, as d0A = ∫ρ0A(r)r dr. (We drop the subscript “mol” for the sake of simplicity.) However, LR-TDDFT and RT-TDDFT spectra differ substantially in the high-frequency region, making difficult to match the absorption peaks and hence the underlying excitations. Therefore, to get the TDM associated with the high-energy transition B, d0B, we follow a different path. To find the magnitude of d0B, we use the fact that ratio between absorption peak strengths iswhere ΔE0A and ΔE0B are the excitation energies corresponding to excitations A and B, respectively, and α indicate Cartesian components x, y, z. We extract ΔE0A and ΔE0B from the average absorption spectrum (see Figure bottom panel), whereas fAB,α is taken from the absorption spectra with the field along α direction (see Figure top and middle panel for α = x, y); both quantities are obtained from RT-TDDFT calculations. Using, in addition, the LR-TDDFT result for d0A, we obtain |d0B|. To get the direction of d0B, that is, d0B/|d0B|, we use RT-TDDFT results directly. We extract the optical absorption cross section matrix Σ at the energy of the transition B, and we diagonalize it. For a nondegenerate band, we expect just one eigenvalue different from zero, and since ,[22,34] the normalized eigenvector associated with the nonzero eigenvalue would set the direction for the TDM, d0B. If in turn the band is almost degenerate, we can still find the TDMs with the largest contributions to the absorption by diagonalizing the cross-section matrix and taking the eigenvectors associated with the highest eigenvalues. In our case, we find at least two possible directions for the B band TDMs that we will call B1 and B2 transitions. Results for A, B1, and B2 TDMs are reported in Table . Figure shows a schematics of trans-azobenzene molecule along with the TDMs corresponding to A, B1, and B2 excited states.

Returning to Figure , first, we observe a red shift in the absorption peaks of the solvated molecule with respect to its gas-phase counterparts. The latter is an expected reaction field effect. Cavity–field effects might induce shifts in the spectrum when the solvent dielectric function dispersion in the absorption frequency region is large. In this case, the red and blue ends of an absorption band are affected differently by the cavity field effect, and the band would be deformed. Such a strong dispersion is however not the case for common solvent, and we do not expect shifts associated with cavity field effects.

Second, we check that the cavity field effects change not only the absolute but also the relative strength of absorption peaks, therefore impacting in a nontrivial way on the absorption spectrum. To assess cavity field effects for trans-azobenzene in solution in a quantitative way, we compare our results for the cavity field factor (CFF) with the prediction by Onsager model, which considers the molecule as hosted by a spherical cavity.[1] We select two representative optically active transitions (called A and B in Figure ) and we compute the CFF directly from the spectra for transitions A and B (we refer to it as “PCM” CFF or just CFF), usingand through the Onsager analytic formula, that is,The results are reported in Table . Notice that not only there is a substantial quantitative departure from Onsager prediction for both A and B transitions, but also there is a strong dependence of the CFF on the specific transition that is overlooked by Onsager model.

Table 2

CFF Computed from Eq Compared with the Prediction by Onsager Model, Eq , for Two Selected Transitions from Figure

transition	wavelength (nm)	CFF	Onsager CFF
A	369	1.03	1.37
B	174	1.24	1.37

Since the Debye dielectric function for water is almost constant in the ultraviolet and visible (UV–vis) spectrum, neither the PCM (eq ) nor Onsager CFF (eq ) are able to bear any frequency dependence through the dielectric constant. The reason for the CFF difference for transitions A and B is instead explored in the next section.

Anistropic Cavity Field Effects

Keeping the dielectric function as a constant, the ultimate cause for the disagreement between PCM and Onsager model comes from the fact that the solvent–molecule interface cannot be adequately modeled by a spherical surface for the molecule in question (see Figure ). In the Onsager (spherical) model, the electric field acting on the molecule (cavity field) is increased by the same amount independently on the polarization direction and the frequency. Therefore, cavity field effects treated in such fashion, provide equal enhancement of the strength of every peak in the average absorption spectrum (cf. lower panel in Figure ). If we consider more realistic shapes for the molecular cavity, the absorption strength corresponding to certain excitations would be significantly more enhanced than others depending on the orientation of their TDMs. Even in the case of transitions that are associated with TDMs that are practically colinear this effect may be substantial.

To illustrate the latter, let us consider a still simple but more physically informed shape for the molecular cavity. We remark that trans-azobenzene is more than 2.5 larger in the direction connecting the two benzene moieties than in the orthogonal direction. We notice as well that the transitions of interest are characterized by in-plane dipole moments (see Figure ). To describe the cavity field influence in these two transitions, we are merely interested in addressing the possible effect of the cavity anisotropy in the plane. Hence, we need a cavity featuring different characteristic length in the “long” and “short” axes of the planar structure, disregarding the chosen (perhaps arbitrary) characteristic length in the third direction orthogonal to the molecular plane. We take advantage of the known analytic expression for the CFF in the case a prolate ellipsoidal cavity.[14,35,36] We choose the prolate ellipsoid with the major axis along the direction of the TDM corresponding to A excited state, that is, d0A, and having a circumference section with a diameter 0.38 times as long as the major axis of the in-plane ellipse. Then, we calculate the CFF as the squared ratio between TDMs with and without cavity field effects.[14] The computation of TDM with cavity field effects within the prolate ellipsoidal model amplifies more the components in the shorter axis.[14,36] This setup is sufficient to explain the transition-dependence of the CFF, yielding a CFF = 1.12 for A, 1.20 and 1.46 for B1 and B2, respectively. From the last numbers, we recover the correct trend of cavity field effects, that is, the largest CFF is obtained for the excited states with largest component of the TDM along the shorter axis of the molecule, namely, B2 transition in our case. Within the ellipsoidal model, we have an overall CFF = 1.33 for the B band (assuming B1 and B2 TDMs to be exactly orthogonal) which match qualitatively the PCM data reported in Table . We do not expect to find a quantitative agreement through the former analysis since the molecular cavity is in any case more complex than an ellipsoid and better described within PCM.

Nonequilibrium Cavity Field Effects

Finally, we comment on the nonequilibrium effects as treated using EOM-TDPCM. In section , we obtained the expression for the polarization charges propagation in the case of infinite relaxation time, that is, eq . Using this limit in both the reaction-field and cavity-field time-dependent polarization charges, we obtained the same absorption spectrum as for the EOM-TDPCM plotted in Figure (not shown). The reason is that the Debye relaxation of water is much longer (τD = 3.37 ps) than the characteristic time for electronic excitations (tenth of fs). The slow degrees of freedom of the solvent (driving the inertial polarization) cannot follow the rapid change of either the electronic density or the Maxwell field in the frequency range of electronic excitations (the UV–vis spectrum). Therefore, the dynamics of the solvent polarization relaxation after the application of the impulsive field is dominated only by its fast degrees of freedom, which are characterized by the high-frequency dielectric constant of the solvent, ϵd (see eq ). As a consequence, the only transition-dependence appearing in the CFF is coming from the anisotropic cavity, and not from the dielectric function frequency dependence, which is the ultimate responsible for nonequilibrium solvent polarization effects. Nevertheless, this situation might be drastically different when the relaxation time of the solvent is small enough to be comparable with the characteristic time of electronic excitations. For these cases, the strong frequency dependence of the dielectric function in the electronic excitation range reflects on the CFF. Such a CFF frequency dependence is a unique feature of nonequilibrium cavity field effects in the absorption spectrum that is unrelated to the anisotropy of the cavity.

To expose nonequilibrium cavity field effects, we compute the CFF for the two selected transitions (namely, A and B) for a set relaxation times and ascertain that (i) the CFF (CFFB and CFFA) and (ii) the CFF difference (CFFB – CFFA) change with the relaxation time. This dependence shows that there is an underlying frequency dependence in the CFF that cannot be accounted for by considering anisotropic cavity field effects.

Therefore, we performed further EOM-TDPCM/RT-TDDFT calculation on the solvated molecule for hypothetical solvents with shorter Debye relaxation times and with the same static and dynamic dielectric constant of water. For the sake of simplicity, we choose only one direction for the application of the impulsive field, namely, the x-axis (see Figure ). The results for the optical absorption spectra of solvents with τD = 3370 (water), 100, 10, and 1 fs are plotted in Figure for comparison.

Figure 3

Absorption spectra of trans-azobenzene in solution featuring cavity field effects or not for solvents with different relaxation times.

In Figure , we observe a great variation in the spectra featuring no cavity field effects for different relaxation times. In particular, we notice a systematic redshift of the first absorption peak (transition A) for diminishing relaxation times. This strong effect is due to reaction field and its discussion is out of the scope of this Article. We rather focus on the variation of CFF with the solvent relaxation time in the case of transitions A and B. A summary of these results is reported in Table .

Table 3

CFF Computed from Eq for Transitions A and B for Solvents with Different Relaxation Timesa

τD (fs)	3370	1000	100	10	1
CFFA	1.04	1.04	1.05	1.12	1.16
CFFB	1.26	1.26	1.26	1.30	1.70
CFFB–CFFA	0.22	0.22	0.21	0.18	0.54

The case of 3.37 ps is the water-like solution.

In Table , we observe a clear CFF relaxation-time dependence for both transitions, our first evidence that nonequilibrium cavity field effects are taken into account. In particular, we detect a monotonic behavior of the CFF with the diminishing relaxation time that can be explained qualitatively on the grounds of Onsager model. For any given frequency, Onsager CFF increases toward a maximum value for small relaxation times, that is, , and decreases toward a minimum value for large relaxation times, that is, . The ellipsoidal cavity model (vide supra) does not change the qualitative picture drawn from Onsager’s model for CFF versus τD; it simply tunes the maximum and minimum CFF at zero or infinite relaxation times and has excitation-dependent minimum and maximum values. The PCM, featuring a more complex molecular-shape cavity, further adjust quantitatively these maximum and minimum CFF values (as seen in Table ). In each model, the maximum and minimum CFF values for zero or infinite relaxation times are associated with the static and dynamic dielectric constants of the solvent, respectively. For the spherical model, the minimum and maximum CFF are independent from the frequency, whereas in the ellipsoidal model they depend on the orientation of the underlying transition dipole moment (anisotropic cavity effects). In the last row of Table , we show the difference CFFB – CFFA. As expected, the difference in CFF is not constant, which is another evidence of CFF dependence on the frequency unaccounted by anisotropic cavity effects. Therefore, nonequilibrium cavity field effects have an impact in the absorption spectral shape.

In Figure , we plot the CFF versus τD for the spherical (Onsager), the prolate ellipsoidal and the molecular-shape (PCM) cavity models. For the ellipsoidal model, we plot the overall CFF corresponding to the B band (before decomposed in B1 and B2), to facilitate its comparison with PCM results. Apart from the quantitative departure of the molecular-shape model from the rest, there is an important qualitative difference: the transient between maximum to minimum CFF values is narrower, occurring in a relaxation-time window of tenth of femtosecond for the molecular-shape cavity, as opposed to a few picosecond for spherical and ellipsoidal cavities (see Figure inset). The spherical and ellipsoidal cavity models do not predict relevant nonequilibrium cavity field effects for τD ≤ 0.5 ps or τD ≥ 3 ps, where the CFF takes approximately its maximum and minimum values, respectively. We recall that the latter might be transition-dependent for the ellipsoidal model because of the cavity anisotropy. It is in the relaxation time window 0.5–3 ps, where we expect to find nontrivial nonequilibrium cavity field effects in the form of strongly frequency-dependent CFF. However, in the case of a molecular-shape cavity, we find that nonequilibrium cavity field effects can only be relevant in the relaxation time scale of 0–100 fs. The latter disagreement highlights the fact that an accurate description of the molecular cavity is also essential to account for nonequilibrium cavity field effects, making impossible to decide when to apply our full EOM-TDPCM treatment (with cavity field effects) based on estimations extracted from simple models (such as Onsager’s).

Figure 4

CFF corresponding to transitions A and B of trans-azobenzene in solution vs the solvent relaxation time τD. Solid and dashed lines represent the spherical (Onsager) and prolate ellipsoidal cavity models, respectively, whereas points are associated with PCM results with a molecular-shape cavity. PCM data points are joined as a guide for the eye. Black (blue) color refers to A (B) transition. The inset shows the spherical and ellipsoidal models curves in a different scale.

Concluding Remarks

In summary, we have extended the time-dependent polarizable continuum model (TDPCM) framework,[2,5] in its equation-of-motion (EOM-TDPCM) version,[13] to account for nonequilibrium cavity field effects within quantum mechanical calculations of molecules in solution and subject to external electrostatic fields and/or electromagnetic radiation (light). Such a scheme (EOM-TDPCM with cavity field effects) provides a more realistic picture for the description of solvated molecule, shortening the gap between simulation/theoretical and experimental results.

Our theoretical contribution has been incorporated in the Octopus code (to be released as version 9.0).[33] This implementation enables us to perform real-time time-dependent density functional theory (RT-TDDFT) simulations coupled with EOM-TDPCM and optical absorption spectrum calculations, including nonequilibrium reaction and cavity field effects.

To show an application of our methodology, we compute the optical absorption spectrum of trans-azobenzene in water employing EOM-TDPCM/RT-TDDFT simulations. These calculations are relevant to highlight the nontrivial role of anisotropic and nonequilibrium cavity field effects in the optical absorption spectrum of solvated molecules, however, they are not intended to provide benchmarking results for trans-azobenzene in water. We underline that the cavity field factor, that is, the ratio between the absorption spectra with and without cavity field effects, changes throughout the UV−vis. spectrum frequencies depending strongly on the orientation of the TDMs of the excited states associated with each absorption peak.

Our cavity field effects description within PCM goes beyond the simplified Onsager’s model, where the molecule is hosted by a spherical cavity.[3,4] At difference with Onsager’s model, PCM give rise to a modification of the optical absorption intensities that is transition-dependent. This effect was rationalized in terms of another simplified model for an anisotropic cavity with ellipsoidal shape, which is able to reproduce qualitatively the PCM results.

Finally, we show that nonequilibrium solvent polarization affects the overall shape of the optical absorption spectrum. In particular, we find evidence of the nonequilibrium effects in the change of absorption peak height due to the cavity field for different relaxation times. We remark that an accurate description of the cavity shape is so important to our treatment of nonequilibrium cavity field effects that we cannot even rely on the typical Onsager’s (spherical cavity) model to conclude whether a nonequilibrium treatment is required or not.