Superabsorption in an organic microcavity: Toward a quantum battery

In a major step toward the development of a quantum battery, superabsorption has been achieved in an organic microcavity.


Introduction
The properties of physical systems can typically be categorised as intensive (i.e. they are independent of the system size, such as density) or extensive (i.e. they grow in proportion to system size, such as mass).However, in some cases, co-operative effects can lead to super-extensive scaling.A well-studied example of this is superradiant emission (1).In its original form, this describes emission from an ensemble of N emitters into free space.Constructive interference in the emission process means that the time for emission scales as 1/N, so that peak emission power is superextensive, scaling as  2 .
Such behaviour has been demonstrated on a number of platforms (low pressure gases (2,3), quantum wells (4,5) and dots (6), J aggregates (7), Bose-Einstein condensates (8), trapped atoms (9), nitrogen-vacancy centres (10)).A less-studied example is superabsorption (11), describing the N-dependent enhancement of absorption of radiation by an ensemble of N two-level systems (TLSs).Only very recently has this been demonstrated for a small number of atoms (12).In principle, superabsorption could have important implications for energy storage and capture technologies, particularly if realised in platforms compatible with energy harvesting, such as organic photovoltaic devices.However, there are challenges in engineering the precise environment in which such behaviour can occur, and in monitoring the ultrashort charging time scales.In this article we show how these can be overcome, by combining organic microcavity fabrication with ultrafast pump-probe spectroscopy.
Superextensive scaling of energy absorption is also a key property of quantum batteries (QB).These represent a new class of energy storage devices that operate on distinctly quantum mechanical principles.In particular, they are driven either by quantum entanglement, that reduces the number of traversed states in the Hilbert space compared to (classical) separable states alone (13)(14)(15)(16)(17)(18)(19)(20)(21), or by cooperative behaviour that increases the effective quantum coupling between battery and source (22)(23)(24).Such effects mean that QBs exhibit a charging time that is inversely related to the battery capacity.This leads to the intriguing idea that the charging power of QBs is superextensive; that is, it increases faster than the size of the battery.For a QB consisting of a collection of  identical quantum systems, superextensive charging rate density (charging rate per subsystem) that scales as  or √ in the thermodynamic limit (20), has been predicted.
Here we experimentally realise a paradigmatic model proposed as a Dicke QB (24), which displays superextensive scaling of energy absorption, using an organic semiconductor as an ensemble of TLSs coupled to a confined optical mode in a microcavity.We also demonstrate how dissipation plays a crucial role; in a closed system, the coherent effects that lead to fast charging can also lead to subsequent fast discharging.As such, stabilisation of stored energy remains an open question: proposed stabilisation methods include continuous measurements (25), dark states (21), and novel energy trapping mechanisms (26,27).In our open noisy system, dephasing causes transitions between the optically active bright mode, and inactive dark modes.This suppresses emission into the cavity mode, so that we have fast absorption of energy but slow decay, allowing retention of the stored energy until it can be used.

Device structure
The structures fabricated consist of a thin (active) layer of a low-mass molecular semiconductor dispersed into a polymer matrix that is deposited by spin-coating and positioned between two dielectric mirrors, forming a microcavity as illustrated schematically in Fig. 1a (see Methods section for fabrication details).Organic semiconductors are particularly promising for many applications as the high oscillator strength and binding energy of molecular excitons means that light can be absorbed efficiently and excitons can exist at room temperature (28).The organic semiconductor used in this study was the dye Lumogen-F Orange (LFO), whose chemical structure is shown in Fig. 1b.The normalised absorption and photoluminescence spectra for LFO dispersed at 1% concentration by mass in a polystyrene (PS) matrix are shown in Fig. 1b.
By diluting the LFO, we reduce intermolecular interactions that lead to emission quenching, producing a high photoluminescence quantum yield of around 60% at low concentration (see Supplemental Fig. S1).The absorption peak at 526 nm and the emission peak at 534 nm correspond to the 0-0 transition, i.e. an electronic transition from and to the lowest vibrational state.Operating around the 0-0 transition, the LFO molecules can reasonably be considered as a TLS.We prepared samples with 0.5%, 1%, 5%, and 10% concentrations, as these are representative of the optimal operating regimes -further increases in concentration lead to quenching, and signals from lower concentrations are indiscernible from noise.The absorption and photoluminescence spectra for the 0.5%, 5% and 10% concentrations are given in Supplemental Fig. S2.The optical microcavities fabricated support cavity modes whose energy is determined by the optical-thickness of the LFO layer and the penetration of the optical field into the cavity mirrors (29).The confined photon field drives coherent interactions with the molecules, which underpin the collective effects that drive superabsorption.The LFO concentration dictates the operating coupling regime, with the 0.5% and 1% LFO cavities operating in the weak coupling regime, the 5% in the intermediate coupling regime, and the 10% in the strong coupling regime (see Supplemental Fig. S2 and discussion in Methods).

Experimental setup
Charging and energy storage dynamics were measured using ultrafast transientabsorption spectroscopy (30), allowing femtosecond charging times to be measured.In this technique, we excite the microcavity with a pump pulse, and then measure the evolution of stored energy (i.e.corresponding to the number of excited molecules) with a second probe pulse, delayed by time  (Fig. 1d).The probe pulse is transmitted through the top distributed Bragg reflector (DBR) of the cavity, and the reflection from the bottom DBR is measured.The differential reflectivity induced by the pump-pulse is given by where   (  ) is the probe reflectivity with (without) the pump excitation.Note that control films (active layers without the microcavities) are measured under differential transmittivity Δ ∕ .The control films will allow us to identify the underlying photophysics of the molecules.
In our experimental setup (shown schematically in Fig. 1e), transient-absorption measurements were performed in a degenerate, almost collinear configuration.Pump and probe pulses were generated by a broadband non-collinear optical parametric amplifier (NOPA) (31) and spanned the wavelength range 500 to 620 nm with a nearly transformlimited sub-20-fs duration (further details in Methods).An optical delay line was used to control the probe delay time, and a mechanical chopper was used to modulate the pump pulse, providing alternating probe-only and pump-probe pulses, allowing us to measure pump-induced absorption changes.Measurements at different molecular concentrations were performed, adjusting the pump fluence in order to maintain an approximately constant photon density (i.e.pump photons per LFO molecule)  =    ⁄ , where N is the total number of molecules in the excitation volume,   is the total number of pump laser photons, and  is the fraction of them that actually reach the active layer of the microcavity.We estimate from the reflectivity data that only 6% to 8% of the initial pump excitation enters the cavity.We conducted our experiment in air at room temperature.

