Linear optical properties of organic microcavity polaritons with non-Markovian quantum state diffusion

2.1 Linear optics

In linear materials the polarization field is proportional to the applied electric field

where ɛ ₀ 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 [45]. 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 z-component 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) [46].

2.2 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 field is used to probe a quantum system. The interaction term between the classical electromagnetic field and the matter is taken to be

where z _m is the location of a point dipole and μ _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) [39], [40], [47]. The susceptibility of the system can be computed as the Fourier transform of 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 [48], [49]. Generalization to anisotropic cases is straightforward. Using Eqs. (2) and (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. Alternative methods for solving Eq. (7) exist such as cumulant expansion techniques [33], [34] or exact stochastic wave function sampling [35].

3.1 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 [50], [51], [52]. 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 b λ , b λ ′ † = δ λ , λ ′ . Generalization to spin baths is possible [53]. 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 [50]

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 J ( ω ) = ∑ λ g λ 2 δ ( ω λ − ω ) , 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 [51], [54]

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 open system states |ψ _t (z*)⟩ = ⟨z‖Ψ _t ⟩ are analytical functionals of the noise process z t * = − i ∑ λ g λ z λ * e − i ω t , where z λ * are the labels of the Bargmann coherent states | | z 〉 = ∏ λ e z λ a λ † | 0 〉 of the environment [52]. By construction the exact open system is recovered by averaging over the stochastic trajectories

The NMQSD approach works in general with an arbitrary bath correlation function. For example, bath correlation functions obtained from surrogate harmonic bath embeddings such as in Ref. [38] could be incorporated into our formalism. The challenging part in using NMQSD is the occurrence 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 [51], [55]. We explain next the weak coupling approximation for the functional derivative term we use later in this work. The weak coupling means that we keep terms up to quadratic order in the coupling strength g _λ . In particular we have α ( t ) ∼ O g λ 2 and z t * , η t * ∼ O ( g λ ) . The HOPS equations are obtained when the BCF is expanded into a sum of exponentials

and defining auxiliary states | ψ t μ ( z * ) 〉 = ∫ 0 t d s G μ e − W μ ( t − s ) δ δ z s * | ψ t ( z * ) 〉 . The NMQSD equation now reads

and the auxiliary states satisfy following equation of motion [55]

We truncate hierarchy to first order, which means that we neglect the last term. We also neglect the term V(t) and z t * as they are proportional to g _λ . With these approximations we find

As G _μ are already proportional to g λ 2 we can replace | ψ s ( z * ) 〉 ≈ e i H S ( t − s ) | ψ t ( z * ) 〉 under the integral sign, which we obtained by integrating the NMQSD equation from backwards in time from t to s and keeping only the zeroth order terms. By considering the sum over all of the terms in the exponential expansion we obtain the weak coupling approximation

We stress that these approximations lead to closed NMQSD equations of motion and correspond to a weak system bath coupling approximation [56]. Further refinements are possible in terms of HOPS, but we could explain the experimental data in this work using the weak coupling approximation only.

3.2 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 [41]. 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 |Ψ _G ⟩ = |g ₁, 0⟩|g ₂, 0⟩ … |g _N , 0⟩, where we consider N 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 [39], [40], [41]

where |ψ ₀⟩ 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 [39].

4.1 Description

The exciton of the TDAF is determined from experimentally measured absorption of a 60-nm-thick film of TDAF on a quartz substrate. The absorption is defined as A = 1 − T − R, where T and 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 [57].

4.2 Model

The dynamics of the molecule is governed by the Holstein model [26]

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 ) . This model is not anymore solvable in closed form. If we neglected the radiative damping term, the model would admit an analytical solution [27]. For computing the susceptibility we can set both noise terms to zero. The spectral density is taken to be super-Ohmic 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 g ( t ) = ∫ − ∞ t d s ∫ − ∞ s d u α ( u ) . 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 σ. As can be seen, the analytical solution is obtained by using Eq. (17). It turns out that this ansatz is exact for the Holstein model and also a good approximation for our model in the parameter regime relevant for this work.

4.3 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. Now the initial state given by (19) is the excited state |e⟩ of the molecule. Inserting this 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.

4.4 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 Figure 2. The radiative lifetime of the TDAF thin film is reported to be 133 ps ( ∼ 1 0 − 5 eV) [58]. 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.

