# Simulations of coherent nonlinear optical response of molecular vibronic dimers

###### Abstract

We have implemented vibronic dynamics for simulations of the third order coherent response of electronic dimers. In the present communication we provide the full and detailed description of the dynamical model, recently used for simulations of chlorophyll-carotenoid dyads, terylene dimers, or hypericin. We allow for explicit vibronic level structure, by including selected vibrational modes into a ”system”. Bath dynamics include the Landau-Teller vibrational relaxation, electronic dephasing, and nonlinear vibronic (to bath) coupling. Simulations combine effects of transport and dephasing between vibronic levels. Transport is described by master equation within secular approximation, phase is accumulated in cumulants and its calculation follows the transport pathways during waiting time period.

## 1 Introduction

Two dimensional optical spectra [mukamel2000multidimensional, jonas2003two, cho2009two] of molecular aggregates shows almost invariantly important role of vibrations on relaxation and electronic transport [perlik2015vibronic] within organic dyes, light harvesters and similar molecular systems of intensive interest. Interplay of vibrational and electronic dynamics is a computational challenge which has been approached on the various level of theory [abramavicius2009coherent].

Description of nonlinear response [mukamel1999principles] of transporting multichromophoric electronic systems was pioneered by Zhang, Meier, Chernyak and Mukamel [zhang1998exciton]. Algorithm is based on partitioning bath induced fluctuations of excitons into ”diagonal part”, that is, frequency (eigenfrequency) fluctuations and off-diagonal, that is coupling fluctuations, the former accounted for by cumulants and responsible for lineshapes, the latter using some kind of reduced dynamics (master equations) and responsible for transport. This strategy after some refinements became standard, the choice of transport dynamics, bath models and further details of simulations, however differs among authors. More recent developments, and other strategies (important, but not relevant for the present communication) for simulating 2D spectra can be found in reviews [abramavicius2009coherent, mukamel2004many, falvo2013quasi].

Calculations along such a strategy work quite well, when the vibrational modulation are broad in frequencies and the spectra are thus without clear signatures of underdamped mode such as vibronic progression in the spectra. A better representation of underdamped mode requires to make them explicit, i.e. exclude important vibrational modes from bath and include them formally into the ”system” together with electronic degrees of freedom and calculate transport rates between the mixed electronic-vibrational levels. While direct implementation is costly and is thus of limited use for extended aggregates, it is accessible for small aggregates. After redefining the system, however, the rest of simulation strategy instituted by Ref. [zhang1998exciton] can remain largely untouched. Various variant of such a strategy has been simultaneously followed recently by several authors [butkus2014vibronic, butkus2013distinctive, polyutov2012exciton, schroter2015exciton].

In the present communication we report a dynamical model of vibronic dynamics affordable for small (dimers, trimers) molecular aggregates with resolved vibrational structure which we have implemented recently to inquire phase relation in vibronic systems [perlik2014distinguishing], and further used it for simulation of 2D spectra of carotenoid-chlorophyll dyad [perlik2015vibronic], and transient grating on hypericin [lincoln2015A]. Few others applications (perylene dimers) are yet in preparation [perlik2016perylene]. While the vibronic structures and system-bath dynamics were reported in the respective publications, there remains couple of interesting details of our implementation to be clarified to public. Last but not least, we have to report the full potential of our code yet.

So far, we have implemented and applied the vibronic dynamics only for molecular dimers. However, we are intending to apply vibronic dynamics to some trimeric systems. In addition, there would be no significant simplification of this manuscript if we limited the introduced formalism just to dimer. We thus define the vibronic dynamics for a general aggregate in sections 2 to 6, keeping in mind that our implementation was actually tested only for dimers so far.

The paper is structured as follows. In section 2 we introduce model of vibronic aggregate. System is separated from bath and diagonalized. In section 3 we introduce bath induced fluctuations. Dynamics include Landau-Teller vibrational relaxation, electronic dephasing and nonlinear bath-to vibration Hamiltonian. Diagonal eigenfrequency fluctuations are distinguished from off-diagonal eigenstate fluctuations. In section 4 we specify effect of eigenstate fluctuations and describe transport dynamics using secular time-convolutionless master equation. In section 5 we take care of eigenfrequency fluctuations and define their correlation and lineshape functions. In section 6 we calculate linear and the third order optical response. We show Feynman diagrams, and calculate second cumulants for eigenfrequency fluctuations along the diagrams. In our implementation we follow to some extent the reorganization of bath modes during transport. Section 7 discuss tensorial character of nonlinear response and rotational averaging of a response from randomly oriented chromophores in an isotropic sample. In section 8 we linked response functions to physical observable of interest: absorption, fluorescence, 2D electronic spectra, frequency resolved transient grating and pump probe are calculated. In section 9 we connect the previous sections with our recent simulations. In section 10 we conclude.

## 2 Vibronic System