Results
We first show that ultrafast transient-absorption spectroscopy can monitor the population of excited molecules, even in a cavity, by comparing the control film and the microcavity spectra as shown in Fig. 2a is challenging, so we instead extract it through fitting to the theoretical model which is discussed below.This fitting is discussed in section S3 of the Supplementary Material.We also note that two of the microcavity spectra show a negative Δ  ⁄ band, which results from the change in the refractive index induced by the pump pulse (33).
Figure 2b shows the experimental values for the time-dependent stored energy density.In all microcavities studied, the energy density undergoes a rapid rise followed by slow decay.The timescale of the rapid rise varies with concentration.We adjust the laser power to fix photon density r across comparable microcavities, and compare behaviour with different LFO concentrations.Details of how r is estimated are provided in the supplementary information.We found that to achieve a sufficiently high signal-tonoise ratio, it was not possible to compare all microcavities at the same r value; instead, a constant r value was maintained for matched structures.Specifically, measurements were made on microcavities with LFO concentrations of 10%, 5% and 1% with approximately constant  ≃ 0.14 (respectively labelled A1, A2, and A3), and 1% and 0.5% with  ≃ 2.4 (labelled B1 and B2).Overlaying the experimental data are the corresponding theoretical predictions (see Theoretical Model section).To account for both the instrument response time (~20 fs) and the cavity photon lifetime (which was estimated as discussed in the supplementary information), the theoretical curves are convolved with a Gaussian response function with a full-width half maximum of ~120 fs.There is good agreement between the experimental data and the corresponding theoretical curves.
To obtain energetic dynamics, we take away the response function from the theoretical fit, as shown in Supplementary Figure S15.

Theoretical Model
The experimental dynamics can be reproduced by modelling, with the Lindblad master equation, the N TLSs in an optical cavity with light-matter coupling strength g, a driving laser with a Gaussian pulse envelope and peak amplitude  0 , and three decay channels corresponding to the cavity decay (), TLS dephasing (  ), and TLS relaxation ( − ).To solve this many-body Lindblad master equation, we make use of the cumulant expansion (34)(35)(36), with model parameters given by a chi-squared minimisation of the experimental data.Experimental uncertainties are estimated from the point-to-point variance of the data.Further details can be found in the Methods and Supplementary Information.
From our cumulant expansion simulations, we show how ,  max , and  max vary as a function of  in Fig. 3a and b.The interplay amongst the decay channels, driving laser, and cavity couplings give rise to a rich set of behaviour.We identify three regimes: decay-dominated at small N, and coupling-dominated at large N, along with a crossover regime between them.The system exhibits superextensive energy density scaling in the decay-dominated regime, and subextensive charging time in the coupling-dominated regime.In the crossover regime, the system exhibits both superextensive energy density scaling and subextensive charging times.Importantly, charging power density is superextensive in all regimes.
Figures 3c and d show the typical time dependence in decay-dominated and coupling-dominated regimes, indicating how the model parameters affect the dynamics.
In particular, the presence of the decay channels gives rise to ratchet states which are capable of absorbing but not emitting light (37), thereby allowing the energy to be stably stored.See Methods and Supplementary Information for further discussion on the operating regimes.Figure 3 is augmented with an animation of how the energetic dynamics changes with N (see Supplementary Video).

Figures 3a and b provide an explanation for the different scaling of the observables
with N shown in Table 1.Specifically, A1 and A2 operate in the coupling-dominated regime, where  scales slightly less than  −1 2 ⁄ ,  max scales slightly more than  0 , and  max scales slightly more than  1 2  ⁄ .For the region between A2 and A3, the average scaling of  falls between  0 and  −1 2 ⁄ ,  max between  2 and  0 , and  max between  2 and  1 2 ⁄ .As A2 is further in the coupling-dominated regime than A3 is in the decaydominated regime, the average scaling values between A2 and A3 are skewed towards the coupling-dominated scalings.B1 and B2 operate in the cross-over regime, with an average scaling with  that is between the decay-dominated and coupling-dominated scalings, as reflected in Table 1.

Summary
We have provided direct experimental evidence of superextensive energy storage capacity and charging in an organic microcavity by using ultrafast optical spectroscopy.
Our realisation of a prototype Dicke QB highlights the fact that purely closed unitary dynamics is insufficient for realising a practical QB.The retention of energy requires finely-tuned decoherence processes, allowing the battery to charge quickly and yet discharge much more slowly.Such stabilisation of stored energy is a key step to exploit superextensive charging.Our observation of dephasing shows that realistic noisy environments can aid implementation and application of useful QBs.A challenge for future work is to explore further how concepts of ratchet states could keep a QB operating in the range of higher-lying energy states that are associated with maximum absorption enhancement, i.e. near the mid-point of the Dicke ladder (37).
We conclude by discussing the potential for future applications based on superextensive charging.One practical challenge noted above, is that quenching limits the performance of the QB at high concentrations.Overcoming this limitation requires careful choice of materials in order to suppress intermolecular quenching.We note that there are classes of materials where quenching is particularly suppressed.For example in proteins such as Green Fluorescent Protein (GFP) (38) the active chromophore is surrounded by a cage, which suppresses exciton-exciton quenching at high intensities.
Such materials might provide a route to allow the study of higher concentrations.Beyond energy storage, the key challenge for practical applications of this effect is its integration in devices where energy can be efficiently extracted and used.While our focus has been on the quantum advantage in charging, there do exist approaches to efficiently extract energy.For example, this may be achieved by including charge transport layers between the active layer and the cavity layers (39).The transport layers allow charge separation of the excitons as well as preventing recombination.This transforms the top cavity layer into a cathode and the bottom cavity into an anode, giving rise to an electric current.As such, our work provides a direct path for the integration of the superextensive energy absorption process in an organic photovoltaic device.The fast dynamics of such a device may also be useful as an optical sensor in low-light conditions, or potentially for energy harvesting applications (40)(41)(42)(43).More generally, the idea of superextensive charging may have widereaching consequences for sensing and energy capture and storage technologies.

