Linear optical properties of organic microcavity polaritons with non-Markovian Quantum State Diffusion

Hybridisation of the cavity modes and the excitons to polariton states together with the coupling to the vibrational modes determine the linear optical properties of organic semiconductors in microcavities. In this article we compute the refractive index for such system using the Holstein-Tavis-Cummings model and determine then the linear optical properties using the transfer matrix method. We first extract the parameters for the exciton in our model from fitting to experimentally measured absorption of a 2,7-bis [9,9-di(4-methylphenyl)-fluoren-2-yl]-9,9-di(4-methylphenyl) fluorene (TDAF) molecular thin film. Then we compute the reflectivity of such a thin film in a metal clad microcavity system by including the dispersive microcavity mode to the model. We compute susceptibility of the model systems evolving just a single state vector by using the non-Markovian Quantum State Diffusion. The computed location and height of the lower and upper polaritons agree with the experiment within the estimated errorbars for small angles ($\leq 30^\circ$). For larger angles the location of the polariton resonances are within the estimated error.


I. INTRODUCTION
The localized nature of molecular excitons imparts two crucial attributes: a large binding energy and a dense population of excitons compared to the photonic density of states.The former feature promotes strong lightmatter coupling at room temperature, making possible the macroscopic study of Bose-Einstein condensation and superfluidity at elevated temperatures [1][2][3][4].The latter feature enabled the observation of ultrastrong coupling [5] and single photon nonlinearity [6].Recently, the concept of polariton chemistry was introduced, which promises to reshape the energy landscape of molecular systems, exerting control over their photochemical and photophysical processes [7][8][9][10][11].In optical microcavities filled with molecular absorbers, polariton modes are the eigenstates resulting from the strong coupling between the optical modes and the exciton resonances.They can be observed experimentally as the avoided crossing of the bare exciton and microcavity photon dispersion through reflectivity measurements [12].
An interesting implication of strong coupling at nonzero temperature is that the exciton might also strongly couple to the surrounding molecular vibrations which can be the catalyst for complex relaxation dynamics of polaritons [13].Currently, there is significant research being conducted to explore the impact of strong collective light-matter interactions on the intersystem crossing (ISC) and reverse intersystem crossing (RISC) rates [14][15][16][17][18][19][20][21].Rabi splitting is proportional to the square root of the molecular density that share the microcavity photon ( N/V ).To achieve strong coupling in planar microcavities with highly delocalized photonic mode, a large number of molecules N , in the order of 10 6 is needed.However, perturbative quantum mechanical calculations predict that possible polaritonic contributions to the RISC or ISC processes, for example, scale with inverse of this number [22][23][24][25].
The purpose of this article is to introduce a new method for investigating the polariton dynamics and pave the way for gaining new insights into the so called "large number of molecules N " problem.In this work we take the first steps of applying non-Markovian Quantum State Diffusion (NMQSD) for computing the linear optical properties of organic microavity polaritons.While a possible drawback of this approach is that it is stochastic and typically requires a large number of trajectories for computing expectation values of observables, for computing linear optical properties using just one trajectory is enough [26][27][28].
In Fig. 1 we present a graphical summary of this article.a) We use NMQSD to compute the susceptiblity of the model systems evolving just a single state vector.From the susceptibility we obtain the refractive index which we use to model the experimental situations studied in this article.b) We compute the absorption of a thin film of TDAF molecules.By fitting the model to the experiment we obtain the parameters of the exciton.c) Lastly, we investigate a situation where the TDAF molecules are placed within the microcavity and we compute the reflectivity of such a system.
The structure of this article is the following.In Sec.II we introduce the linear response theory and how the susceptibility and refractive index can be computed.In Sec.III we introduce the NMQSD method.Furthermore, we show how the susceptibility can be computed using the NMQSD approach.Then in Sec.IV we fit the model parameters to the single molecule data obtained from experimental measurements and density functional calculations.We will focus on 2,7-bis [9,9di(4-methylphenyl)-fluoren-2-yl]-9,9-di(4-methylphenyl) fluorene (TDAF) as it is a model system for strong lightmatter studies [21,29,30].In Sec.V we construct the model for microcavity polaritons using the same TDAF molecule.We discuss in detail how the susceptibility can FIG. 1. a) We model the linear optical properties of a slab of material in terms of the dielectric function ε, which can be computed from the susceptibility χ.The susceptibility is computed from a quantum mechanical model evolving just a single state vector.This leads to a computationally efficient scheme in contrast to approaches where density matrix evolution is needed.b)Susceptibility of a thin film can be modeled as the susceptibility of a single quantum absorber (red box) multiplied with the the number of such absorbers, when spatial disorder can be neglected.We include energetic disorder as indicated by the distribution of exciton energies εs, and a coupling to vibration models when the system is excited indicated by the spring.c)In the case of microcavity polaritons we consider a slice of emitters in the z-direction (red box), whereas the cavity mirrors are located at x = ±L/2 (top and bottom).
be computed in this case when the dispersive microcavity mode is also included.We also present the results of the theoretical calculations.In Sec.VI we present our conclusions and outlook.Experimental details are presented in Sec.VII II.LINEAR RESPONSE a. Linear optics In linear materials the polarization field is proportional to the applied electric field where ε 0 is the vacuum permittivity and χ is the susceptibility.We set ℏ = 1 in all subsequent equations.Due to causality polarization can depend only on the fields applied on earlier times and may have non-local spatial dependency [31].In general, χ is a second rank tensor.
In this work we focus only on situations where the polarization is aligned with the applied electric field making χ a scalar.By using the convolution theorem and after dropping the vector notation this relation is where ω is the angular frequency and k z is the zcomponent of the wave vector.The dielectric function (or relative permittivity) is The dielectric function determines the refractive index by n = √ ε.Once the refractive index is known, the linear optical properties of planar systems can be computed using the transfer matrix method (TMM) [32].
b. Dipole density Polarization corresponds to the dipole density of the medium, which can be computed from a microscopical model.We consider a situation where weak classical field is used to probe a quantum system.The interaction term between the field and the matter is taken to be where z m is the location of the dipole and mu m is the dipole operator.The other degrees of freedom are described by a time independent Hamilton operator H.The dynamics generated by H is given by the unitary operator The linear response can be computed from the dipole correlation function M (t) [26,27,33].The susceptibility of the system can be computed as the Fourier transform of the dipole correlation function M (t) where and |Ψ G ⟩ is the ground state of the Hamiltonian H.We assume that the systems we investigate do not have permanent dipole moment.Under this assumption the dipole correlation function is obtained from perturbation theory and keeping only the positive energy terms and the terms linearly proportional to the applied field [34,35].Generalization to anisotropic cases is straightforward.Using Eqs. ( 2), (6) gives the polarization density of the system.The macroscopic polarization is obtained by multiplying the polarization density with the sample volume.In this work we neglect any spatial disorder.Finite temperature results can be computed using the ground state initial condition in the NMQSD approach as we show later.

