Tuning the Coherent Propagation of Organic Exciton‐Polaritons through Dark State Delocalization

Abstract While there have been numerous reports of long‐range polariton transport at room‐temperature in organic cavities, the spatiotemporal evolution of the propagation is scarcely reported, particularly in the initial coherent sub‐ps regime, where photon and exciton wavefunctions are inextricably mixed. Hence the detailed process of coherent organic exciton‐polariton transport and, in particular, the role of dark states has remained poorly understood. Here, femtosecond transient absorption microscopy is used to directly image coherent polariton motion in microcavities of varying quality factor. The transport is found to be well‐described by a model of band‐like propagation of an initially Gaussian distribution of exciton‐polaritons in real space. The velocity of the polaritons reaches values of ≈ 0.65 × 106 m s−1, substantially lower than expected from the polariton dispersion. Further, it is found that the velocity is proportional to the quality factor of the microcavity. This unexpected link between the quality‐factor and polariton velocity is suggested to be a result of varying admixing between delocalized dark and polariton states.


S1: Optical modelling of microcavities
We describe the optical properties of our microcavities using a transfer matrix model, the basis of which is described in detail in the "Basics of optics of multilayer systems". [1] We treat the microcavity as a series of isotropic, homogeneous layers, each defined by a unique complex refractive index n and thickness d. As input to the model, the absorption spectrum of the organic active layer was fitted using a series of Lorentzian functions, with the amplitude of the functions describing the excitonic oscillator strength. From fitting to the absorption spectrum, we determine the complex refractive index of the exciton absorption and superimpose this on the background refractive index of the polystyrene host. The incident beam is assumed to be propagating from the air towards the substrate at an angle of incidence θ, suffering multiple reflections/transmissions at each interface, interfering either constructively or destructively at the output, which can be either at the reflected side or the transmitted side. The model separately determines the characteristic matrix of each layer based on its material properties, and due to the continuity of the tangential components of the electric and magnetic fields at each interface we can determine the effective characteristic matrix of the total structure as the product of the single-layer characteristic matrices. This model was used to simulate the angle-dependent reflectivity and PL emission from the planar microcavity structure.
In order to minimize optical artefacts and oxygen-related degradation processes, the samples for TA microscopy were prepared on thin microscope coverslips and immediately encapsulated in inert atmosphere. This prevented them from fitting into our reflectivity goniometer to directly characterize the dispersion, so instead we describe their optical properties based on the measured properties of organic thin films made from the same solution batch and individual DBRs taken from the same deposition runs. The spectra of these components are illustrated in main-text Figure 1, and they can be well described by our model. We can thus fully model the reflectivity dispersions of our cavity Qfactor series, as illustrated in Figure S1 for an active layer thickness of 215 nm. The positions of the modulations detected in transient absorption, as reported in main-text Figure 1e, confirm the general picture suggested by this model. We note that penetration into the DBRs contributes to the cavity length in such structures; the invariance of the polariton dispersions through this series indicates that this penetration effect is the same for n=3.5 through n=6.5. To define the Q-factor of our structures, we turn to an equivalent optical model for each cavity with the dye absorption strength set to zero. The Q-factor is extracted from these "empty" cavity models as ⁄ at the mode minimum. The corresponding photonic lifetime quoted in the main text is determined through , where λ is the bare photonic mode wavelength at normal incidence. [2] Figure S1. Reflectivity dispersion as a function of cavity Q-factor. Transfer-matrix modelled reflectivity, benchmarked by the experimental data in main-text Figure  Aside from the expected changes in linewidth and total reflectivity, there are no systematic changes as a function of Q-factor. The Rabi splitting, detuning, and curvature of the UP and LP branches are fully invariant, suggesting there should be no alteration to transport properties of polaritons along our microcavity series.
Because the solution-processed thin films used as an active layer in our structures inevitably exhibit thickness variation, [3,4] we use the transfer matrix model to establish what range of cavity thicknesses would be consistent with the chief experimental constraint of our TAM measurements, namely the position of the derivative-like feature at 640 nm. The sensitivity of the polariton dispersions to the cavity thickness means this condition strongly constrains the cavity thickness. Our model results indicate this behaviour is consistent with an active-layer thickness of 200-215 nm; as shown in Figure  S2, this range entails a roughly 10-nm shift in the LP minimum. We use the thickness range in Figure S2 as the basis to determine the expected group velocities in our samples. We first use a simple two-state coupled oscillator model to describe the dispersions predicted in our transfer matrix approach. We determine the probability of our wide-angle, broadband excitation source to photoexcite a particular polariton state based on the product of the angle-dependent photonic Hopfield coefficient squared and the intensity of the pump pulse at the corresponding wavelength. The range of k-vectors considered spanned from 0 to a k value corresponding to the intersection of a free photon at 64 degrees (the limit of our high-NA objective) with the UP or LP dispersion. The integration was also limited by the wavelength range spanned by our excitation pulse (550-630 nm), which reduced the degree of excitation of the lower portion of the LP dispersion. This approach yields polariton population distributions along the LP and UP; these are typically heavily weighted towards the bottom of the LP and the top the UP, where the states are most photonic. At each k-vector we determine the expected group velocity and use these population distributions to produce a weighted average velocity for UP versus LP excitation; we keep these separate because the precise contribution of the (much more rapidly moving) UP on the short timescales of our measurement is not well determined. The process was run for structures with 3.5, 4.5, 5.5, and 6.5 mirror pairs. The same procedure was repeated for cavity active layer thicknesses 200-215 nm, which were averaged together at each Q-factor to produce the mean and standard deviation results in main-text Figure 2e. Within the limits of the transfer matrix and coupled oscillator models, there is no reason to anticipate any impact of Q-factor on polariton group velocity. Figure S3. Rise time of LP PIA band as a function of Q-factor. The PIA band slightly redshifted from the steady-state LP position exhibits increasing rise time with increasing Q-factor, in a similar manner to the bleach peaks discussed in Figure 1 of the main text. This strong Q-factor dependence confirms assignment to a polaritonic species. Figure S4. Decay dynamics of R-BODIPY films embedded in microcavities. a-b. In general, the decay time of the signals observed at the spectral position of the upper polariton branch (UP, based on steady-state measurements) increases with an increasing number of SiO 2 -Nb 2 O 5 mirror pairs (n). This is further reflected in the increased rise and decay times for the signals at the lower polariton branch spectral position (LP, again referred to the steady-state assignment) with an increasing n. We note that on longer timescales (10 ps -500 ps) a transient signal persists but its rate of decay does not show a Q-factor dependence, i.e. the decay curves have a constant gradient after the log-break in a and b. c. Decay dynamics for n = 3.5 mirror pairs microcavity as a function of the intensity of the exciting laser pulse. The dynamics appear independent of the laser fluence in the range applied. d. Decay dynamics of bare (~150 nm thick) R-BODIPY film as a function of laser fluence, revealing clear power dependence.