Device fabrication
The microcavities constructed consist of a thin layer of LFO (Kremer Pigmente) dispersed in a PS (Sigma-Aldrich, average molecular weight ~ 192,000) matrix.The bottom DBR consisted of 10 pairs of SiO2/Nb2O5 and were fabricated using a mixture of thermal evaporation and ion-assisted electron beam deposition by Helia Photonics Ltd.Solutions of LFO dissolved in 25 mg/mL PS in dichloromethane were prepared at 0.5%, 1%, 5%, and 10% concentration by mass.Each LFO solution was then spin-coated on top of the bottom DBR to produce a thin film with an approximate thickness of 185 nm.An 8-pair DBR was then deposited on top of the LFO layer using electron beam deposition.With this pair of mirrors, the reflectivity was >99% in the spectral region of interest (44).
The diluted molecules are expected to be isolated at low concentration 0.1 -1%, but at higher dye concentrations, the 0-0 emission transition red-shifts by a few nm and the second peak increases in intensity due to aggregation of the dye molecules.This is evident in Supplemental Fig. S2a and b, with additional broader features observed at longer wavelengths, which we assign to intermolecular states such as excimers.
The 0.5% and 1% cavities lie in the weak-coupling regime, i.e. no polaritonic splitting could be seen in the cavity reflectivity spectrum, as shown in Supplemental Fig. Supplementary Fig. S3 shows a transfer matrix simulation of the electric field distribution of the 1% cavity (the cavities exhibit similar distributions).

Pump-probe spectroscopy
Probe and pump pulses were generated by a broadband non-collinear optical parametric amplifier (NOPA).The NOPA was pumped by a fraction (450 μJ) of the laser beam generated by a regeneratively amplified Ti:Sapphire laser (Coherent Libra) producing 100 fs pulses at 800 nm at a repetition rate of 1 kHz.A pair of chirped mirrors were placed at the output of the NOPA to compensate for temporal dispersion, and by using 7 'bounces' we were able to generate pulses with a temporal width below 20 fs.The laser beam was then split by a beam-splitter, with the probe being delayed via a translation stage and the pump being modulated mechanically using a chopper at 500 Hz.

Lindblad master equation
As noted above, we find the experimental behaviour is well reproduced by the dynamics of the Dicke model, a model of a microcavity photon mode coupled to two-level systems representing the molecules.As further discussed in the supplement, such a model is generally an approximation for organic molecules, but for some systems can become a very accurate approximation in the limit of low temperatures (45).
The open driven nature of the experimental system is modelled with the Lindblad master equation, where () is the density matrix and , and a carrier frequency ω  .We work in the frame of the laser carrier frequency, and so write where Δ = ω − ω  is the detuning of the cavity frequency from the laser driving frequency.The LFO molecules are initially in the ground state, and the laser is onresonance (Δ = 0).

Cumulant expansion
The energy density of the cavity containing identical molecules with transition energy depends on both the first order moments ⟨ ,, ⟩ and ⟨⟩ as well as higher order moments, leading to a hierarchy of coupled equations.Within mean field theory, the second order moments are factorised as ⟨⟩ = ⟨⟩⟨⟩ which closes the set of equations at first order.
This approximation is valid at large N, as corrections scale as 1  ⁄ .To capture the leading order effects of finite-sizes we make a second-order cumulant expansion (34)(35)(36), i.e. we keep second-order cumulants ⟨⟨⟩⟩ = ⟨⟩ − ⟨⟩⟨⟩ and assume that the third-order cumulants vanish, which allows us to rewrite third-order moments into products of first and second-order moments (46).In our experiments, the number of molecules in the cavity is large (>10 10 ) and we find higher order correlations are negligible.We give the equations of motion up to second order in the Supplementary Information.

Operating regimes
The decay-dominated (purple region in Fig. 3a and b) regime occurs when the collective light-matter coupling is weaker than the decay channels, √′ < {,   ,  − }, where ′ = max(1, ).In this regime, the time scale of cavity dynamics is slow relative to the decay rate.Fig. 3c shows a typical time dependence in this regime, indicating how the model parameters affect the dynamics.In this regime, the increase in the effective coupling relative to the decay strength sees an  2 superextensive scaling of the energy and power density, while rise time remains constant.Experiment A3 operates near the boundary of this regime (Fig. 3a).
In the coupling-dominated (green region in Fig. 3a and b) regime, the effective collective light-matter coupling √′ > {  ,  − ,  }, dominates over the decay channels.In this regime, the time scale of cavity dynamics is fast relative to the decay rate, and we observe √-superextensive power scaling and 1 √ ⁄ dependence of rise time, while the maximum energy density remains constant.While power scaling is superextensive in both regimes, the origin of this differs: for the decay-dominated regime this is the result of the superextensive energy scaling, while for the coupling-dominant regime it is the result of a superextensive decrease in the rise time.Experiments A1 and A2 operates in this regime (Fig. 3a).
In the crossover between the regimes (purple-green), the collective coupling falls between the cavity decay rate and the TLS dephasing rate, {,  − } < √′ <   .In Fig. 3a and b,  − is small such that √′ ≫  − for all values of N, and so there is no boundary labelled for this decay rate.In this case, capacity and rise-time can simultaneously scale super-and subextensively, but at a rate slower than in the decay and coupling-dominated regimes, respectively.Experiments B1 and B2 operate in this regime (Fig. 3b).