We simulate nonlinear optical response of molecular aggregates with significant vibronic structure. Electronic structure of each molecule () of the aggregate is modelled by a two level chromophore (with ground and excited levels, and transition (gap) frequency ). For the third order response only a limited part of the composed electronic Hilbert space of aggregate is relevant. The relevant part consists of a ground state , one-exciton states where a single chromophore is excited and a doubly excited states . The molecules are further resonantly coupled; in a standard approximation the number of exciton is conserved and excitation of -th chromophore is associated with deexcitation of some -th chromophore. The Frenkel exciton (electronic) Hamiltonian is thus

(1) |

While the determination of transition frequencies does not pose significant problems, strategies for estimate of intermolecular coupling differs. It is approximated by dipole-dipole force, calculated using Quantum chemistry libraries or fitted from the spectra.

Quantum dynamics (defined by Eq. (1)) within the three electronic manifolds forms standard approaches for calculations of third order response. In the present work we have included certain vibrations as part of the system. Vibrations are assumed local. Each vibration is attached to some chromophore, is the coordinate of -th mode on -th chromophore. We assumed electronic potential surface to be harmonic with respect to nuclear coordinates , thus is the ground state’s potential and is the excited states’ potential surface, where is the vibrational frequency, is the mass and is displacement. Vibrational Hamiltonian reads

(2) |

where are creation (annihilation) operators of harmonic vibrational mode.

Vibrons are complex electronic-vibrational excitations. To describe them we have to define system as composed from electronic and vibrational degrees of freedom. We thus works in composed Hilbert space and the system molecular Hamiltonian is

(3) |

We next introduce convenient basis. We start with basis of vibrational states over site. It consist of ground state electronic wavefunctions , with a well-known wave functions of harmonic oscillators

(4) |

and excited state electronic wavefunctions with a wave functions of displaced harmonic oscillators

(5) |

where is the vibrational ground state on electronic ground state and is the (shifted) vibrational ground state on electronic excited state. In absence of coupling Hamiltonian (3) is diagonal in the product basis where is some vector listed in Eq.(4) or Eq. (5) .

In a general situation molecular Hamiltonian (Eq. (3)) is no longer diagonal in this basis, however the basis still forms a convenient starting point for numerical implementations. Indeed, the matrix element for coupling terms are composed of the widely-known Franck-Condon factors and can be readily implemented along with diagonal part of the rest of Hamiltonian, and subjected to a standard routines for numerical diagonalization.

The diagonalization of the full vibronic Hamiltonian must be thus, in general, diagonalized numerically.

(6) |

Eigenstates and eigenfrequencies will be hereafter indexed by Greek letters. Hamiltonian Eq. (3) conserves number of electronic excitations. Each exciton manifold can be thus diagonalized separately.

On the electronic ground state manifold, there is no resonant coupling, thus the eigenfunctions are direct product of ground state electronic wave function with a well-known wave functions of harmonic oscillators introduced in Eq (4). Corresponding eigenenergies are . The diagonalization of the single excited electronic manifold is more complex, eigenstates are always obtained by a numeric diagonalization of the corresponding block of the Hamiltonian. For double excited electronic manifold the situation is different for dimer and longer aggregates. In a special case of dimer, doubly excited states does not allow transport (and resonance coupling) and eigenstates are direct product of doubly excited electronic wave function with a wave functions of displaced harmonic oscillators with eigenenergies . For larger aggregates the diagonalization of the second manifold is again numerical. Third and higher manifolds, which may appear in larger aggregates, do not enter calculation of the third order response [abramavicius2009coherent].

Interaction with the probing laser fields will be treated in a dipole approximation and in Condon approximation described using interaction Hamiltonian

(7) |

where ; and is a transition dipole moment between the ground and the excited state of -th molecule. Its matrix elements

will be used hereafter. For the section 3-6 we neglect the vector structure and consider dipole moment as a real number. This approach is sufficient for describing the optical response from most common aggregates with (anti-) parallel dipoles. The neglect will be cured in section 7, where the tensorial structure will be calculated together with rotational averaging of the response.

## 3 Interaction with bath

The environmental fluctuations modulating system Hamiltonian are mainly of solvent origin or come from less important vibrations. These are modelled by a dense set of bath harmonic oscillators with bath Hamiltonian

(8) |

where is a frequency of a -bath mode, and are standard bosonic creation and annihilation operator, respectively.

System-bath interaction is responsible for damping of vibrational modes and electronic dephasing. Three different couplings are included into

(9) |

where the three terms are responsible for Landau-Teller vibrational relaxation, non-linear vibronic-bath couplings, and for electronic dephasing, respectively.

(10) |

Note that and are constant and linear, respectively in vibrational coordinate and are thus within the standard spin-boson model. Such a dynamics can thus be directly compared to standard simulations, they use the same Hamiltonian, but define system and bath differently. In contrast, the nonlinear coupling takes us beyond the standard spin-boson model, and there is no correspondence to standard simulations with all vibrations included in bath.