III. NON-MARKOVIAN QUANTUM STATE DIFFUSION
a. General theory The aim of the NMQSD approach is to solve the time evolution of the full Schrödinger equation for the open system and the environment [36][37][38].A typical model consists of an open system with Hamiltonian H S and a coupling operator L. These operators are arbitrary at this point.The environment is assumed to consist of quantum harmonic oscillators with a Hamiltonian where Generalization to spin baths is possible [39].The interaction is taken to be linear in the coupling operator L of the open system and the creation and annihilation operators of the environment The initial state of the bath is the thermal state ρ β and the system and the bath are initially uncorrelated.In the interaction picture with respect to (8) the Schrödinger equation is The finite temperature NMQSD equation corresponding to the Schrödinger equation ( 10) is [36] where z * t is a zero mean Gaussian stochastic process with correlations The hermitian autocorrelation of the process corresponds to the zero temperature bath correlation function (BCF) We have introduced the spectral density which controls the properties of the environment.
Finite temperature is incorporated by a "stochastic potential" V (t) where η t is a zero mean Gaussian stochastic process with correlations [37,40] and n λ = (e βω λ − 1) −1 is the thermal photon number.
The NMQSD equation is a stochastic differential equation containing the Hamiltonian term, stochastic driving term and a memory term with a functional integral.The states |ψ t (z * )⟩ are analytical functionals of the noise process , where z * λ are the labels of the Bargmann coherent states of the environment [38].By construction the exact open system is recovered by averaging over the stochastic trajectories The challenging part in using NMQSD is the occurence of the functional derivative.In recent years a powerful hierarchy of pure states (HOPS) approach has emerged as a general numerical approach to solve the NMQSD equations [41].However, in this work we use a different approach.Namely, we approximate the functional derivative with the following expression This approximation is obtained from the general HOPS approach by truncating the hierarchy after the first level and corresponds to a weak coupling approximation (we neglect any terms higher than second order in the coupling strength) [37,41,42].b.Susceptibility using NMQSD.Using the NMQSD we can compute the susceptibility of the system evolving pure states only.In case the coupling operator is hermitian we need to evolve only one pure state, otherwise we need to average over the thermal noise [28].In the case that we have multiple transition dipole moments we need to extend the NMQSD to many systems which all couple to their individual environments.This simply means that each subsystem has their own coupling operator L m , noise term z * t,m and bath correlation function α m (t).We assume that the Hamiltonian is such that the global ground state is a product from the system ground state and the vacuum of the bath , where we consider M systems and |0⟩ is the vacuum state of the environment of the respective open system.In the NMQSD approach the dipole correlation function M (t) can be computed from the time-evolution [26][27][28] where |ψ 0 ⟩ is the initial state for the NMQSD evolution and |ψ t (z * = 0)⟩ is the solution to the NMQSD equation where the driving noise is set to zero, i.e. z * t = 0.The initial condition for the evolution is chosen as If the coupling operator is hermitian, we can replace the stochastic potential with the finite temperature bath correlation function Otherwise we need to average over different realizations of the thermal noise [26].
IV. MOLECULAR THIN FILM a. Description.The exciton of the TDAF is determined from experimentally measured absorption of a 60nm-thick film of TDAF on a quartz substrate.The absorption is defined as A = 1 − T − R, where T, R are the fractions of the transmitted and reflected light, respectively.To accomplish the measurement of the reflected and transmitted light from the film without increasing the optical path length, the sample was excited at a 15 • angle.We will model this process by computing the refractive index from a microscopic model and then the reflected and transmitted light using TMM [43].
b. Model.The dynamics of the molecule is governed by the Holstein model We denote by K = σ + σ − from now on.The system Hamiltonian, and the coupling operator are in this case where ζ is a disorder parameter.The coupling operator is hermitian.The NMQSD equation in this case is where α(t − s) is the thermal BCF.We model the radiative damping by and additional white noise process ξ * t with correlation M [ξ t ξ * s ] = γδ(t − s).For computing the susceptibility we can set both noise terms to zero.The spectral density is taken to be superohmic with exponential cut-off [23] where a is parameter additionally controlling the coupling strength.In the limit that the radiative damping is small compared to other parameters of the system the model admits a solution where We set the integration limits in this way to remove boundary terms and the reorganization energy term, which we absorb to the singlet energy.We do not highlight explicitly in the notation that the noises ξ * t = z * t = 0 anymore.The disorder ζ is distributed according to a Gaussian distribution with zero mean and standard deviation σ.
c. Susceptibility.The dipole operator for this system is with the abuse of notation we use the same symbol for the transition dipole operator and the transition dipole moment.The initial state to use is the ground state of the molecule |g⟩.Inerting this ground state and computing the average over the disorder gives the following expression for the dipole correlation function In the case that g(t) = 0 this corresponds to the Voigt lineshape.The susceptibility is obtained by taking the Laplace transform from the dipole correlation function (27) and the refractive index can be readily computed.d.Thin film absorption.The first singlet excited state is of the system is at ε S ≈ 3.6 eV as can be seen from the experimental trace in Fig. 2. The radiative lifetime of the TDAF thin film is reported to be 133 ps (∼ 10 −5 eV) [44].When fitting the model to the data we ensure that the standard deviation of the disorder parameter is larger than the radiative life time so that the solution (25) remains valid.For the above mentioned parameter values the thermal contributions to the dipole correlation function are insignificant and we use the zero temperature BCF when fitting the model to the experimental data.
The parameter values found in the fitting process are ε s = 3.6 eV, σ = 0.14 eV, and ξ = 0.09.We kept the coupling strength parameter at a fixed value a = 1 and γ = 5 × 10 −5 eV and included a background refractive index n bg = √ 1.5 + 0.015i to model the residual absorption at small energies.The reorganization energy for the fitted parameters is λ s ≈ 0.19 eV which is agreement with the values reported in the DFT calculations [21].Absorption of a 60 nm thick thin film of TDAF molecules.The exciton energy is approximately 3.6 eV.Our model (NMQSD) fits well to the measured data, whereas a blind fit with a Voigt lineshape performs weaker in line with reported χ 2 values for each fit.The errorbars in the fit are smaller than the linewidth.
The reorganization energy is absorbed in ε s .The measured data, the fit using our model and a blind fit to a Voigt lineshape are shown in Fig. 2. The χ 2 N M QSD is significantly smaller than the χ 2 V oigt .We obtain the Voigt lineshape by neglecting the molecular vibrations contained in the term g(t) in Eq. ( 27).This term is responsible for asymmetric broadening of the lineshape on higher energies.e. Summary.We model the thin film of TDAF molecules as two level systems which each are coupled to their respective molecular vibrations and are probed by a weak classical electromagnetic field.We assume that there is no spatial disorder but the energies of the excitons are distributed according to Gaussian distribution with mean ε s and variance σ.The macroscopic polarization is then obtained by multiplying the induced dipole moment with the number of molecules in the sample.The inclusion of the disorder is motivated by the likeness of the experimental absorption lineshape to a Voigt profile and the fact that we found that thermal effects in our model were negligible in this parameter regime.The asymmetry in the absorption lineshape is explained by the coupling to the vibrational bath with spectral density given in Eq. (24).We find the parameter values by fitting the model to the experimentally measured absorption using TMM.The fitted parameters are the excitonic energy (ε s ), energy disorder (σ), coupling strength a and the cut-off ξ.We keep the radiative lifetime of the exciton and the background refractive index constant during the fit.The fit of our model is very good and the obtained parameter values are in agreement with what is obtained from DFT calculations [21].