Figure 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 χ ² values for each fit. The errorbars in the fit are smaller than the linewidth.

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⁻⁵ eV and included a background refractive index n b g = 1.5 + 0.015 i to model the residual absorption at small energies. The reorganization energy for the fitted parameters is λ _s ≈ 0.19 eV which is in agreement with the values reported in the DFT calculations [21]. 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 Figure 2. The χ NMQSD 2 is significantly smaller than the χ Voigt 2 . 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.

4.5 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]. The quality of the fit justifies the use of approximate Eq. (25), which excludes for example possible direct interactions and cooperative effects between the molecules themselves.

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

5.2 Model

We follow [22], [59] 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 ω ₀ = 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 ₁ and u 2 = k × u 1 k . The Hamiltonian for the cavity modes is then

where a k z 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, ɛ ₀ 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, thus we consider a situation described by Eq. (2). The coupling between the molecules and the cavity modes is given by the Tavis–Cummings interaction

where e ± i k z z j 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 Section 4

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 we denote it by κ. The cavity damping is described in terms of additional zero mean white noise processes w t , k z 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 paper. Namely, we read off the terms of the NMQSD equation from the HTC Hamiltonian (Eq. (35)) and 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.

5.3 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 k z + μ c * a k z † 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. Our expression allows for the possibility that the classical electromagnetic field probing the system couples directly to the molecules as it is transmitted through the mirrors as well as indirectly after being first absorbed by the cavity mode. This is different than what is considered for example in [38], [42], where only cavity absorption is considered, which corresponds to the μ _m = 0 limit of our results. We use the refractive index to calculate the reflectivity of the cavity system by using TMM for incoming light with the angular frequency ω and z-component 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.

5.4 Reflectivity

We show the computed and measured reflectivity in Figure 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 Figure 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 Figure 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 Section 6. 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 Section 4. However, we set molecular disorder to zero and discuss this point in Section 6. To cover the full range of angles and energies used in the experiment, we use 21 cavity modes. Then the number of molecules is chosen to be high enough where additional increase does not change the reflectivity spectrum significantly. We choose N = 84. 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.15 for the refractive index that determines the cavity dispersion relation (30), E(0) = 3.41 eV for the energy of the cavity mode with k _z = 0, Ω = 0.91 eV for the Rabi splitting, and μ _c = 3.5μ _m for the magnitude of the cavity dipole moment.

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

Figure 4:

Error estimates for the polariton model. For small angles θ ≤ 30° the location and amplitude of the computed and measured upper and lower polariton peaks are within the estimated error bars. For larger angles the model is in good qualitative agreement with the measured reflectivity and predicts the locations of the polariton peaks.

With this model we can consider the classical field interacting with the system in three different ways. The first is that the light is absorbed by one of the cavity modes and there is no direct molecular absorption (μ _m = 0). This makes the upper polariton absorb more and more as the angle increases because it becomes increasingly photonic. The second is that the light can only be absorbed by the molecules and μ _c = 0. This results in both the upper and the lower polariton absorbing roughly the same amount. The third option is to include both cavity and molecular absorption which we found to match the experimental results the best. The vibrational coupling reduces the absorption of the upper polariton and this effect is greater at small angles where the polariton is mostly excitonic. Excluding this coupling 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 both the molecular and the cavity mode absorption leads to an agreement between the experiment and our theory.

5.5 Summary

We extend the thin film model from Section 4 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 [60], [61], where the coupling to the cavity mode diminishes the exciton broadening effect. There, however, the cavities had higher Q factors than 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 Figure 3. We consider the following deviations: The input energy has a standard deviation σ _E = 10⁻⁶ eV, the input angle σ _θ = 0.001°, and the thickness of the system σ _d = 10⁻⁹ m. We assume that the computed refractive index has energy dependent standard deviation σ n = 1 0 − 2 e − ( ω − ω 0 ) / 2 for the real and imaginary part, where ω ₀ is the smallest energy used. The deviation is smaller for higher energies, which is necessary for obtaining meaningful errorbars around the polariton energies. The variations in the input parameters lead to two kinds of errors: the amplitude of the polariton peak and the location of the polariton peak may be wrong. 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.

7.1 Fabrication

The samples were fabricated using thermal evaporation at a base pressure below 10⁻⁷ Torr (Angstrom Engineering physical vapor deposition system). We used 15 × 15 mm² 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 25-nm-thick aluminium layer as a top polariton microcavity mirror.

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