Following the strategy of [mukamel2004many] system-bath Hamiltonian shall be divided

into the diagonal and off-diagonal part of fluctuations responsible for lineshapes and exciton transport respectively. To that end we shall transform system-bath Hamiltonian Eq. (9) into eigenbasis (Eq. (6)) and isolate the diagonal part representing eigenenergy fluctuations to yield

(11) |

Similarly the off-diagonal system-bath Hamiltonian representing eigenstate fluctuations reads

(12) |

where we defined the matrix elements

The third-order nonlinear optical response for a system with only diagonal system-bath Hamiltonian could be obtained by using a second cumulant expression. The off-diagonal fluctuations shall be treated by means of a master equation. We followed the strategy that combines the two. Before we give the formula to be implemented we first prepare the relevant master equation.

## 4 Eigenstate fluctuations: master equations

The off-diagonal fluctuations will be accounted for by using time convolution-less master equation

(13) |

where is reduced (system) density matrix and stands for a relaxation tensor. The relaxation tensor shall be evaluated to the second order in system-bath coupling. As argued by Redfield [redfield1957theory] only certain terms called secular and obeying

(14) |

contribute significantly. For aggregates, the eigenfrequencies are not degenerated or systematically built, so that Eq. (14) can only be obeyed when either and which term describes population transport, or and which term describe coherence decay. We adopt this secular approximation on relaxation tensor in (13) leaving only terms representing population transfer or coherence dephasing and neglect coherence transfers, coherence to transfer terms etc. In addition to Redfield argumentation based on relevance, we note that such a choice significantly reduces number of Feynman diagrams involved in calculation of coming sections.

The relaxation tensor is calculated in Tokuyama-Mori formalism [tokuyama1976statistical]. We choose form of Ref [vcapek2001violation, vcapek1994interplay] in asymptotic limit and evaluate it to the second order in .

(15) |

where represent secular approximation and where the following superoperator notation is introduced:

(16) |

(17) |

(18) |

(19) |

where

(20) |

is canonical bath density matrix.

Evaluating (15) yields :

(21) |

is relaxation tensor element responsible for populations intra-manifold dynamics between states and . When we have . Element represents decoherence

(22) |

where is the Bose-Einstein distribution

(23) |

and where we defined spectral densities

(24) |

(25) |

(26) |

In Eqs (21) and (22) we have standardly neglected cross terms . Spectral densities Eq. (24)-(25) can be arbitrary positive functions defined on positive semi-axis. Successful estimation of spectral densities from microscopic foundations, i.e. molecular dynamic simulations are rare [olbrich2010time]. In practice there is handful of popular forms for spectral densities which are parameterized by fitting them to spectra or transfer rates.

For our simulations we used spectral densities of overdamped Brownian oscillator (exponential decay)

(27) |

(28) |

(29) |

where is a relaxation rate and is a reorganization energy. denotes Heaviside’s step function.

## 5 Eigenfrequency fluctuations: Lineshape functions

We next inquire the effect of diagonal part of system-bath interaction Hamiltonian. These are responsible for eigenfrequency represented by bath-space operator

(30) |

Our -linear coupling introduces Gaussian fluctuations, which can be fully characterized in terms of matrix of correlation function

(31) |

where with define in (20). Inserting (30) into (31) yields

(32) |

We can calculate (32) in terms of spectral densities

(33) |

where

(34) |

The effect of Gaussian noise on lineshapes is infamously given by line-broadening function which can be obtained by a double integration of

(35) |

Now we recall our spectral functions , , posses form of a Brownian spectral density and we obtain well-known result for overdamped Brownian oscillator. For :

(36) | ||||

(37) | ||||

(38) |

where are Matsubara frequencies. To get cumulant at negative times, we use the relation .

## 6 Liouville space pathways

Description of exciton population transport by using master equation requires to switch formally into the Liouville space [mukamel1999principles]. To that end we introduce superoperators (L) operating on ket (e.g. ) and (R) operating on bra () index of density matrix. Expanding commutator as the linear response function reads

(39) |

and six contributions for photon echo and free induction decay signals of third order response are

(40) |

where is the evolution superoperator. Equations (40) are still quite formal since the evolution superoperator must include both the master equation and the dynamics of diagonal fluctuations, what is complex task for pathways changing its exciton index in population transfer during a waiting time period. We approximate as follows

(41) |

where is evolution of bath in eigenstate written in interaction picture and is Green function of master equation (13).

With use of Eq. (39) linear response function reads

(42) |

where is total density in state and we merged bath phase factor into

(43) |

The equilibrium density matrix is defined as

(44) |

where we assumed that the process starts at chromophore’s electronic ground state with vibrational level populated according to Boltzman distribution. Using Gell-Mann’s theorem [gell1951bound] we can obtain the equilibrium density matrix by switching the interaction at time where we start from uncorrelated density matrix , i.e.

(45) |

where the index runs over vibrational levels of electronic ground state , was approximated . That leads us to