V. MICROCAVITY POLARITON REFLECTIVITY
a. System.The system is a 80 nm TDAF film in a cavity formed by two aluminium mirrors with thicknesses 25 nm and 100 nm.The reflectivity of this system can be measured experimentally as a function of the angle and energy of the incoming light.We again compute the refractive index from a microscopic model and then compare the reflectivity calculated using TMM with the experimental data.
b. Model.We follow [22,45] in the construction of the model.The system consists of N molecules in a planar microcavity.The Hamiltonian for the molecules is Photons in the cavity have the energy where c is the speed of light in the vacuum and n r is the refractive index of the propagation medium.k is the wave vector of the light which is assumed to be in the x-z plane.Mirrors at x = ± L 2 confine the light in the x-direction so that the wavevector k x = mπ/L, where m ∈ N. The cavity frequency at zero incidence is ω 0 = mπc/(n r L).We restrict ourselves to m = 1 case, so that the cavity dispersion relation is Each mode has two orthogonal transverse polarizations u 1 and u 2 = k×u1 |k| .The Hamiltonian for the cavity modes is then where a kz is the annihilation operator for the cavity mode with the wave vector z-component k z .The cavity modes couple to the transition dipole moment µ j with coupling strengths where u indicates the direction of the electric field of the confined mode, ε 0 is the vacuum permittivity and V the cavity mode volume.We assume that the polarization of the cavity modes u, the dipole moments µ j and the polarization of any incoming light are all aligned in the y-direction.Considering a system with polarization parallel to the applied electric field, see Eq. ( 2).The coupling between the molecules and the cavity modes is given by the Tavis-Cummings interaction where e ±ikzzj describes the phase of the electric field of the cavity mode at the position z j of the molecule j.
The coupling of each molecule to local vibrational modes is the same as in Sec.IV where we assume that each molecule couples to its own bath of vibrational modes.The total Hamiltonian for the system is then where H E is the free Hamiltonian for the vibrational modes In addition the cavity modes are damped with a rate that does not depend on k z and denote it by κ.The cavity damping is described in terms of additional zero mean white noise processes w t,kz with correlations Similarly, the excitons are damped with the rate γ, which is described by white noise processes ξ t,j .The evolution of the system is given by the NMQSD equation which is trivially extended from the ones given earlier in the manuscript.Namely, we read off the terms of the NMQSD equation from the HTC Hamiltonian (Eq.( 35)), add independent excitonic and cavity dampings described by the white noise terms (Eq.( 37)).Differently to the previous case, the thermal effects are significant and we use the non-zero temperature bath correlation function (Eq.( 20)) in the subsequent calculations.c.Susceptibility.We write the transition dipole moments for the molecules similarly as in Eq. (26) for each individual molecule.The dipole correlation function of this system can be calculated using (18).In this case the initial state ( 19) is used.We take into account the possibility of light being absorbed by the cavity mode by an additional dipole moment µ N +1 (k z ) = µ c a kz + µ * c a † kz with z N +1 = 0 reflecting the fact that the cavity mode is delocalized in the whole cavity volume.We can then calculate the linear susceptibility of the system from equation ( 6) and the refractive index.We use the refractive index to calculate the reflectivity of the cavity system by using TMM for incoming light with the angular frequency ω and zcomponent of the wave vector k z .These are related to the angle θ of the incoming light by where θ c is the angle inside the cavity and we have used the Snell's law.d.Reflectivity.We show the computed and measured reflectivity in Fig. 3 as a density plot.The quantitative agreement is good at small angles.For larger angles the locations of the polariton modes are in good qualitative agreement.In Fig. 4 we show the computed and measured reflectivity as a function of the energy for different angles.There the disagreement at larger angles is more prominently visible.In Fig. 4 we also present the estimated errorbars.We can conclude that for angles up to approximately 30 • the model and the experiment are within the estimated errorbars.We discuss the errors involved in Sec.VI.In the TMM calculations the front mirror and the TDAF film is replaced with a material that has the computed refractive index.We use the parameters found in the fitting process in Sec.IV.However, we set molecular disorder to zero and discuss this point in Sec.VI.The number of molecules N = 30 and the number of cavity modes is 21.The cavity decay rate κ = 0.21 eV is estimated from experimental photoluminescence of the cavity system.We fitted the positions and depths of the reflectivity minima to the experimental data and got the values n r = 2 for the refractive index that determines the cavity dispersion relation (30), E(0) = 3.42 eV for the energy of the cavity mode with k z = 0, Ω = 0.92 eV for the Rabi splitting, and µ c = 2µ j for the magnitude of the cavity dipole moment.At small angles the lower polariton is mostly photonic and the upper polariton is mostly excitonic.Without vibrational coupling and cavity mode absorption this would lead to pronounced absorption of the upper polariton in comparison to what is observed experimentally.Including quantum mechanical coupling to the vibrational modes and taking in account the possibility for cavity mode absorption leads to an agreement between the experiment and our theory.If the cavity mode absorption is not included in (38), the upper and lower polaritons absorb roughly the same amount, which is not what we see in the experiment.
e. Summary.We extend the thin film model from Sec. IV by including the dispersive cavity mode and spatial locations of the molecules.We have excluded the energy disorder of the excitons from the model.This is done because the effect of the energy disorder led to poor agreement with the model and the data.It may be that the approximations we do in our modeling do not correctly take into account the disorder and the coupling to the cavity modes.On the other hand, neglecting the effects of the exciton broadening are in line with the observations [46,47], where the coupling to the cavity mode diminishes the exciton broadening effect.There, however, the cavities had higher Q factors as in our case.We estimate the error in our modeling by Monte Carlo sampling.We sample the computed reflectivity with randomly choosing the input parameters.We assume that all of the deviations are independent and distributed normally around the parameter values used in Fig. 3.We consider the following deviations: The input energy has a standard deviation σ E = 10 −6 eV, the input σ θ = 0.001 • and the thickness of the system σ d = 10 −9 m.We assume that the computed refractive index has energy dependent standard deviation σ n = 10 −2 e −(ω−ω0)/2 for the real and imaginary part, where ω 0 is the smallest energy used.The deviation is smaller for higher energies, which necessary for obtaining meaningful errorbars around the polariton energies.The error bars in the location of the polaritons (x-axis) is taken to be the average variance of the estimated variances of the upper and lower polariton peak locations.