Decay and coupling rates
The parameters needed in the theory calculations are the cavity leakage rate , the dephasing rate   , the non-radiative decay rate  − , the interaction strength , and the temporal width of the instrument response function,   .Note that the temporal width of the pump pulse is fixed at  = 20 fs.For the dephasing rate, we note that as one enters the strong-coupling regime, exciton delocalisation suppresses the effect of dephasing (47).
To approximately capture this effect, we assume that the dephasing rate scales with the number of molecules as   =  0  (  5%  ⁄ ) where  0  is taken to be constant and  5% is the number of molecules in the 5% cavity.The experimental uncertainty in N is estimated to be 10%.The cavity lifetime  comes into the model in both   =  and the cavity leakage rate  = 1/.From transfer matrix modelling on the 1% and 0.5% cavities (where polariton effects are small) we estimate that  ≈ 306 fs.However, based on the measured finesse of the cavities we estimate that  = 120 fs.Transfer matrix modelling assumes perfectly smooth mirrors, whilst measured finesse includes inhomogeneous broadening effects, neither of which we want to include in  and   .In the following optimisation, we therefore assume that  ∈ [120, 306] fs, with lower values more likely due to transfer matrix calculations being prone to error.
For  values within this range, the remaining three parameters in the model This section presents further details of the properties of the fabricated microcavities, and measurements used to calibrate the results in the main text.
a. Quenching at large concentrations Figure S1 shows the photoluminescence quantum yield (PLQY) as a function of dye concentration.At large concentration the yield drops to zero; this provides an upper limit on the concentration that can be studied in experiment.The photoluminescence measurements were taken using a Coherent Mira 900 laser operating at 400 nm with a repetition rate of 80 MHz.The laser beam was focused onto the surface of a sample placed at the centre of an integrating sphere.The laser light and the emission from the sample were scattered by the diffuse interior of the sphere and collected by an optic fibre, which was coupled to an Andor Shamrock SR-303i-A CCD spectrometer.Spectra were taken for different concentration LFO films, as well as a blank glass substrate in order to calculate the proportion of laser light absorbed by the LFO samples.

c. Transfer matrix calculations
To provide bounds on the cavity lifetime, separate from the measured cavity linewidth-which contains effects of inhomogeneous broadening-we make use of transfer matrix calculations of the cavity structure.These calculations, shown in Fig. S3, also enable one to visualise the electric field intensity in the microcavity structure.These calculations give a designed cavity lifetime of 306fs, which serves as an upper bound of the actual cavity lifetime.Figure S4 compares the dynamics at 525 nm (ground state bleaching) and 571 nm (stimulated emission region) for the 1% concentration film.Unfortunately a coherent artefact, due to the degenerate pump-probe configuration, masks the time dynamics in the first 100 fs.However, the same rise and decay times are seen at both probe wavelengths, indicating that we are probing the same exciton population.
We also explored pump-fluence dependence of the bare films.No dependence on the excitation fluence was detected in any control film.This indicates the absence of bimolecular effects or multi-photon excitation.Figure S5 shows such data for all concentrations, at two different pump fluences.As in Fig. S4, a coherent artefact is present at zero delay due to the degenerate pump and probe spectra.e. Estimating number of molecules in film.To determine the number of LFO molecules in the microcavities we study, we first determined the absorption cross-section of a single LFO molecule, σ LF O .The transmission spectrum of a 0.1% solution of LFO in 25 mg/mL PS/dichloromethane in a 1 mm thick cuvette was measured using a Horiba Fluoromax 4 fluorometer with a xenon lamp.The absorption coefficient (α = nσ LF O ) of the 0-0 transition was then calculated using the relation T /T 0 = e −αd , where T /T 0 is the fractional transmission of the xenon lamp at the 0-0 transition, d is the cuvette thickness, n is the number density of absorbing molecules in solution per unit volume, and σ LF O is the absorption cross-section of a sin-LFO molecule [49].Using the known value of n for this solution, σ LF O was calculated as 3.3 × 10 −16 cm 2 .The transmission of the 10% LFO concentration in film was then measured to obtain α and hence n (number density of molecules in the cavity active layer), using the measured value of σ LF O , with d (film thickness) measured using a Bruker DektakXT profilometer.This value was then multiplied by the area of the laser beam and d to obtain N .Here we assume a uniform distribution in the active layer.N for other concentrations were scaled accordingly.
f. Estimating number of photons in cavity.To estimate the number of photons entering the cavity in each different cavity, we consider the overlap between the pump spectrum and the cavity transmission.The number of photons in the cavity is given by multiplying the number of pump photons by 1 − R, where R is the reflectivity of the cavity: n = N γ (1 − R).An example of this is shown in Fig. S6, for the 1% cavity.Starting from the pump spectrum (yellow curve), by considering the reflectivity spectrum of the cavity (black line), we calculate the fraction of photons entering into the cavity (purple line).Table S1 shows the resulting estimates of photon numbers for each experiment.