S4: Extraction of
We can define as the distance between the centre of the signal and the point in the image plane (in any direction) where the signal strength is of that at the centre. Mathematically this is expressed by where and are the coordinates of the signal centre and the points where the signal strength at time drops to of that at the centre at angle relative to the positive xaxis, respectively. We will refer to as the "sigma points" from now on. Extracting from the experimental data allows monitoring of the spatial movement of population as a function of time.
To extract from our experimental data we must first extract the centre of our signal. It is tempting to use the pixel with the strongest signal in each frame as the centre, as one might expect its position to be constant over time. However, since one pixel on the EMCCD camera corresponds to 55.5 × 55.5 nm 2 in real space and the pump spot has a FWHM of ~270 nm, the central four pixels often receive similar numbers of photons. In practice, the position of the maximum signal pixel fluctuates between 1-2 pixels in either direction. This is undesirable as physically the centre of the signal is defined as the centre of the pump pulse, which has a well-defined position that stays constant throughout the measurement window. To overcome this we utilise a number of pixels around the maximum to estimate the centre through a weighted average. It works almost identically to calculating the centre of mass (CoM) of discrete masses on a plane, which can be expressed algebraically using the following equations: is the coordinates of the centre of mass, is the total mass, and are the mass and coordinates of the i-th individual mass respectively. To calculate the centre of our signal, we substitute the ΔT/T data points for the individual masses in a 10 × 10 square window surrounding the maximum signal pixel (i.e. ). We can then average the CoM over several frames and calculate its uncertainty in a similar fashion.
After we have estimated the signal centre, we need to extract the sigma points (i.e. coordinates where the signal strength is of that at the centre), at any angle and time with sub-pixel resolution. This is achieved by mapping the data onto a fine polar grid centred at followed by selecting in each direction the point which has an amplitude closest to of that at the centre. can then be calculated.
Finally, we need to define time zero (i.e. the moment when the pump pulse arrives at the sample), for which three common definitions exist: Here, t 0 is placed when the ΔT/T signal is at the maximum. This is where is found to be closest to the diffraction-limited pump spot size ( and avoids any potential artefacts (coherent artefact, pump scatter, etc) from the arrival of the pump pulse. S5: Pump-probe spot size at t 0 Figure S5. fs-TAM of bare BODIPY-R films and half-cavities. a. fs-TAM decay kinetics from bare BODIPY-R films and half-cavity (BODIPY-R film on a 3.5-pair SiO 2 /Nb 2 O 5 stack). The lifetime is taken at the GSB of the signal ~640 nm. The MSD traces for both remain flat over the first 1.5 ps confirming no spatial transport occurs in this regime in the absence of polaritons. b. Spot size at t 0 (σ) increases with the number mirror pairs due to refractive index mismatch between microscope oil and cavity which increases with n (see below). This systematic artefact is directly accounted for in our analysis and does not affect the conclusions drawn. Figure S6. Variation of excitation spot size with mirror thickness. TAM images of the tightest observable reflected pump spot from microcavities of varying numbers of mirror pairs n. The size of the reflected spot increases with n, due to the refractive index effects discussed above and in the main text. Scale bar in each panel is 500 nm. Figure S7. Independent characterization of excitation spot size. Photoluminesence intensity profile obtained on scanning pump pulse over a fluorescent bead (200 nm) on half-cavities of varying number of mirror pairs n. Blue shows raw data and orange a Gaussian fit. By deconvolving the spatial profile from the known diameter of the bead, the size of the pump spot can be obtained. The size of the imaged pump spot (varying between 130 ± 10 nm and 190 ± 10 nm in σ (see above)), though nominally fixed by our objective N.A. to a diffraction limit of ~135 nm in σ, varies precisely as described above, confirming the validity of our assignment of the polariton population size at t 0 in each cavity.

S6: Modelling coherent polariton propagation
At long times an initial Gaussian population will follow diffusive dynamics due to multiple scattering events effectively randomising the velocities of the carriers and giving the overall profile the dynamics of a Markovian random walk. It can be shown in this t >> τ limit that the characteristic diffusion constant in a 1D system appears as, It is desirable to have a similar σ(v,t) for the coherent transport regime in order to reliably extract the coherent group velocities v from the σ(t). This is derived below.
We assume our microcavity has the simplest possible 1D band structure -a two band model for the ground state and the polariton branch. We further assume that we are exciting with a delta function like excitation. We assume the excitonic bands are sufficiently flat that the overall dispersion obeys inversion symmetry and the two k states we are exciting are at exactly +k" and -k". As we are exciting states of a determinate k, they have a welldefined group velocity v g .
Our initial excitation is a Gaussian with a known σ 0 with half the excitation having a velocity +v g and the other half a velocity − v g . Hence we can decompose our initial excitation"s time evolution as, Examining the above equation it is conceivable that in the limit v g t/σ 0 << 1 the profile still appears Gaussian and we still have an effective σ. As coherent transport velocities are typically < 10 6 m s -1 = 1 nm fs -1 , σ 0 ≈ 100 nm for a diffraction limited spot and the times range in which this physics is relevant is ≤ 50 fs at room temperature. For these parameters v g t/σ 0 < 1 is well justified. Experimentally it is also well known that in this regime our profiles are Gaussians and we can find σ(t). This expression therefore requires a Taylor expansion in .
This expansion"s terms are zero for odd b, hence setting b = 2b we collect the relevant terms Knowing that is small we expand to the lowest relevant order. It is also important to note here that a,2b,c Defining a Gaussian that is not normalised as, we can write, A is simply the unperturbed Gaussian sitting at µ = 0 which defines the overall TAM profile as Gaussian even in the coherent transport regime. It is expected that we would have such a term as we are expanding for small changes around this profile. B is an inverted Gaussian that is responsible for reducing the maximum value of our Gaussian by a small amount. This results in a "dip" in the centre of the profile that we usually cannot resolve. C is what is of interest to us as it represents a Gaussian that is being "stretched out" by a term in x 2 . One naturally expects C to capture to lowest order in vt/σ 0 the changes in σ of the overall profile.
We begin to try and find an effective σ for this distribution by starting with the ansatz, We set µ = 0 for simplicity and write σ(t) = σ 0 + f(v,t), where f(v,t) will capture the effective small changes in σ due to coherent transport. We can naturally also make the approximation f(v,t)/σ 0 << 1 that results from the approximation v g t/σ 0 < 1.
Comparing the expansion of the n ansatz with n, it is clear that the term A and D are the same Gaussian by setting A 0 = 2A as expected. It is not surprising that we were not able to capture the term B with the Ansatz as it is a single Gaussian and cannot have any "dips" in the middle caused by the term B -we also note that this term is a Gaussian of with a maximum amplitude lower than the main Gaussian. Finally and most importantly, comparing terms C and E we see clearly that they both are of the form differing by a factor of 2. Absorbing the factor of 2 into the gaussian and setting these terms as equal we get,