VI. CONCLUSIONS AND OUTLOOK
We have computed the susceptibility for organic microcavity polaritons from the Holstein-Tavis-Cummings model using non-Markovian Quantum State diffusion.This approach is very efficient since we can compute susceptibility from evolving just a single state vector, so that density matrix computations are not needed.We have shown that our model can explain the absorptive properties of TDAF thin films and the reflectivity of the TDAF microcavity polaritons.In the thin film case the model explains the asymmetric shape of the absorption line very well.In the polariton case we do have good quantitative agreement around the polariton peaks at small angles (≤ 30 • ) and good qualitative agreement at larger angles.We verify the modeling results with estimated error bars.
We have identified several open questions that will be the focus of our forthcoming investigations.We have seen that energetic disorder leads to a good fit in the thin film case but needs to removed in microcavity polariton case: i) Possible explanations used in the literature are the cavity filtering effect [46,47] or motional narrowing [48,49].The former may not be the right explanation in this case since the cavities used in this work have such poor Q factors.The latter explanation fails as the spectroscopy we use conserves the planar wave vector of the cavity [31].
ii) In this work we have relied on perturbative solutions to the NMQSD equation as they provided reasonable fits to the experiment.It will be interesting to investigate situations where such perturbative approaches fail.In such cases, there may be more complex intramolecular dynamics which may involve also the spin orbit coupling between the singlet and triplet states [23].
iii) Lastly, we point out that state vector based approaches, such as NMQSD, open up a new way to study delocalization degree of the polaritons and polariton transport as it is possible to observe dynamically how the localization due to coupling to vibrational degrees of freedom and delocalization due to cavity coupling compete.