Experiment
N dye (×10 10 ) N photon (×10 g.Derivative features in the differential transmitivity As seen in Fig. 2 of the main text, the transient signal shows both positive and negative features in the differential transmitivity.This can be understood from the existence of a derivative feature in the spectrum.Such a feature occurs if the pump causes an absorption peak to move in energy.In that case the change of absorption will see a decrease in absorption where the feature used to be, and an increase where the feature now is [33,50].This leads to a feature that takes the form of the derivative of the original absorption peak with respect to energy, and which thus contains both negative and positive contributions.This structure is exactly what we observe, indicating that the feature corresponds to a resonance which moves with excited state population.This feature though coexists with other positive features in the transient reflectivity, arising from mechanisms such as ground state bleaching or stimulated emission.Whether a negative feature is seen depends on how these contributions compete and whether the feature is above the detection threshold. The presence of this derivative feature does not influence our determination of the energy stored in the molecules and the charging time.As shown in the Fig. S7, the time dependence of the negative and positive peaks is the same.(The negative feature does show more noise, as this signal is weaker compared to background noise.)Since both features show the same time evolution, they would both recover the same theoretical fits.

S2. THEORETICAL MODELLING
In this section we describe our approach to modelling the system.We first provide a more detailed discussion of the model we use, and the physical origin of the terms involved.We then derive the equations of motion for the expectation values using a second-order cumulant approach.

A. Details of theoretical model
As discussed in the main text, we model the experiments through a Dicke model-which describes N two-level systems coupled to a photon mode.Describing organic molecules as two-level systems is an approximation, as it neglects both the role of rotational and vibrational molecular modes, as well as the existence of higher excited states of the molecules.It is known that this approximation can become valid in some limiting situations, such as at low temperatures [45].While our experiments are performed at room temperature, we nonetheless expect the model to provide a reasonable approximation.This is because, as noted in the main text, the molecules we consider show a small Stokes shift, indicating vibrational dressing of optical transitions is weak.Moreover, as seen in the main text and discussed further below, our model matches the experimental results well without such additional features.
The Hamiltonian describing our system, including the external pump laser, takes the form (setting = 1): where a (a † ) is the photon annihilation (creation) operator, σ α i for α = x, y, z are the Pauli matrices.We write the Hamiltonian in the rotating frame of the pump laser, so ∆ a (∆ c ) is the energy detuning of the laser from the cavity (molecules), g is the coupling strength of each molecule to the photon mode, η(t) is the Gaussian profile of the pump laser.In the following we assume that the molecules and cavity are resonant (∆ a = ∆ c ≡ ∆) and that ∆ = 0 unless specifically noted.
To account for dissipative processes, we consider the time evolution of the density matrix, including Lindblad terms for various dissipative processes: where The first term, with rate κ, describes loss of photons, due to the finite reflectivity of the cavity mirrors.The second term, with rate γ z , describes dephasing of the molecular excitations.This term describes the effect of coupling between the molecular electronic state and vibrational degrees of freedom-vibrations both of the molecule and of the polystyrene matrix in which it is contained.The final term, γ − , describes the decay of electronic excitation, due to emission into non-cavity modes, such that 1/γ − would be the excited state lifetime in the absence of the cavity.
As noted in the main text, dephasing γ z plays a crucial role in the dynamics, with quite distinct behaviour occurring with and without dephasing.In particular, dephasing introduces transitions between the "bright" and "dark" molecular excited states.To understand these states, let us first consider states with a single excitation.The form of Eq. (1) shows that the cavity photon couples only to the totally symmetric molecular excited state, i.e. a state with equal weight and phase of excitation on all molecules.However, for N molecules, N excited states exist.The remaining N − 1 states are orthogonal to the coupling to light, and are known as "dark states".While we have introduced this for the space with a single excitation, a generalization to higher excited states exists.In this case, the language of superradiant and subradiant states is often used [1], with superradiant states referring to those that can be reached using the collective raising and lowering operators, j σ ± j .Because the dephasing term acts on individual molecules, it describes a process that allows loss of phase coherence between different molecules.As such, this causes a transition from the optically bright state created by the laser, to an incoherent mixture of bright and dark states.Since the dark states do not couple to the cavity, this process is responsible for the asymmetry between collectively enhanced absorption, and the lack of collective enhancement of emission.
For N identical molecules of energy ω a (in the lab frame), the energy density stored on the molecules is given by As such, in the following, our aim is to predict the time evolution of this quantity.

B. Cumulant equations
To determine the time evolution of σ z (t) = Tr [σ z ρ(t)] we begin by writing down the first order expectation values of the system.We adopt the notation C a (t) ≡ a(t) for photon operators, and C α=x,y,z (t) ≡ σ α (t) for spin operators, along with a similar notation for higher order expectations, e.g., C ax (t) ≡ aσ x (t) .The equations of motion for the first order expectation values are where ∂ t is short for ∂ ∂t , γ tot = 2γ z + 1 2 γ − , and for notational ease we have dropped the explicit time dependence of observables.As described in the main text, in mean field theory we would now set the second order cumulants to zero.These are defined as This would result in the usual decomposition of second order expectation values into products of first order ones, C AB = C A C B , which is the assumption that molecule-molecule, molecule-photon, photon-photon and all higher order correlations are negligible.However, we instead derive equations of motion for the second order expectation values, capturing the leading order 1/N corrections to mean field theory.The second order photon correlations obey: while molecule-photon correlations follow: These now depend on third order expectation values, some of which contain multiple Pauli operators.We must note that these terms indicate Pauli operators representing different molecules and so commute -we have already taken into account the cases where the Pauli operators correspond to the same molecule by using the Pauli algebra σ α σ β = 1δ αβ + iσ γ αβγ .The molecule-molecule expectation values for the same Pauli operator acting on different molecules are Finally, the molecule-molecule expectation values for different Pauli operators acting on different molecules are In principle one can continue to write equations of motion for increasingly higher orders of expectation values, however, at large N , most essential physics is obtained at second order.We therefore truncate the cumulant expansion by setting third order cumulants to zero.These are defined as and so setting ABC = 0 closes the system of differential equations, allowing us to write ABC in terms of first and second order correlations.
C. Behaviour in the thermodynamic limit In Fig. S8 we present the theoretical N -dependence of the charging time τ , maximum energy density E max , and maximum power density P max over a wider range of N than shown in Fig. 3 of the main text.This shows that in addition to the decay-dominated (purple) and coupling-dominated (green) behavior described in the main text, a third region occurs at even larger N , which we discuss below.
In the main text we discussed the energetic dynamics around the decay-to-coupling dominated crossover regime, as this was the experimental operating region.Moving deeper into coupling-dominated regime does not necessarily improve the energy storage properties.This is illustrated in Fig. S8(b) which shows the simulated four points, corresponding to the circles in Fig. S8(a).Within the coupling-dominated regime, energy stored within the microcavity rapidly oscillates which is not a desirable feature.This occurs because the light and matter degrees of freedom hybridise to form polaritons with upper and lower branches split by Rabi frequency ±g √ N , leading to beating between these modes.These oscillations are not present in the experimentally studied crossover region.In this region, dephasing is strong enough to prevent oscillation in energy, yet weak enough to warrant superextensive charging.Therefore, this is the optimal region to produce a QB.Going deeper into the coupling-dominated regime would only be advantageous if energy was extracted on a shorter timescale than the period of oscillations, or additional mechanisms were in place to stabilise the oscillations.
At even larger N (red region) the stored energy falls with increasing N .This can be understood as arising from a condition where the polariton energy splitting exceeds the bandwidth of the pump (set by its finite pulse duration), suppressing energy absorption.Numerically, we find this occurs when N > N σ where g √ N σ = (2/5) 1 4 (1/σ), which signifies the onset of this non-resonant regime.The prefactor (2/5) will be explained in Section S2 D. To build an efficient QB in this regime, one should tune the frequency of the laser to match the polariton energies.Additionally, the time dynamics of energy absorption here change significantly, with the second half of the laser pulse causing stimulated emission, reducing the stored energy -such dynamics arises naturally from a toy model of strongly coupled modes with a splitting larger than the pulse bandwidth, and can be seen in the form of the red line in Fig. S8(b).
In Figure S9, we show the theoretical N -dependence of τ , E max and P max when the frequency of the driving laser is tuned resonant to the lower polariton, i.e. ∆ a = ∆ c = g √ N in the cumulant equations given in Section S2 B. We emphasise that this is not the condition under which the experiments were performed, but of theoretical interest.By comparison of Figures S9 and S8 one can see that the behaviour of the measures with the different driving frequencies are the same until N > N σ .In Figure S9 when N > N σ , the laser frequency continues to drive at the frequency of the lower polariton, instead of at the molecular energy as in Figure S8.In this case, the total energy and power in the cavity continue to grow linearly with N , and so the energy and power densities are constant.

D. The boundary between decay and coupling dominant regimes
There are two timescales in this system: the vacuum Rabi splitting (i.e.polariton detuning) g √ N and the Rabi splitting g √ rN where rN is the number of photons in the cavity.In Fig. 3(a), we show that the microcavity charges super-extensively once g √ N is greater than all decay channels.However, this is only true if r ≤ 1, as is true in experiments A1, A2 and A3.In Fig. 3(b), the boundaries N κ and N γ z are instead determined by g √ rN being equal to the decay rates.This is because r ≥ 1 in experiments B1 and B2.More generally, the important timescale is the larger of the polariton detuning and the Rabi splitting, and so the coupling dominant regime occurs when g Max(1, r)N is larger than all decay channels.
In Figure S10 we plot the charging time τ as a function of N and r.Here, we set κ = γ − = γ z ≡ Γ = 2 meV (note that γ z is independent of N ) so that there is only one boundary between the decay dominant and coupling dominant regimes.The green, red and dashed-black lines show the boundaries between the decay dominant and coupling dominant regimes (N = N Γ ) if g √ N , g √ rN or g Max(1, r)N are used as the relevant coupling scale respectively.Clearly, the boundary is determined by g Max(1, r)N for all values of r.We also show the boundary between the coupling dominant and non-resonant regimes (N = N σ ) as the cyan line.When r > 1, we find that N σ becomes linearly dependent on r.The prefactor (2/5) 1 4 is necessary for N σ to align with the contours of increased charging time for r > 1.

E. Dependence on laser intensity
Figure S11 shows how capacity, charging time and power vary as a function of laser intensity r at fixed number of molecules N .For small r, we find that the maximum energy and power densities vary linearly with r, while charging time is constant.This simply reflects the total energy in the cavity.The charging time is constant because decay channels still dominate over coherent dynamics.As r is increased beyond r = 1, the important timescale g max(1, r)N begins to scale with r, and so the boundaries separating  the coupling dominant and decay dominant regions N κ and N γ z are pushed to smaller N .When these boundaries become smaller than the number of molecules in the cavity, the charging time begins to scale as 1 / √ r.Additionally, the energy density begins to saturate because there are already many more photons than there are molecules within the cavity.In Figure S11(b) we also plot the experimentally measured energy densities, and we see there is good agreement to the theoretical curve.The coloured points in Figure S11(

S3. FITTING OF MODEL PARAMETERS
As outlined in the main text, we used a reduced chi-squared optimisation procedure to determine the light matter coupling g, dephasing constant γ z 0 and non-radiative decay rate γ − , as well as to estimate uncertainties on these parameters.As these quantities represent molecular properties, we would expect them to be the same in all the different experiments.For this reason, we performed a global fit rather than performing the procedure individually for each experiment.In this section we give further details on this calculation.As also discussed in the main text, our fitting process is key in estimating the charging time, stored energy and peak charging power.Extracting these quantities requires a smooth curve of charge vs time.Any smoothing process implicitly introduces an effective fitting function (such as e.g.fitting the data to piecewise cubic splines).Since our best guess for such a fitting function is in fact the theoretical model described above, we use this continuous function, fit to the experimental data, to extract charging times, energies, and powers.
In performing our fitting, it is necessary to use the values of the molecule numbers, estimated as discussed above.Because of this, errors in the estimates of N affect the fitting parameters, and thus the resulting charging timescales and powers.These correlated errors make it challenging to estimate errors in the powerlaw scaling of charging time vs N .
The cavity lifetime (or equivalently the cavity linewidth κ) can also be considered as a fitting parameter, however the value of this parameter is more strongly constrained by other measurements.As noted above, transfer matrix simulations on the designed cavity give an upper bound of 306fs.An alternate estimate is provided by comparing the theoretical and measured reflectivity spectra of the cavities.As discussed below, this implies a cavity lifetime of 120fs.In the following, we first present fitting results for a cavity lifetime of 120fs, and then in Sec.S3 B we discuss how the results change with alternate cavity lifetimes in the range 120fs to 306fs.

A. Fitting procedure
The steps for our fitting procedure are as follows: 1. Calculate the theoretical E(t) curves for a grid of parameter values g, γ z , γ − , along with the values of N relevant for all five experiments, A1, A2, A3, B1 and B2.Based on previous observations, we chose the search region of the parameter space as g ∈ [0.1, 5000] neV; γ z 0 ∈ [0.1, 5000] meV and γ − = [0.001,1] meV.Subsequent refinements of this search region were made to give higher resolution near the optimal point.2. We estimate uncertainties, σ i on each experimental data point (transient reflectivity vs time), by considering the point-to-point variation.Because the uncertainty is higher near t = 0, when the pump arrives, we use different error estimates in different time windows.Specifically, we divide experiments A1 and A2 into five windows t < −300fs; −300fs < t < 300fs; 300fs < t < 700fs; 700fs < t < 1000fs; t > 1000fs.For experiments A3, B1 and B2 we found that four windows t < −300fs; −300fs < t < 300fs; 300fs < t < 1000fs; t > 1000fs was sufficient.In each window, the uncertainty estimate for each experiment is taken from the variance over a narrow range of points (typically 150fs) where there is no strong time dependence.

3.
For each set of parameters, we performed an "internal" chi-squared minimisation to find the optimal scaling factor S between the stored energy E(t) and the measured differential reflectivity ∆R/R, and a time shift between the theory and experiment T 0 .That is, we minimise with respect to S and T 0 .We treat the result of this minimisation as the chi-squared value which we use in the following steps to estimate the meaningful parameters g, γ z , γ − and their uncertainties.
Estimating the scaling factor S from first principles is difficult because of reflections by the cavity mirror, hence this factor is found by the best fit value.The time shift reflects uncertainty of delays in the optics, so that it is not a-prior clear when the peak of the pump pulse arrives.After this shift, we define t = 0 as the moment the pump arrives.This is important when calculating the charging time τ , which we defined as the time from the arrival of the pump until reaching half maximum energy.
4. We then use the chi-squared value described above, and divide by the total number of degrees of freedom k eff = k − 3 (where k is the total number of data points), to arrive at the final reduced chisquared χ2 map.A slice of this three dimensional reduced chi-squared map for the 120 fs lifetime is shown in Figure S12 for γ − = 0.0263 meV, which is the optimal non-radiative decay rate given in the main text.The optimal parameter set used in the main text that optimises χ2 is shown as the red point in Figure S12.We find χ2 min = 3.048, suggesting our estimated measurement uncertainties on χ are reasonable, but likely underestimates.5. Finally, the 68% confidence interval for each parameter was estimated by considering the contour for which χ2 = χ2 min + 1 k eff ∆ * where ∆ * = 3.51 is extracted from the reduced chi-squared distribution for 3 parameters and error tolerance (68%), see [48].In the right panel of Fig. S12 we show the contour as a white line, and the actual parameter values which lie within this 68% contour as black points.

B. Fitting cavity lifetime
Figure S13 shows the optimal reduced chi-squared as a function of cavity lifetime, along with the corresponding best-fit values of the parameters g, γ z 0 and γ − , following the fitting procedure described above.The minimal reduced chi-squared is 2.803 occurring for a lifetime of 185 fs.
As noted earlier, the cavity lifetime is also constrained by the measured reflectivity spectrum, shown in Figure S2(d).To check this consistency, Fig. S14 shows the calculated absorption spectra for the 0.5%, 1%, 5% and 10% cavities using the optimal parameter sets given by the 185 fs and 120 fs lifetimes.These are where Ω eff = g 2 N − (κ − 2γ tot ) 2 /4 is the effective Rabi splitting.From the measured spectra in Absorption spectra for the 0.5%, 1%, 5% and 10% cavities.The spectrum is calculated using Eq. ( 22) for the best-fit parameters (see Figure S13) for both 120 fs and 185 fs cavity lifetimes.In the 5% and 10% cavities, there are clear polariton peaks forming at ∆ν = ±g √ N .
It is clear from Figure S14 that although the 185 fs cavity lifetime gives a smaller reduced chi-squared, these parameters predict significantly stronger coupling than is seen in experimental reflectivity spectra.In contrast, the 120 fs cavity lifetime parameters reproduce the polariton splittings across all cavities more accurately, while showing a reduced chi-squared that is not significantly larger.We therefore conclude that 120 fs cavity lifetime, with g, γ z 0 and γ − given in the main text provides the best fit to the experimental data.In Table S2, we summarise the scaling S and time shifts T 0 that relate the measured differential reflectivity ∆R/R to the theoretically calculated stored energy E(t), as used to calculate the fits plotted in Figure 2(b) of the main text.For a formal defintion of S, T 0 , see Eq. ( 21).The theoretical time evolution arising from using the above fitting procedure is shown in Figure 3 of the main text.In plotting that figure, the results of the theoretical fit are convolved with an instrument response function, as required to match the experimental data.Figure S15 shows the same data but without convolution by the instrument response, thus providing a more direct picture of the intrinsic dynamics of the system.From the theoretical curves, one can extract the rise time of stored energy τ , the peak stored energy E max , and maximum charging power P max .These values are summarised in Table 1 in the main text.
As well as the scaling to the observables with N shown in Figure 3 in the main text, we could also estimate an effective power-law scaling of the observables q i ∈ {τ, E max , P max } directly from pairs of experiments i, j by the relation q i /q j = (N i /N j ) fq .As all our observables are intensive (i.e.densities), f q > 0 indicates superextensive behaviour, f q = 0 indicates extensive, and f q < 0 subextensive behaviours.Table S3 gives the observed values of f q .

D. Residuals of the best fit
To check whether systematic errors arise from our fitting procedure, Fig. S16 shows the residual errorsi.e.difference between the theoretical curves and the experimental data-for the five experiments shown in

Fig. 1
Fig. 1 Schematics of the LFO microcavity and experimental setup.(a) Microcavity consisting of Lumogen-F Orange (LFO) dispersed in a polystyrene (PS) matrix between distributed Bragg reflectors (DBRs).(b) Normalised absorption (red) and photoluminescence (blue) spectra for 1% concentration LFO film, with the molecular structure shown in the inset.We operate near peak absorption/photoluminescence.(c) Angle-dependent reflectivity of the 1% cavity, with a fit for the cavity mode shown by the blue dashed line.(d) A laser pump pulse excites the LFO molecules.The energetics of the molecules are then measured with probe pulses delayed by time t, from which we can ascertain the peak energy density ( max ), rise time (), and peak charging power ( max ).(e) The experimental setup for ultrafast transient reflectivity measurements.The output of a non-collinear optical parametric amplifier (NOPA) is split to generate pump (dark green) and probe (light green) pulses.A mechanical chopper is used to modulate the pump pulse to produce alternating pump-probe and probe-only pulses.
Fig.1b, we attribute the 530 nm band to ground state bleaching (GSB)i.e.suppression

Fig. 3 :
Fig. 3: Charging dynamics as function of the number of molecules.(a) and (b) show S2.For the 5% cavity, we see a weak anti-crossing feature in the reflectivity spectrum (a small kink near the crossing), indicating operation in the intermediate coupling regime.The 10% cavity operated in the strong-coupling regime, showing a Rabi splitting of around 100 meV around the 0-0 transition (along with intermediate-coupling between the cavity mode and the 0-1 transition).
is the Lindbladian superoperator. † and  are the cavity photon creation and annihilation operators, and   ,, are the Pauli spin matrices for each molecule, with the raising and lowering spin operators defined as   ± = (   ±    )/2 .There are three decay channels corresponding to the cavity decay (), dephasing (  ), and relaxation rate ( − ) of the individual TLSs.The Hamiltonian for the LFO molecules in cavity is modelled as a collection of non-interacting TLSs with characteristic frequency  equal to that of the cavity mode, and resonantly coupled to the cavity with strength .The molecules are driven by a laser described by a Gaussian pulse envelope ()

(
0  ,  − and ) were found through a global chi-squared optimisation, simultaneously optimising over all experiments.Uncertainties in these fitting parameters were then estimated by using the reduced  ̃2 distribution to find the 68% confidence interval of the model parameters.This corresponds to the range  ̃2 ≤  ̃ 2 + ∆, where for a three parameter optimisation and  total data points ∆≈ 3.51 ( − 3) ⁄ .(48)In the Supplementary Information we present a figure showing minimum reduced chi-squared value as a function of , and for each point we show the optimal set of parameters ( 0  ,  − , ) along with the 68% confidence intervals.From this, and by comparison of the experimentally measured and theoretically calculated reflectivity for each parameter set, we concluded that the lifetime most representative of the data was  = 120 fs with  − = S1.CHARACTERISATION OF SAMPLES AND CALIBRATION MEASUREMENTS FIG. S1.Photoluminescence quantum yield as a function of LFO concentration.

1 FIG
FIG. S2.Absorption, photoluminescence of the LFO films and reflectivity spectra of the microcavities.(a) Absorption and (b) photoluminescence spectra for the 0.5%, 1%, 5%, and 10% LFO-concentration films.(c) Reflectivity spectra for 0.5%, 1%, 5%, and 10% LFO-concentration microcavities.UPB, MPB, and LPB label the upper, middle, and lower polariton branches, respectively.Also indicated are the 0-0 and 0-1 transition wavelengths.(d) is a slice of the reflectivity spectra at 37 • .The single dip in the 0.5% and 1% concentration spectra indicate the weak-coupling regime.The double dip seen in the 10% concentration spectra, represent the polaritonic states, indicating the strong-coupling regime.The 5% concentration spectrum represents a situation intermediate between a single and double dip, indicating an intermediate-coupling regime.

FIG. S3 .
FIG. S3.Transfer matrix simulation of the electric field distribution for the 1% cavity.(a) shows the cavity reflectivity, (b) the spectrally-resolved electric field amplitude, and (c) the electric field amplitude at the cavity mode wavelength.The shaded sections of (c) indicate the different materials which make up the cavity, with Nb2O5 in purple (refractive index, n = 2.25), SiO2 in blue (n = 1.52), and LFO in PS in orange (n = 1.60).All simulations were made using transfer matrix modelling at an angle of 20 • to the cavity normal to maintain consistency with the transient reflectivity measurements.
FIG. S4.Dynamics of the control 1% control film at different wavelengths.The two panels show two different time windows (top panel until 10 ps, bottom panel until 1.4 ps.The same behavior is observed for other film concentrations. FIG. S6.Calculating photon number.Yellow curve (right axis): pump spectrum.Black curve (left axis) reflectivity spectrum of 1% cavity.Purple curve (right axis): photons transmitted into the cavity.
FIG. S7.Time evolution of positive and negative features in the differential transmitivity.
FIG. S8.Charging dynamics vs number of molecules N .(a) Charging time, peak stored energy, and maximum power as a function of N .This figure is identical to Fig. 3(a) in the main text, but extended to larger range of N .At large N the stored energy can reach half maximum before the laser pulse finishes, in which case the charging time τ becomes negative.(b) Examples of dynamics in each regime.The values of N used in these dynamics corresponds to the circles of the same colour in (a).
FIG. S9.Charging dynamics when pumping at lower polariton.Charging time, peak stored energy, and maximum power as a function of N .This figure is identical to Fig. S8 except that the frequency of the laser is tuned to the lower polariton energy, ∆a = ∆c = g √ N , rather than the molecular energy.
FIG. S10.Charging time as a function of number of molecules N and laser intensity r.All parameters are equivalent to the Q1% cavity (see main text) with the exception that the dephasing and non-radiative decay rates are equal to the cavity leakage rate (set to γ z = κ = Γ = 2 meV) and note that the dephasing rate is independent of N .
a) indicate the charging time, maximum capacity and maximum power of the temporal dynamics of the same colour in Figure S11(b).
FIG. S11.Charging dynamics vs pump intensity r.(a) Capacity, charging time and maximum power vs r for a fixed number of molecules N = 1.62 × 10 10 (the 1% cavity).All other parameters align with those used to model the 1% cavity.(b) Comparison of theoretical and experimental charging dynamics for four values of r.These are indicated by circles in the charging time panel of (a), showing the experimentally measured values.

3.21 3 .
FIG. S12.Reduced-chi square map to find the optimal parameters for the theoretical model and their 68% confidence intervals.The chi-squared contour plots shown in this figure are slices of the full three-dimensional map at the optimum non-radiative decay rate γ − = 0.0141 meV used in the main text for a cavity lifetime of 120 fs.In the yellow region of the bottom right corner χ2 > 5, which we do not show to emphasise smaller variations in χ2 .In the right panel, the highlighted contour shows the 68% confidence interval, found using ∆ = ∆ * /k eff .
Temporally resolved energy density of the microcavities shows that rise time decreases

Table 1 : Summary of the experimental results. In
≈ 0.105 and 2.4 respectively).The rise time , is defined by the time to reach   /2, where   is the peak stored energy per molecule or energy density.The charging rate

TABLE S1
. Estimated photon number and molecule number for each experiment.

TABLE S2 .
The optimal scaling factors and time shifts used to calculate the theoretical curves in Figure2(b) in the main text.
C. Results of fitting procedure