A. Fabrication
The samples were fabricated using thermal evaporation at a base pressure below 10 −7 Torr (Angstrom Engineering physical vapor deposition system).We used 15×15 mm 2 quartz substrates that were cleaned by sonication for 10 min in soapy water (3 % Decon 90), acetone, and isopropanol, respectively and dried with nitrogen.A 100-nm-thick aluminium was deposited on top of the substrate as the bottom mirror, followed by the deposition of 80 nm TDAF as the active layer, 1 nm of LiF, and a 25nm-thick aluminium layer as a top polariton microcavity mirror.

B. Characterization
The TDAF absorption and polariton angle-resolved reflectivity were measured with a spectroscopic ellipsometer (J.A. Woollam VASE) in reflectivity and transmission configuration.To extract the absorption of TDAF film, we measured transmitted and reflected light at a 15 • excitation angle, which represents the minimum angle our setup can measure reflectivity and adds only a 2% increase in the optical path.
FIG. 2.Absorption of a 60 nm thick thin film of TDAF molecules.The exciton energy is approximately 3.6 eV.Our model (NMQSD) fits well to the measured data, whereas a blind fit with a Voigt lineshape performs weaker in line with reported χ 2 values for each fit.The errorbars in the fit are smaller than the linewidth.

FIG. 3 .FIG. 4 .
FIG.3.Experimental and calculated reflectivities for the cavity system.The qualitative agreement of the experiment and the theory is good.The reduced reflectivity of the upper polariton is due to the coupling to the vibrational modes.The agreement between our theory and the experiment is better at smaller angles.