Thomas F Allard

Guillaume Weick

Disorder-enhanced transport in a chain of lossy dipoles strongly coupled to cavity photons

We study the interplay between disorder and light-matter coupling by considering a disordered onedimensional chain of lossy dipoles coupled to a multimode optical cavity, through a microscopically derived Hamiltonian. Such a system, hosting polaritonic excitations, may be realized experimentally in a wide range of platforms under strong light-matter coupling. By analyzing both the eigenspectrum and the driven-dissipative transport properties of our system, we find that in the strong-coupling regime, increasing disorder leads almost uncoupled dark states to acquire a photonic part, allowing them to inherit polaritonic long-range transport characteristics. Crucially, we show that this disorder-enhanced transport mechanism is increasingly noticeable when the considered dipoles are lossier.

I. INTRODUCTION

Over the past several years, strong coupling of matter excitations with confined electromagnetic modes has been shown to significantly modify material properties [START_REF] Garcia-Vidal | Manipulating matter by strong coupling to vacuum fields[END_REF]. In particular, it has been proposed both theoretically [START_REF] Feist | Extraordinary exciton conductance induced by strong coupling[END_REF][START_REF] Schachenmayer | Cavity-enhanced transport of excitons[END_REF][START_REF] Gonzalez-Ballestero | Uncoupled dark states can inherit polaritonic properties[END_REF][START_REF] Hagenmüller | Cavity-enhanced transport of charge[END_REF][START_REF] Botzung | Dark state semilocalization of quantum emitters in a cavity[END_REF][START_REF] Chávez | Disorder-enhanced and disorderindependent transport with long-range hopping: Application to molecular chains in optical cavities[END_REF][START_REF] Dubail | Large random arrowhead matrices: Multifractality, semilocalization, and protected transport in disordered quantum spins coupled to a cavity[END_REF][START_REF] Gera | Effects of disorder on polaritonic and dark states in a cavity using the disordered Tavis-Cummings model[END_REF][START_REF] Engelhardt | Unusual dynamical properties of disordered polaritons in microcavities[END_REF][START_REF] Patton | Allphotonic quantum simulators with spectrally disordered emitters[END_REF][START_REF] Sun | Dynamics of disordered Tavis-Cummings and Holstein-Tavis-Cummings models[END_REF][START_REF] Ribeiro | Multimode polariton effects on molecular energy transport and spectral fluctuations[END_REF] and experimentally [START_REF] Lerario | Room temperature Bloch surface wave polaritons[END_REF][START_REF] Orgiu | Conductivity in organic semiconductors hybridized with the vacuum field[END_REF][START_REF] Zhong | Energy transfer between spatially separated entangled molecules[END_REF][START_REF] Lerario | High-speed flow of interacting organic polaritons[END_REF][START_REF] Rozenman | Long-range transport of organic exciton-polaritons revealed by ultrafast microscopy[END_REF][START_REF] Hou | Ultralong-range energy transport in a disordered organic semiconductor at room temperature via coherent exciton-polariton propagation[END_REF][START_REF] Nagarajan | Conductivity and photoconductivity of a p-type organic semiconductor under ultrastrong coupling[END_REF][START_REF] Pandya | Tuning the coherent propagation of organic exciton-polaritons through dark state delocalization[END_REF][START_REF] Xu | Ultrafast imaging of coherent polariton propagation and interactions[END_REF][START_REF] Balasubrahmaniyam | Unveiling the mixed nature of polaritonic transport: From enhanced diffusion to ballistic motion approaching the speed of light[END_REF] to use the strong light-matter coupling between emitters and cavity photons to modify transport properties and induce long-range transport, especially in disordered, molecular systems. This cavity-enhanced transport effect can be understood from the fact that in the strong-coupling regime, the matter excitations hybridize with cavity photons, leading to collective polaritonic excitations that are delocalized throughout the whole system.

Notably, the interplay between strongly coupled lightmatter systems and disorder has been shown to lead to unexpected phenomena, such as a cavity-protection effect, namely that polaritonic states are more robust than others against disorder [START_REF] Houdré | Vacuum-field Rabi splitting in the presence of inhomogeneous broadening: Resolution of a homogeneous linewidth in an inhomogeneously broadened system[END_REF][START_REF] Michetti | Polariton states in disordered organic microcavities[END_REF]. Such a behavior goes beyond the well-known theory of Anderson localization with short-range interaction, which states that in one-dimensional (1d) systems, disorder makes all the eigenstates exponentially localized, thereby suppressing transport [START_REF] Anderson | Absence of diffusion in certain random lattices[END_REF][START_REF] Evers | Anderson transitions[END_REF]. For long-range interaction, however, the situation becomes highly nontrivial, since the 1d system can host extended states [START_REF] Levitov | Delocalization of vibrational modes caused by electric dipole interaction[END_REF][START_REF] Mirlin | Transition from localized to extended eigenstates in the ensemble of power-law random banded matrices[END_REF][START_REF] Rodríguez | Anderson transition in low-dimensional disordered systems driven by long-range nonrandom hopping[END_REF][START_REF] De Moura | Localization properties of a one-dimensional tight-binding model with nonrandom longrange intersite interactions[END_REF][START_REF] Deng | Duality in power-law localization in disordered one-dimensional systems[END_REF]. Since polaritons are partly mediated by cavity photons, hence through an effective long-range coupling, they are less affected by the suppression of transport arising from Anderson localization. Interestingly, dark states, i.e., matter states that are only weakly hybridized to photonic degrees of freedom, can also show nontrivial transport properties induced by the strong-coupling regime [START_REF] Gonzalez-Ballestero | Uncoupled dark states can inherit polaritonic properties[END_REF][START_REF] Botzung | Dark state semilocalization of quantum emitters in a cavity[END_REF][START_REF] Dubail | Large random arrowhead matrices: Multifractality, semilocalization, and protected transport in disordered quantum spins coupled to a cavity[END_REF][START_REF] Engelhardt | Unusual dynamical properties of disordered polaritons in microcavities[END_REF][START_REF] Pandya | Tuning the coherent propagation of organic exciton-polaritons through dark state delocalization[END_REF].

Alongside the above-mentioned cavity-protection and cavity-enhanced transport effects, recent works [START_REF] Chávez | Disorder-enhanced and disorderindependent transport with long-range hopping: Application to molecular chains in optical cavities[END_REF][START_REF] Dubail | Large random arrowhead matrices: Multifractality, semilocalization, and protected transport in disordered quantum spins coupled to a cavity[END_REF][START_REF] Engelhardt | Unusual dynamical properties of disordered polaritons in microcavities[END_REF] have shown that the photonic cavity can also lead to an improvement of the transport characteristics when increasing the disorder strength, instead of an expected suppression. This counterintuitive phenomenon was theoretically unveiled by considering a disordered single-mode Tavis-Cummings model, namely a chain of emitters coupled to a single cavity mode through a spatially homogeneous coupling constant. Such enhanced transport induced by disorder has been very recently observed in experiments [START_REF] Cohn | Vibrational polaritons in disordered molecular ensembles[END_REF][START_REF] Sauerwein | Engineering random spin models with atoms in a highfinesse cavity[END_REF][START_REF] Baghdad | Spectral engineering of cavity-protected polaritons in an atomic ensemble with controlled disorder[END_REF].

The interplay between disorder and strong light-matter cou-

L z L x x ŷ ẑ ! 1 ! 2 ! i ! N d N 1 d i d 2 d 1 FIG.
1. Sketch of a disordered chain of N oscillating dipoles polarized along the x axis and arranged along the z direction of a cuboidal cavity of linear sizes Lx, Ly, and Lz. Each dipole of frequency ωi, depicted by a red arrow, is located at a distance di from its neighbor to the right. pling is therefore highly nontrivial, and its understanding is of primary importance, since disorder is always present in experimental setups. Here, we address such interplay by considering a finite disordered 1d chain of N oscillating dipoles placed inside a mirror cavity (see Fig. 1). We keep our model general such that these dipoles represent a wide range of physical systems whose main coupling mechanism is of dipolar nature, and which are not dominated by quantum effects. These range from plasmonic, dielectric or SiC nanoparticles [START_REF] Krenn | Squeezing the optical near-field zone by plasmon coupling of metallic nanoparticles[END_REF][START_REF] Maier | Local detection of electromagnetic energy transport below the diffraction limit in metal nanoparticle plasmon waveguides[END_REF][START_REF] Koenderink | Experimental evidence for large dynamic effects on the plasmon dispersion of subwavelength metal nanoparticle waveguides[END_REF][START_REF] Markel | Propagation of surface plasmons in ordered and disordered chains of metal nanospheres[END_REF][START_REF] Crozier | Experimental measurement of the dispersion relations of the surface plasmon modes of metal nanoparticle chains[END_REF][START_REF] Apuzzo | Observation of near-field dipolar interactions involved in a metal nanoparticle chain waveguide[END_REF][START_REF] Slobozhanyuk | Subwavelength topological edge states in optically resonant dielectric structures[END_REF][START_REF] Wang | Topological phonon polaritons in one-dimensional non-Hermitian silicon carbide nanoparticle chains[END_REF][START_REF] Zhang | Second-order topological photonic modes in dipolar arrays[END_REF] to magnonic microspheres [START_REF] Pirmoradian | Topological magnon modes in a chain of magnetic spheres[END_REF][START_REF] Rameshti | Cavity magnonics[END_REF] or macroscopic microwave helical antennas [START_REF] Mann | Manipulating type-I and type-II Dirac polaritons in cavityembedded honeycomb metasurfaces[END_REF][START_REF] Mann | Tunable pseudomagnetic fields for polaritons in strained metasurfaces[END_REF]. Our model can also describe molecular and semiconductor excitons [START_REF] Feist | Extraordinary exciton conductance induced by strong coupling[END_REF][START_REF] Yuen-Zhou | Plexciton Dirac points and topological modes[END_REF], for which most of the recent experimental approaches on disordered polaritonic systems were performed [START_REF] Lerario | Room temperature Bloch surface wave polaritons[END_REF][START_REF] Orgiu | Conductivity in organic semiconductors hybridized with the vacuum field[END_REF][START_REF] Zhong | Energy transfer between spatially separated entangled molecules[END_REF][START_REF] Lerario | High-speed flow of interacting organic polaritons[END_REF][START_REF] Rozenman | Long-range transport of organic exciton-polaritons revealed by ultrafast microscopy[END_REF][START_REF] Hou | Ultralong-range energy transport in a disordered organic semiconductor at room temperature via coherent exciton-polariton propagation[END_REF][START_REF] Nagarajan | Conductivity and photoconductivity of a p-type organic semiconductor under ultrastrong coupling[END_REF][START_REF] Pandya | Tuning the coherent propagation of organic exciton-polaritons through dark state delocalization[END_REF][START_REF] Xu | Ultrafast imaging of coherent polariton propagation and interactions[END_REF][START_REF] Balasubrahmaniyam | Unveiling the mixed nature of polaritonic transport: From enhanced diffusion to ballistic motion approaching the speed of light[END_REF], as well as cold atoms [START_REF] Browaeys | Experimental investigations of dipole-dipole interactions between a few Rydberg atoms[END_REF][START_REF] Barredo | An atom-by-atom assembler of defect-free arbitrary two-dimensional atomic arrays[END_REF][START_REF] Perczel | Topological quantum optics in two-dimensional atomic arrays[END_REF][START_REF] Wang | Topological photonic states in onedimensional dimerized ultracold atomic chains[END_REF] where controllable frequency disorder has been achieved [START_REF] Sauerwein | Engineering random spin models with atoms in a highfinesse cavity[END_REF][START_REF] Baghdad | Spectral engineering of cavity-protected polaritons in an atomic ensemble with controlled disorder[END_REF], and which, in the single excitation manifold, behave as classical dipoles [START_REF] Asenjo-Garcia | Optical waveguiding by atomic entanglement in multilevel atom arrays[END_REF].

We consider dipoles that are coupled to each other through an all-to-all short-range quasistatic dipole-dipole interaction, and, importantly, we go beyond the usually considered single mode Tavis-Cummings model by studying a microscopically derived minimal coupling Hamiltonian [START_REF] Downing | Topological phases of polaritons in a cavity waveguide[END_REF][START_REF] Downing | Polaritonic Tamm states induced by cavity photons[END_REF] which couples the dipoles to a multimode cavity. Since the majority of materials used in experimental polaritonics, such as Fabry-Pérot or plasmonic cavities, are intrinsically multimodal, it is of fundamental importance to go beyond single-mode models [START_REF] Ribeiro | Multimode polariton effects on molecular energy transport and spectral fluctuations[END_REF][START_REF] Tichauer | Multi-scale dynamics simulations of molecular polaritons: The effect of multiple cavity modes on polariton relaxation[END_REF]. We take into account disorder both in the individual resonance frequencies of each dipole, modeling possible inhomogeneities and in the interdipole distances, accounting for positioning uncertainties. Crucially, unlike the recent studies of Refs. [START_REF] Chávez | Disorder-enhanced and disorderindependent transport with long-range hopping: Application to molecular chains in optical cavities[END_REF][START_REF] Dubail | Large random arrowhead matrices: Multifractality, semilocalization, and protected transport in disordered quantum spins coupled to a cavity[END_REF][START_REF] Engelhardt | Unusual dynamical properties of disordered polaritons in microcavities[END_REF], we study transport through an out-ofequilibrium driven-dissipative scenario, which permits us to analyze the propagation of polaritonic excitations along the chain as a function of the driving frequency, and to study the effect of dissipation, taking into account losses both in the dipoles and in the mirror cavity.

Our model allows us to show that disorder-enhanced transport can be readily explained as a result of an hybridization between dark and polaritonic states, which is induced by the bandwidth increase led by frequency disorder. This disorderinduced hybridization enables dark states to inherit polaritonic properties, namely to take advantage of cavity-enhanced transport. Importantly, we show that several different transport regimes are found for chains of large-enough size, so that cavity-and disorder-enhanced transport only occur at intermediate distances. Indeed, the short-range propagation of the dark states, which results from the nearest-neighbor quasistatic dipole-dipole coupling only, cannot be enhanced through polaritonic hybridization, and only the long-range propagation can profit from cavity-and disorder-enhanced transport. Moreover, no enhancement is found over very long distances. Crucially, the consideration of dipolar losses is proved to be of great importance, since the latter emphasizes the aforementioned effects of cavity-and disorder-enhanced transport.

The paper is organized as follows: Section II is dedicated to the presentation of our microscopic model Hamiltonian of a disordered chain of dipoles coupled to a cuboidal photonic cavity. In Sec. III, we study the localization properties of the system from an eigenspectrum analysis, both in the case of ordered and disordered chains, where we discuss the mixing between dipolar and photonic degrees of freedom induced by the interplay between strong light-matter coupling and disorder. In Sec. IV, we study the transport along the chain in a drivendissipative scenario, considering both lossy dipoles and lossy cavity mirrors. We first study the effects of strong coupling in an ordered chain, and then the interplay between disorder and cavity-enhanced transport in a disordered chain. Finally, we summarize our results and draw conclusions in Sec. V.

In Appendix A we reproduce for the sake of selfcontainedness results for the case of an ordered chain [START_REF] Downing | Polaritonic Tamm states induced by cavity photons[END_REF], while in Appendix B we discuss some approximations adopted in our model. In Appendixes C, D, E, and F, we complement the discussion of the main text by studying the precise nature of the semilocalized eigenstates when increasing the system size, giving technical details of our transport simulations, investigating the influence of cavity losses, and presenting absorption spectra obtained from transport calculations, respectively.

II. DISORDERED CHAIN OF DIPOLES COUPLED TO CAVITY PHOTONS

To describe a disordered chain of N dipoles coupled to cavity photons, we rely on the model recently developed in Refs. [START_REF] Downing | Topological phases of polaritons in a cavity waveguide[END_REF][START_REF] Downing | Polaritonic Tamm states induced by cavity photons[END_REF] for an ordered chain, and supplement it by including disorder both in the interdipole spacings d i and in the individual resonance frequencies ω i (see Fig. 1). Our model is described by the Hamiltonian

H = H dp + H ph + H dp-ph . (1)
The first term on the right-hand side (r.h.s.) of Eq. ( 1) accounts for the dipolar degrees of freedom and reads [START_REF] Downing | Retardation effects on the dispersion and propagation of plasmons in metallic nanoparticle chains[END_REF]

H dp = N i=1 ω i b † i b i + N -1 i=1 N j=i+1 Ω ij b † i b j + b i b † j , (2)
where counter-rotating terms in the quasistatic dipole-dipole coupling have been ignored. A treatment of the Hamiltonian (1) without using the rotating wave approximation is discussed in Appendix B. In the above dipolar Hamiltonian (2), the bosonic operator b i (b † i ) annihilates (creates) a dipolar excitation, polarized along the x axis and with resonance frequency ω i on site i ∈ [1, N ]. The latter individual resonance frequencies are uncorrelated random variables distributed uniformly within an interval [ω 0 -W/2, ω 0 +W/2], hence with an average resonance frequency being ω 0 . The dipolar coupling strength Ω ij between two dipoles on sites i and j separated by a distance r ij = j-1 l=i d l is given by

Ω ij = Q 2 2M √ ω i ω j r 3 ij , (3)
where the interdipole distances d i are also uncorrelated random variables distributed uniformly within an interval [d 0 -∆/2, d 0 + ∆/2], the average spacing being d 0 . The two dimensionless parameters W/ω 0 and ∆/d 0 hence govern the disorder strength on the dipole frequencies and on the spacings between two neighboring dipoles, respectively. In our microscopically derived model, the disorder parameters W/ω 0 and ∆/d 0 are restricted, respectively, to values smaller than 2 and smaller than 2 -6a/d 0 , in order to maintain positive resonance frequencies ω i and spacings d i 3a which allow us to safely neglect multipolar effects and describe purely dipolar excitations [START_REF] Park | Surface-plasmon dispersion relations in chains of metallic nanoparticles: An exact quasistatic calculation[END_REF]. In Eq. ( 3), Q < 0 and M are the effective charge and effective mass of a dipole. For convenience, we recast the two latter parameters in a single one by defining the dipole length scale a = (Q2 /M ω 2 0 )1/3 . 1 The average dipolar coupling between two neighboring dipoles is then Ω 0 = Ω i,i+1 = (ω 0 /2)(a/d 0 ) 3 , where the bar denotes averaging over the disorder realizations. In the remainder of our paper, we choose an average spacing d 0 = 4a, so that the dipoles are coupled in the near-field regime, the one most studied both experimentally [36-38, 40, 41, 60, 61] and theoretically [START_REF] Markel | Propagation of surface plasmons in ordered and disordered chains of metal nanospheres[END_REF][START_REF] Downing | Retardation effects on the dispersion and propagation of plasmons in metallic nanoparticle chains[END_REF][START_REF] Koenderink | Complex response and polariton-like dispersion splitting in periodic metal nanoparticle chains[END_REF][START_REF] Brandstetter-Kunc | Nonradiative limitations to plasmon propagation in chains of metallic nanoparticles[END_REF]. Such a small interdipole spacing leads to Ω 0 /ω 0 = 1/128.

The second term on the r.h.s. of the Hamiltonian (1) describes the photonic degrees of freedom and reads

H ph = k, λk ω ph k c λk k † c λk k , (4)
where the bosonic ladder operator c λk k (c λk k †

) annihilates (creates) a photon with wavevector k and transverse polarization λk , i.e., k • λk = 0. We consider these photonic degrees of freedom inside a mirror cavity of cuboidal geometry with perfectly conducting walls on all sides, with linear sizes L x , L y , and L z , as sketched in Fig. 1. This amounts to consider hardwall boundary conditions, and the cavity photon wavevector along the direction σ = x, y, z takes the form k σ = πn σ /L σ , with n σ ∈ N. The photonic dispersion relation then reads

ω ph k = ω ph nxnynz = c πn x L x 2 + πn y L y 2 + πn z L z 2 , (5)
with c being the speed of light in vacuum.

Finally, the last term on the r.h.s. of Eq. ( 1) represents the minimal coupling Hamiltonian between dipolar excitations and cavity photons. In the point-dipole approximation, where we consider a dipole length scale a much smaller than the average inverse wavevector k -1 0 = c/ω 0 associated with an individual dipolar mode, 2 we have [START_REF] Craig | Molecular Quantum Electrodynamics[END_REF]

H dp-ph = Q M c N i=1 Π i • A(r i ), (6)
where r i = xx + y ŷ + z i ẑ corresponds to the location of the center of the dipole on site i. By placing the chain at the center of the cavity and separated from the walls in the z direction by a distance d 0 , we have x = L x /2, y = L y /2 and

z i = i-1 l=0 d l .
The momentum associated with the ith dipole excitation, polarized along the x axis, reads

Π i = i M ω i 2 b i † -b i x, (7)
while the vector potential quantized in a cuboidal cavity with hard-wall boundary conditions on all sides reads in the Coulomb gauge [START_REF] Kakazu | Quantization of electromagnetic fields in cavities and spontaneous emission[END_REF] A

(r i ) = k, λk 2π c 2 ω ph k f λk k (r i )c λk k + f λk k * (r i )c λk k † , (8)
2 For this purpose, we choose to fix the value of the dimensionless dipole strength k 0 a to 0.1 in the remainder of this paper.

where the mode functions are

f λk k (r i ) = C   cos(k x x) sin(k y y) sin(k z z i )(x • λk ) sin(k x x) cos(k y y) sin(k z z i )(ŷ • λk ) sin(k x x) sin(k y y) cos(k z z i )(ẑ • λk )   , (9)
with C = 2/ √ V if k x , k y , or k z = 0, and C = 2 2/V otherwise, V = L x L y L z being the volume of the cavity. Note that in our modeling of the light-matter interaction, we do not consider dipole-photon Umklapp processes, the latter being irrelevant in the near-field regime, when the dipole spacings d i λ 0 = 2π/k 0 [START_REF] Downing | Topological phases of polaritons in a cavity waveguide[END_REF][START_REF] Allard | Quantum theory of plasmon polaritons in chains of metallic nanoparticles: From near-to far-field coupling regime[END_REF]. In order to simplify the minimal coupling Hamiltonian (6), we use the fact that we consider a number of dipoles N 1, such that the longitudinal size of the cavity L z = z N + d 0 L x , L y , and choose a cavity aspect ratio L y = 3L x . We further consider only the resonance between the photonic mode (n x , n y , n z ) = (0, 1, 0) and the average dipole resonance frequency, i.e., ω ph 010 ∼ ω 0 , so that the y and z components of the mode functions (9) vanish. This condition is fulfilled for small cavity heights 3a L x 20a and it allows us to neglect all of the higher frequency modes such that ω ph nxny0 3ω 0 , which we consider being out of resonance and which, due to the above choice of geometry, are all of the (n x , n y ) = (0, 1). Hence, as justified in greater details in the supplementary materials of Refs. [START_REF] Downing | Topological phases of polaritons in a cavity waveguide[END_REF][START_REF] Downing | Polaritonic Tamm states induced by cavity photons[END_REF] (in the case of periodic boundary conditions in the z direction), this resonance condition allows us to keep only the low frequency photonic modes (0, 1, n z ), where n z is chosen in [1, N ], and only one of the two photon polarizations λk is relevant, so that we no longer display the polarization label.

Within the above-discussed hypotheses and using in addition the rotating wave approximation, which we discuss in Appendix B, the Hamiltonian (1) becomes

H = H dp + N nz=1 ω ph nz c † nz c nz + i N i=1 N nz=1 ξ inz b † i c nz -b i c † nz , (10)
where the photonic dispersion (5) simplifies to

ω ph nz = c π L y 2 + πn z L z 2 , (11)
the light-matter coupling writes

ξ inz = ω 0 4πa 3 ω i Vω ph nz sin πn z L z z i , (12)
and with the dipolar Hamiltonian H dp being given in Eq. ( 2).

In the following, we will first study the properties of the spectrum and eigenstates of the approximate Hamiltonian [START_REF] Engelhardt | Unusual dynamical properties of disordered polaritons in microcavities[END_REF], and then study the propagation along the chain in a driven-dissipative transport scenario, considering losses both in the dipoles and in the cavity mirrors. Furthermore, we ±π/d 0 ,-. The blue and green dashed lines represent the bare dipolar dispersion ω dp q [see Eq. (A2)] and the bare photonic one ω ph q [see Eq. (A3)], respectively. The orange, red, and blue symbols point out eigenstates that will be studied in detail in the sequel. As in the remaining of this paper, the cavity aspect ratio is fixed to Ly = 3Lx, the interdipole spacing is d0 = 4a, and the dimensionless dipole strength k0a = 0.1.

-π -π/2 0 π/2 π 0.8 1 1.2 1.4 1.6 qd 0 ω pol qτ /ω 0 ω ph q /ω 0 ω dp q /ω 0 UP LP L x /a = 7 (a) (b) -π -π/2 0 π/
emphasize that even if we model the system using a fully quantum theory, the basic ingredients, namely dipolar excitations interacting with confined electromagnetic modes, are of purely classical nature and no quantum effects are considered here. Hence, the system could also be modeled using a classical formalism relying on Maxwell's equations, which is however, in our view, less transparent and less straightforward to implement [START_REF] Markel | Propagation of surface plasmons in ordered and disordered chains of metal nanospheres[END_REF][START_REF] Mann | Tunable pseudomagnetic fields for polaritons in strained metasurfaces[END_REF][START_REF] Koenderink | Complex response and polariton-like dispersion splitting in periodic metal nanoparticle chains[END_REF][START_REF] Allard | Quantum theory of plasmon polaritons in chains of metallic nanoparticles: From near-to far-field coupling regime[END_REF][START_REF] Fernique | Plasmons in two-dimensional lattices of near-field coupled nanoparticles[END_REF][START_REF] Mann | Topological transitions in arrays of dipoles coupled to a cavity waveguide[END_REF].

III. EIGENSPECTRUM ANALYSIS AND LOCALIZATION PROPERTIES

A. Ordered chain in the thermodynamic limit

To gain insight on the model and for clarity, we begin by summarizing and further discussing results recently obtained in Ref. [START_REF] Downing | Polaritonic Tamm states induced by cavity photons[END_REF] for the simple case of an ordered chain in the thermodynamic limit, which can be solved analytically. Within the framework developed in Sec. II, this implies to consider the disorder parameters W = 0 and ∆ = 0, such that ω i = ω 0 and d i = d 0 . By going to the limit N -→ ∞, we can use periodic boundary conditions for both the dipole chain and the cavity in the z direction. This allows one to move into wavevector space and to Fourier diagonalize the periodic boundary condition version of the Hamiltonian [START_REF] Engelhardt | Unusual dynamical properties of disordered polaritons in microcavities[END_REF]. For consistency, we give in Appendix A some details about the Fourier diagonalization, which can also be found in Ref. [START_REF] Downing | Polaritonic Tamm states induced by cavity photons[END_REF].

The result of the diagonalization procedure with periodic boundary condition is shown is Fig. 2, where we plot the polaritonic dispersion ω pol qτ [see Eq. (A6)] in the first Brillouin zone (BZ), for cavity heights L x = 7a [Fig. 2(a)] and L x = 12a [Fig. 2(b)]. The colormap represents the photonic part Ph qτ of each eigenmodes, given in Eq. (A7b), and we hereafter refer to the high (low) frequency band, labeled by the index τ = + (-), as the upper (lower) polariton UP (LP). We also plot by blue and green dashed lines the bare dipolar ω dp q and bare photonic ω ph q dispersions, given respectively in Eqs. (A2) and (A3).

We first notice that increasing the cavity height L x (and hence L y , since the aspect ratio is fixed to L y = 3L x ) leads to an increasing renormalization of the bare dispersions ω dp q and ω ph q , with an increasingly pronounced avoided crossing, typical of the strong-coupling regime. When L x becomes sufficiently large, namely for L x ω 0 πL x /ω dp q=0 k 0 L y which corresponds to L x 10.3a with the parameters used in Fig. 2, the light-matter detuning ∆ q = (ω ph qω dp q )/2 can be smaller than 0, meaning that there is a crossing between the bare dipolar and photonic dispersions at a given wavenumber q res . From this point on, we can thus identify ξ qres [see Eq. (A4)] with half of the Rabi splitting frequency. In the case of L x /a = 12 shown in Fig. 2(b), the Rabi splitting frequency Ω RS 0.12ω 0 .

We note that increasing the cavity height L x reduces the light-matter coupling ξ q in units of ω 0 , as can be seen from Eq. (A4). This can be understood by noting that the lightmatter detuning ∆ q is reduced when the cavity height increases, due to the cavity size dependence of the bare photonic dispersion ω ph q [see Eq. (A3)]. This behavior will be reversed in the limit of large cavity size L x a, where our low-frequency approximation (10) breaks down, and where one recovers the regime of a dipolar chain coupled to vacuum electromagnetic modes only, studied, e.g., in Refs. [START_REF] Downing | Retardation effects on the dispersion and propagation of plasmons in metallic nanoparticle chains[END_REF][START_REF] Allard | Quantum theory of plasmon polaritons in chains of metallic nanoparticles: From near-to far-field coupling regime[END_REF].

From the colormap of Fig. 2, one can see that the hybridization between light and matter degrees of freedom is enhanced by increasing the cavity height L x , in agreement with the entrance in the strong-coupling regime. The most hybridized states are the ones near the center of the BZ, for which the detuning ∆ q is around and below 0 [see expression (A7b)]. In the following, we will refer to these highly hybridized states as polaritons. On the other hand, states that are apart from the center of the BZ, for which the detuning ∆ q is positive, remain almost unhybridized. The states in the LP branch that remain almost fully dipolar, namely with a photonic weight Ph qτ 0.1, will be referred to as dark states3 (visible in black in Fig. 2). Unlike polaritons, these dark states all have an eigenfrequency very close to that of the isolated dipole frequency ω 0 .

Finally, we indicate by a grey dotted line in Fig. 2(b) the eigenfrequency ω pol ±π/d0,-of the darkest, least coupled eigenmode of the LP branch, sitting at the edge of the BZ. This highlights the fact that due to the dipolar coupling between the emitters, which leads to the collective dispersion relation given in Eq. (A2), there are both LP eigenstates with a higher and lower eigenfrequency than that of the least coupled one. In the next subsection, by studying the localization properties of the disordered chain we will show that this at first sight trivial information becomes important, since the latter eigenstates presenting an eigenfrequency lower than ω pol ±π/d0,-will be found to be particularly robust against disorder, while the ones above the dotted line constitute what will later be called the dark state band.

B. Finite disordered chain

We now move to a numerical study of the eigenstates of the finite disordered dipole chain, considering the real space Hamiltonian [START_REF] Engelhardt | Unusual dynamical properties of disordered polaritons in microcavities[END_REF], that can be written in a 2N × 2N matrix form using the basis vector

ϕ † = (b † 1 , . . . , b † N , c † 1 , . . . , c † N ).
In this context, the dipolar and photonic parts of a given eigenstate n ∈ [1, 2N ], that is, the real space counterparts of D qτ and Ph qτ given in Eq. (A7), are given by [START_REF] Ribeiro | Multimode polariton effects on molecular energy transport and spectral fluctuations[END_REF] where the 2N -component normalized eigenvector Ψ(n) = (Ψ 1 (n), . . . , Ψ 2N (n)), the first N entries corresponding to the dipolar subspace.

D(n) = N i=1 |Ψ i (n)| 2 and Ph(n) = 1 -D(n),
To characterize the localization of an eigenstate n, we use the participation ratio (PR) defined as [START_REF] Bell | Atomic vibrations in vitreous silica[END_REF][START_REF] Thouless | Electrons in disordered systems and the theory of localization[END_REF]

PR(n) = N i=1 |Ψ i (n)| 2 2 N i=1 |Ψ i (n)| 4 , ( 14
)
where the summations are taken over the dipolar subspace only. The quantity [START_REF] Lerario | Room temperature Bloch surface wave polaritons[END_REF] gives information about the typical number of sites i occupied by an eigenstate n. In 1d systems it is then proportional to the localization length. Extended states are characterized by a PR scaling with the total number of sites N , while the PR of localized states is size-independent. We insist on the fact that in this paper, when we refer to an extended or delocalized state, we mean delocalized at the scale of the system size, namely, that the localization length is larger that the system. Indeed, no genuine finite size scaling argument has been performed to confirm that such states are still extended at the thermodynamic limit, as states that can be found, e.g., in three-dimensional disordered systems featuring an Anderson transition. A discussion about this possible ambiguity, including the scaling of the PR with the system size for N ∈ [100, 5000], is proposed in Appendix C.

Here and hereafter, we fix the interdipole spacing disorder ∆ to 0 and focus only on the effects of the frequency disorder W . Due to the fact that the individual frequencies ω i appear both on-and off-diagonal in the dipolar Hamiltonian [see Eq. ( 2)], the frequency disorder affects the spectrum more than the interdipole spacing one ∆, which is purely offdiagonal. All the results shown here are then qualitatively the same with a spacing disorder ∆ = 0. The particular case of off-diagonal disorder only which, remarkably, yields to the Dyson singularity in the density of states [START_REF] Dyson | The dynamics of a disordered linear chain[END_REF], would however not lead to important disorder-induced effects as we observe here. It is thus beyond the scope of this study.

Disorder-induced mixing of dipolar and photonic weights

To begin our study of a disordered chain of dipoles in a photonic cavity, we place ourselves in the strong-coupling regime, namely with a cavity height L x /a = 12 [see Fig. 2(b) for the ordered counterpart], and we consider a chain of N = 1000 dipoles. We compute the average of the photonic part Ph(n) and of the participation ratio PR(n) of the eigenmodes n of the system over 100 disorder realizations, for increasing frequency disorder strength W . The result is shown in Figs. 3(a) and 3(b), as a function of both the eigenfrequencies ω pol n which are ordered in ascending order, and the disorder strength W .

We first discuss the spectrum for weak disorder strengths W/ω 0 10 -2 , which could correspond typically to experimental uncertainties obtained in the fabrication of plasmonic nanoparticles, namely inhomogeneities in their sizes [START_REF] Koenderink | Experimental evidence for large dynamic effects on the plasmon dispersion of subwavelength metal nanoparticle waveguides[END_REF][START_REF] Apuzzo | Observation of near-field dipolar interactions involved in a metal nanoparticle chain waveguide[END_REF][START_REF] Mueller | Deep strong light-matter coupling in plasmonic nanoparticle crystals[END_REF], resulting in different resonance frequencies ω i [START_REF] Sönnichsen | Drastic reduction of plasmon damping in gold nanorods[END_REF][START_REF] Sönnichsen | Plasmon resonances in large noble-metal clusters[END_REF][START_REF] Berciaud | Observation of intrinsic size effects in the optical response of individual gold nanoparticles[END_REF]. In Fig. 3(a), we recognize the same behavior as with an ordered chain (see Fig. 2), namely that the bottom of the LP branch is mainly photonic, visible in orange, whereas the top is almost purely dipolar, visible in black around ω pol n /ω 0 = 1.0. Furthermore, the bottom of the UP branch is mainly dipolar, while the rest of the band goes from predominantly photonic to almost purely photonic states. By comparing these results to the PR in Fig. 3(b), we observe that mainly photonic states have a very high PR, i.e., they are delocalized along the chain, while the almost purely dipolar dark states show a lower PR, hence being localized on a small number of sites. This expected behavior shows the cavity-protection effect, namely that polaritonic states are more robust against disorder than purely dipolar states [START_REF] Houdré | Vacuum-field Rabi splitting in the presence of inhomogeneous broadening: Resolution of a homogeneous linewidth in an inhomogeneously broadened system[END_REF]. We note that mainly dipolar states with a photonic weight Ph(n) 0.25, visible in purple in Fig. 3(a), already benefit from this effect.

By now increasing disorder, we observe in Figs. 3(a) and 3(b), as expected, an increase of the bandwidth of the polariton bands. This amounts to the closing of the bandgap such that we cannot unambiguously distinguish between LP and UP branches anymore, so that the system is no longer strictly in the strong-coupling regime. Interestingly, the growing portion of the spectrum is the one containing the dark states, namely the top of the LP branch, which in this sense behaves as an effective band. This observation agrees with Ref. [START_REF] Sun | Dynamics of disordered Tavis-Cummings and Holstein-Tavis-Cummings models[END_REF], where similar behavior has been encountered. All the states outside this effective dark state band profit from the cavityprotection effect. We note that the dark state band is bounded from below by the least coupled, darkest mode, which, in the ordered chain picture, was marked as the dotted grey line in Fig. 2. As shown by the grey dashed lines in Figs. 3(a) and 3(b), the rate of expansion of the dark state band around ω 0 is W/2, that is, half the width of the rectangular distribution in which the individual dipole frequencies ω i are randomly chosen. This is due to the fact that dark states have eigenfrequencies which are almost unchanged by cavity photons, so that they remain of the order of ω i .

When this broadening of the eigenfrequencies approximately reaches the bottom of the LP branch, that is, when the disorder strength W 0.3, fully dipolar eigenstates which are totally localized with a PR(n)

1 appear in the spectrum [see the black region at the bottom right of Figs. n /ω 0 = 0.9, 1.0, and 1.06, plotted as orange, red, and blue solid lines, respectively. These eigenstates are also indicated, respectively, by an orange diamond, a red triangle, and a blue square in the ordered bandstructure of Fig. 2(b). At weak disorder strength, they correspond, respectively, to a mainly photonic polariton, a dark state, and a mainly dipolar polariton. Crucially, we observe in Fig. 3(c) that by increasing the disorder strength, the dark state at ω pol n /ω 0 = 1.0 (red line) has a photonic weight which increases linearly with the frequency disorder, at a rate around 0.25W/ω 0 , as exemplified by the black dashed line. This disorder-induced gain of photonic weight leads to the fact that no more dark states are present at the frequency where they were present without disorder, i.e., around ω 0 , since disorder hybridizes all of them to polaritons. Furthermore, the mainly photonic polariton at ω pol n /ω 0 = 0.9 (orange line) changes drastically in nature when the disorder strength is increased, becoming a mainly dipolar polariton when its frequency enters the dark state band, that is when

W 2|ω pol n -ω 0 | = 0.2.
Once in the dark state band, its photonic weight also follows a linear increase with disorder strength. Finally, the same mechanism occurs for the mainly dipolar polariton at ω pol n /ω 0 = 1.06 (blue line). The corresponding PR values in Fig. 3(d) reveal that once in the dark state band, the states show an intermediate value of the PR, here PR(n) 6, which, as opposed to the photonic weight, remains constant when the disorder strength increases. Importantly, this means that the PR of these states does not fall to 1 with strong disorder strength, in contrast to what is usually the case in 1d disordered systems. This can be understood Site i from the fact that the disorder-induced hybridization of the dark states into polaritons allows them to inherit the polariton robustness against localization, that is, the cavity-protection effect.

|Ψ i (n)| 2 ω pol n /ω 0 = 1.
10 0 |Ψ i (n)| 2 ω pol n /ω 0 = 0.9 10 -6 10 -3 10 0 |Ψ i (n)| 2 ω pol n /ω 0 = 1.0 (a) (b) (c) (d) (e) (f) (g) (h) (i) W/ω 0 = 0.01 W/ω 0 = 0.1 W/ω 0 = 1 Ph = 0.
In order to better understand the mixing of dipolar and photonic weights induced by disorder, we display in Fig. 4 the probability density |Ψ i (n)| 2 along the sites i of the chain for the same three eigenstates, for which the same color code is used. It should be noted, however, that the results of Fig. 4 correspond to a given disorder realization. Increasing values of disorder strength W/ω 0 = 0.01, 0.1, and 1 are considered from left to right panels.

In Figs. 4(a)-4(c), we present our results for the mainly photonic polariton with eigenfrequency ω pol n /ω 0 = 0.9. By increasing the disorder strength, the state goes from delocalized and mostly photonic in panels (a) and (b), to localized on multiple nonadjacent sites with a very small PR in panel (c). In the following we term these states "semilocalized", as recently proposed by Botzung et al. in Ref. [START_REF] Botzung | Dark state semilocalization of quantum emitters in a cavity[END_REF], where they unveiled similar phenomena in a disordered Tavis-Cummings model. As studied in Refs. [START_REF] Botzung | Dark state semilocalization of quantum emitters in a cavity[END_REF][START_REF] Chávez | Disorder-enhanced and disorderindependent transport with long-range hopping: Application to molecular chains in optical cavities[END_REF][START_REF] Dubail | Large random arrowhead matrices: Multifractality, semilocalization, and protected transport in disordered quantum spins coupled to a cavity[END_REF], semilocalized states can be seen as localized states, that is, an exponential peak with a sizeindependent PR, with long tails that have increased so that they are of the same order of magnitude as the original peak. In Appendix C, we discuss the nature of the semilocalized states present in our model, and notably the differences with the ones found in Refs. [START_REF] Botzung | Dark state semilocalization of quantum emitters in a cavity[END_REF][START_REF] Chávez | Disorder-enhanced and disorderindependent transport with long-range hopping: Application to molecular chains in optical cavities[END_REF][START_REF] Dubail | Large random arrowhead matrices: Multifractality, semilocalization, and protected transport in disordered quantum spins coupled to a cavity[END_REF]. In our model, these semilocal-ized states correspond to the purple region in Figs. 3(a) and 3(b), i.e., they belong to the dark state band with eigenfrequencies ω ph q=0 ω pol n ω 0 + W/2. They present a constant value of PR of the order of 10, as well as a photonic weight increasing linearly with the disorder strength.

In the second row of Fig. 4, we show the same quantities for a dark state with eigenfrequency ω pol n /ω 0 = 1.0. As one can see from panel (d), this fully dipolar state already becomes affected by a weak disorder strength W/ω 0 = 0.01, but the exponential localization occurs only at a stronger disorder strength W/ω 0 = 0.1, visible in panel (e). We note that in addition to the exponential peak, exponentially localized states also include algebraic tails with very small amplitude. By increasing further the disorder strength to W/ω 0 = 1 in panel (f), the state acquires a photonic weight and becomes semilocalized. The same mechanism is visible in the third row of Fig. 4 [panels (e)-(g)] for a mainly dipolar polariton with eigenfrequency ω pol n /ω 0 = 1.06, which, due to its polaritonic nature, remains unaffected at weak disorder strength [see Fig. 4(g)].

Influence of the light-matter coupling

Before going to transport computations, here we conclude the eigenspectrum analysis by studying the influence of the light-matter coupling on the eigenstate properties discussed previously. For that purpose, we show in Fig. 5 the photonic weight [Eq. ( 13)] and the PR [Eq. ( 14)] of the eigenmodes as a function of both the eigenfrequencies ordered in ascending order, and the cavity height L x that controls the light-matter coupling.

In the left panels, Figs. 5(a)-5(b), a weak disorder strength W/ω 0 = 0.01 is considered. As can be seen from panel (a), in the weak-coupling regime (L x /a 9) the eigenstates are almost not hybridized and the LP branch corresponds only to dark states. The corresponding PR in panel (b) shows that these dark states are already localized, especially the modes near the band edges, whereas the ones around the middle of the band present a higher PR. This behavior of sharply localized band edges is usual in all 1d Anderson-like disordered systems [START_REF] Evers | Anderson transitions[END_REF]. In the insets of Figs. 5(a) and 5(b), we show a zoom around the modes which interact the most with cavity photons, namely the ones at the center of the BZ in the ordered chain picture (around q = 0 in Fig. 2). The inset of panel (a) shows that theses modes have a higher photonic weight than the surrounding ones, being the most hybridized states. However, the inset of panel (b) demonstrates that these states have a smaller PR than the surrounding ones, hence being more localized, despite their larger photonic weight. This counterintuitive behavior highlights the particularly nontrivial interplay between localization and light-matter hybridization.

At larger cavity height (L x /a 9), we see in Figs. 5(a)-5(b) the emergence of polaritons at the bottom of the LP branch, with a high PR value. The existence of these eigenstates is permitted by the strong enough band modification due to the avoided-crossing, which allows eigenfrequencies smaller than the one of the darkest mode [which was marked as the grey dotted line in Fig. 2(b)]. 4 In contrast, the top of the q=0 , which we represent as a function of L x /a as a grey dashed line. This demonstrates that in our model, semilocalization can be understood as eigenstates being mixed to photonic states by the increase of the bandwidth induced by the frequency disorder. This permits new coupling between dipolar and photonic degrees of freedom. States with an eigenfrequency ω pol n ω ph q=0 , visible in black in Figs. 5(c)-5(d), are the fully dipolar exponentially localized states, namely, dark states subject to Anderson localization. Due to their eigenfrequencies being lower than the ones of the photonic band, they cannot benefit from the disorder-induced mixing between dipolar and photonic degrees of freedom. On the other hand, states not already reached by the dark state band, with ω pol n ω 0 + W/2, remain mainly photonic and delocalized [see the orange/yellow originates solely from the 1/r 3 quasistatic dipole-dipole coupling that leads to a dipolar band which is not flat. Indeed, a flat dipolar dispersion, as considered in usual strong-coupling toy models, would lead to eigenfrequencies that are smaller than the most uncoupled one for any nonzero value of the light-matter coupling. regions in Figs. 5(c)-5(d)], confirming the previous analysis of Figs. 3(a)-3(b). Finally, we note that with such a large value of the disorder strength (W/ω 0 = 1), we are able to see semilocalized states already for a small cavity height L x /a 7, a value from which we considered the system being in the weakcoupling regime without disorder. The transition from dark to semilocalized states when increasing the disorder strength is smooth, but with the above explanation, we understand that for a given value of W , semilocalized states would be present only if 2|ω ph q=0ω 0 | W , the light-matter coupling competing with disorder [START_REF] Botzung | Dark state semilocalization of quantum emitters in a cavity[END_REF].

To conclude this section, by studying the eigenspectrum of our system, we have shown that the bandwidth increase led by frequency disorder tends to incorporate photonic states to the dark state band, leading to states that increase their photonic weight and become semilocalized. Almost fully dipolar dark states, such as the ones shown in panels (d)-(f) of Fig. 4, can thus undergo first a transition from being delocalized along the chain to being exponentially localized on a few sites, and then another transition to being semilocalized and turning into hybridized polaritonic states. As we will show in the next section by performing transport simulations in a driven-dissipative scenario, the latter hybridized states inherit long-range polaritonic transport properties. This remarkable phenomenon thus leads to a disorder-enhanced propagation along the chain, that is, to the decrease and then the re-increase of the long-range transport properties for increasing disorder.

IV. DRIVEN-DISSIPATIVE TRANSPORT SCENARIO

We now move to transport simulations in order to elucidate the fate of the cavity-protection effect, as well as the disorderinduced mixing between dipolar and photonic degrees of freedom and semilocalization in the propagation characteristics along the chain. Especially, we here take into account inevitable losses both in the dipoles and in the cavity mirrors. For that purpose, we consider a driven-dissipative transport scenario by adding to the polaritonic Hamiltonian H given in Eq. ( 10) a driving term

H drive (t) = Ω R f (t) b † 1 + b 1 , (15)
which corresponds to a transversely polarized electric field with amplitude E 0 acting on the first dipole, with the Rabi frequency Ω R = E 0 Q 2 /2M ω 0 , and where f (t) is a timedependent function that depends on the driving frequency ω d .

In the following, we consider a drive f (t) = e iω d t , corresponding to the first site of the chain being continuously illuminated by a monochromatic electric field. We assume that the propagation dynamics can be approximated by the Lindblad master equation for the density matrix

ρ = i ρ, H + H drive (t) - i γ i 2 b † i b i , ρ -2b i ρb † i - nz κ nz 2 c † nz c nz , ρ -2c nz ρc † nz , (16)
where we consider two different phenomenological Markovian baths in order to take into account both the damping rates γ i of the dipole excitations, typically coming from radiative and Ohmic losses, and the damping rates κ nz of the photons, arising from the imperfect cavity mirrors. In the remainder of our paper, we assume that the dipolar losses are independent of the dipole site, such that γ i = γ, and that the cavity losses are independent of the photon mode, i.e., κ nz = κ. To study the transport properties along the chain of dipoles, we introduce the dimensionless dipole moment

p i = b i + b † i
bared by the dipole i. As detailed in Appendix D, we obtain from the master equation ( 16) the steady-state solution for the amplitudes |p i |.

In the following, we fix the cavity loss to κ/ω 0 = 0.001. This is motivated by the fact that all the qualitative propagation features found here are independent of the value of κ, as long as it remains small enough, as discussed in Appendix E. The value of the dipolar damping rate γ, however, is of crucial importance and in order to illustrate its impact, we consider two different cases. The first one, γ/ω 0 = 0.001, can be achieved experimentally in platforms with small or highly controllable losses such as, e.g., microwave resonators or dielectric and SiC nanoparticles. The second one, γ/ω 0 = 0.02 represents the case of more lossy dipoles, which is naturally achieved in, e.g., nanoplasmonic setups where Ohmic losses are significant.

A. Cavity-enhanced transport in an ordered chain

In order to understand the underlying transport mechanisms of the system, we begin our study by considering the ordered case (W/ω 0 = 0), which, due to the presence of polaritonic excitations, already features interesting transport properties. Indeed, cavity photon excitations having an intrinsically collective and delocalized character, they are naturally propagating at longer distances. Moreover, the light-matter coupling Hamiltonian (6) acts as an effective long-range coupling between the dipoles. Being in the Coulomb gauge, it can be seen as containing the retardation effects of the dipole-dipole interaction. In the limit of an infinite cavity, i.e., in vacuum, it thus amounts to consider the usual 1/r long-range dipoledipole coupling term. Such 1/r transport has been studied in the past in ordered [START_REF] Downing | Retardation effects on the dispersion and propagation of plasmons in metallic nanoparticle chains[END_REF][START_REF] Brandstetter-Kunc | Nonradiative limitations to plasmon propagation in chains of metallic nanoparticles[END_REF][START_REF] Fung | Analytical properties of the plasmon decay profile in a periodic metalnanoparticle chain[END_REF] and disordered [START_REF] Markel | Propagation of surface plasmons in ordered and disordered chains of metal nanospheres[END_REF] systems, showing the above-mentioned effect of long-range transport enhancement by the light-matter coupling. In a finite cavity allowing for the strong-coupling regime, it has been shown in the case of a single mode Tavis-Cummings model [START_REF] Botzung | Dark state semilocalization of quantum emitters in a cavity[END_REF][START_REF] Chávez | Disorder-enhanced and disorderindependent transport with long-range hopping: Application to molecular chains in optical cavities[END_REF][START_REF] Dubail | Large random arrowhead matrices: Multifractality, semilocalization, and protected transport in disordered quantum spins coupled to a cavity[END_REF] that under some approximations, the light-matter coupling can be expressed as a distance-independent hopping term.

Here, to understand the effects of the strong-coupling regime, we begin our analysis by comparing transport properties in systems with different cavity heights. 5 To this end, in Fig. 6 we show on a log-log scale the steady-state amplitude of the dipole moment |p i | along the sites i of a chain of N = 2500 dipoles, for cavity heights L x /a = 7 (red lines) and L x /a = 12 (blue lines). The amplitudes are given in units of Ω R /ω 0 , since |p i | ∝ Ω R /ω 0 . Two different dipole losses are considered, γ/ω 0 = 0.001 (solid lines) and γ/ω 0 = 0.02 (dotted lines), and the result corresponds to the driving of a dark state at frequency ω d /ω 0 = 1.0. In the figure, the dashed and dotted grey lines show the result in the case of a quasistatic dipole chain, without any light-matter coupling, and only with nearest-neighbor quasistatic dipole interaction. Such decay has been computed analytically in Ref. [START_REF] Brandstetter-Kunc | Nonradiative limitations to plasmon propagation in chains of metallic nanoparticles[END_REF] for the root-mean-square of the dipole moment. Here, for the amplitude of the dipole moment it translates into

|p i | = 2 π Ω R Ω 0   1 + γ 4Ω 0 2 - γ 4Ω 0   i , (17)
such that it follows an exponential decay |p i | ∝ e -id0/ζ , with i being the dipole site and ζ = d 0 /arcsinh (γ/4Ω 0 ) the propagation length. Remarkably, we observe from Fig. 6 that such an exponential decay is in very good agreement with the propagation of both the weak (red lines) and strong (blue lines) coupling regimes at short distances (i 100) for small losses (solid lines), and at very short distances (i 10) for larger losses (dotted lines). Therefore, strong coupling does not modify the short-range propagation properties of the dark states, which coupling of a dipole chain to vacuum electromagnetic modes, recovered here for an infinite cavity.

are entirely attributable to the nearest-neighbor quasistatic dipole-dipole coupling. The effect of the strong light-matter coupling becomes visible only at larger distances and takes the form of a second, less steep exponential decay, occurring after sharp oscillations that signal a change of regime.

After the exponential laws, the propagation follows an algebraic decay with slope ∼ 1/i 3 , visible as straight lines on the log-log scale of Fig. 6. Such power law results from the 1/r 3 quasistatic dipole-dipole coupling of the dipolar Hamiltonian [START_REF] Feist | Extraordinary exciton conductance induced by strong coupling[END_REF]. It is hence also visible in the quasistatic case, without light-matter interaction, when going beyond the nearestneighbor approximation performed in Eq. [START_REF] Lerario | High-speed flow of interacting organic polaritons[END_REF]. Moreover, this behavior has also been noticed in dipole chains coupled to vacuum electromagnetic modes [START_REF] Markel | Propagation of surface plasmons in ordered and disordered chains of metal nanospheres[END_REF][START_REF] Brandstetter-Kunc | Nonradiative limitations to plasmon propagation in chains of metallic nanoparticles[END_REF][START_REF] Fung | Analytical properties of the plasmon decay profile in a periodic metalnanoparticle chain[END_REF] and it also appears without quasistatic dipole-dipole but with nonzero lightmatter coupling, the photons inducing an effective all-to-all dipolar coupling. Importantly, such an algebraic decay leads the cavity-enhanced transport induced by strong coupling, that is, the second exponential regime, to be effective only for specific system sizes. At too large distances, namely here with dipole losses γ/ω 0 = 0.02 for i 300 (see the crossing between the blue and red dotted lines), the algebraic decay becomes dominant and the propagation is again similar for systems in the weak-and strong-coupling regimes. Now that we have studied the transport properties of a dark state in Fig. 6, we move to the propagation of polaritons. To that purpose, in Fig. 7 we fix the cavity height to L x /a = 12 and we compare the propagation along the chain between driving frequencies ω d /ω 0 = 0.9, 1.0, and 1.06, which correspond, respectively, to a mainly photonic polariton, a dark state, and a mainly dipolar polariton. Note that states with the same eigenfrequencies were already studied in the previous section, see Figs. 3 and4, where the same color code was used.

By looking at Fig. 7(a), where small dipole losses γ/ω 0 = 0.001 are considered, it is clear that the amplitude of the dipole moment on the first site, |p 1 |, is much smaller for polaritons (orange and blue lines) than for the dark state (red line). This can be understood from the fact that the drive [START_REF] Orgiu | Conductivity in organic semiconductors hybridized with the vacuum field[END_REF], as well as the amplitude that we evaluate, are dipole-related quantities. The dark eigenmodes can thus be more easily excited by the drive than the polaritons, and the maximum amplitude resulting from a drive at a given frequency is then related to the dipolar part of the associated eigenstate. The first exponential is very steep for the polaritons, especially for the one which is mainly photonic, in orange, leading to an inefficient short-range transport. Indeed, the first exponential regime arising solely from the quasistatic dipole-dipole interaction, it is less efficient for a state whose photonic part is important, and for the polaritons it occurs on the first three sites only. It is quickly supplanted by the second, photon-induced exponential which is less and less steep when the photonic part of the driven state increases (from red to blue to orange lines). Such second exponential decay allows for the cavityenhanced transport, and the long-range propagation of the polaritons becomes better than the one of the dark state from i 1000. We note that the slight increase between the two exponentials arises from the change of regime. By increasing the dipole losses to γ/ω 0 = 0.02 in Fig. 7(b), we observe that, crucially, it is mainly the first, dipole-induced exponential and thus the dark state propagation that suffers from larger losses. Thereby, it leads to a strong dominance of the polaritons at long distances, thanks to the second, photoninduced exponential regime. The mainly photonic polariton driven at ω d /ω 0 = 0.9, in orange, here leads to a dipole moment amplitude at the end of the chain, |p N |, which is 10 4 times larger than the one of the dark state in red. Due to the algebraic decay, however, the propagation of the dark state catches up the ones of the polaritons at a long enough distance, leading, importantly, to a size-dependency of the cavity-enhanced transport. For example, the mainly dipolar polariton (blue curve) shows a cavity-enhanced transport only from sites i ∼ 50 to i ∼ 1000.

To conclude this subsection, we have seen that the shortrange transport is dominated by the nearest-neighbor quasistatic dipole-dipole coupling, that leads to a first exponential decay which provides efficient dark state propagation, but very poor polaritonic one. The medium and long-range transport are, on the other hand, highly influenced by the strong light-matter coupling, which leads to the appearance of a second exponential decay regime showing the effect of cavityenhanced transport. The steepness of this photon-induced de- cay is flattened when the photonic part of the driven eigenstate increases, leading polaritons to propagate efficiently at long distances. At larger distances, however, an algebraic decay independent of the nature of the eigenstate cancels this transport enhancement. Crucially, the increase of dipole losses lowers essentially the short-range propagation, and has less impact on the photon-induced second exponential decay regime. This allows polaritons to better dominate medium and long-range transport when highly lossy dipoles are considered.

In the next subsection, we study the propagation along the chain in the presence of disorder. We show that the same transport mechanisms are present, with however, the crucial addition of the disorder-induced mixing between dipolar and photonic weights, which we unveiled in Sec. III. Indeed, the dark states, that have been turned into semilocalized polaritonic states by disorder, will inherit the polaritonic propagation properties, that is, cavity-enhanced transport, that we discussed in this subsection.

B. Disorder-enhanced transport

The propagation along a disordered chain is presented in Fig. 8, where the cavity height L x /a = 12. We show on a log-log scale the steady-state amplitude of the dipole moment |p i | along the sites i of the chain in units of the reduced Rabi frequency, resulting from a drive at frequency ω d /ω 0 = 1.0, for increasing disorder strength W . In Fig. 8(a), small dipole losses γ/ω 0 = 0.001 are considered. In such a case, we observe that increasing the disorder strength from W/ω 0 = 0 (blue line) to W/ω 0 = 0.01 (orange line) and W/ω 0 = 0.1 (green line) suppresses more and more the transport, especially affecting the first exponential decay and the amplitude on the first site. Increasing further the disorder to W/ω 0 = 1 (black line), we observe an even weaker short-range transport, but a stronger long-range one. Crucially, this long-range propagation arises from a flatter second exponential decay. The steepness of the latter decay being related to the photonic weight of the driven eigenstate, this corroborates the disorder-induced gain of photonic weight of the dark states around ω 0 that we unveiled in Fig. 3(c). This demonstrates the link between the observed disorder-enhanced transport and the disorder-induced mixing between dark and photonic states, such that at strong disorder strength (black line), the state driven at ω d /ω 0 = 1.0 is not anymore an almost uncoupled dark state but an hybridized polariton. It thus propagates very similarly to polaritons in the absence of disorder (cf. the orange and blue lines in Fig. 7). We note that it also supports the observations made when looking at the localization of the eigenstates for increasing value of W in Fig. 4, where panels (d), (e), and (f) correspond, respectively, to the orange, green, and black curves of Fig. 8.

When increasing the dipoles losses to γ/ω 0 = 0.02 in Fig. 8(b), the dark state transport is strongly reduced while the polaritonic one remains almost unaffected. This enables an increasingly effective disorder-enhanced transport when the considered dipoles are lossy. Here, it leads to a long-range propagation that is stronger with disorder than without, up to four orders of magnitudes around the sites i ∼ 300-400.

To highlight the effect of disorder on long-range transport, we show in Fig. 9 the evolution of the steady-state amplitude of the last dipole, |p N |, of a chain of N = 500 dipoles when the frequency disorder is increased. We normalize the value of the amplitude |p N | by the one at the end of an ordered chain, |p N (W = 0)|. Small dipole losses γ/ω 0 = 0.001 (solid line) as well as larger ones γ/ω 0 = 0.02 (dotted line) are considered, and the driving frequency is ω d /ω 0 = 1.0, just as in Fig. 8. While the long-range transport is only slightly affected at small disorder strength, it becomes suppressed when W increases, in agreement with the eigenspectrum analysis of Sec. III, where exponentially localized states were observed (see Fig. 4). Then, increasing further the disorder strength leads to the rise of the amplitude on the last dipole site, that is, to the disorder-enhanced transport regime. With small dipole losses (solid line), this increase does not exceed the value without disorder, remaining about 10 3 times smaller. But crucially, by considering more lossy dipoles (dotted line), an amplitude of the last dipole moment about 10 3 times larger than without disorder can be observed with a disorder strength W/ω 0 = 1.5. We note that from W/ω 0 = 1, the curves become essentially flat, indicating a regime where the amplitude is no longer affected by disorder. This is in line with what has been termed disorder-independent transport in Ref. [START_REF] Chávez | Disorder-enhanced and disorderindependent transport with long-range hopping: Application to molecular chains in optical cavities[END_REF], where a study comparable to ours was carried out for a disordered Tavis-Cummings model in a two-terminal, lossless transmission setup. In our model, we explain this regime from the fact that once a state gained enough photonic weight to inherit an almost flat, efficient polaritonic transport up to the end of the chain, increasing further the disorder has no effect.

However, as well as the cavity-enhanced transport unveiled in Fig. 7, in our model, the size of the chain N is critical to observe the disorder-enhanced and disorder-independent transport regimes. This is due to the complexity of the propagation which is made of different regimes of decay, and especially due to the algebraic 1/i 3 tail appearing at very large distances. For example, in Fig. 8(b), the highly disordered case with W/ω 0 = 1 (black line) shows an enhanced propagation only between the sites i 100 and i 2000. To improve the transport at longer distances, a larger disorder strength should be considered, such that the photonic weight of the driven eigenstate further increases, leading to a flatter second exponential decay which can still dominate the algebraic decay at the end of the chain. Remarkably, however, since the disorder parameter W/ω 0 cannot be arbitrarily increased [see the discussion after Eq. ( 3)], this points to the fact that the disorder-enhancement cannot be reached for too long chains by tuning only the frequency disorder.

We illustrate this effect of size-dependence in Fig. 10, by computing the steady-state amplitude of the dipole moment normalized by the one without disorder, as a function of both the site i along the chain and the disorder strength W . The dipole losses γ/ω 0 = 0.02 and the driving frequency ω d /ω 0 = 1.0, as in Fig. 8(b) and in Fig. 9 (dotted line). In such density plot, one can distinguish between the different regimes of transport induced by disorder, according to the size of the chain, i.e., to the site i. The white phase at the bottom of the plot shows that small disorder has no effect even on large chains, while the blue phase shows the Anderson localization regime, with a reduced transport. The red phase, at the top right of the plot, indicates the disorder-enhanced transport, especially where the amplitude of the dipole moment is larger with disorder than without. At short distances, such red phase is absent, showing that no disorder-enhanced transport is achievable for small chains with the parameters considered here. For larger distances, one remarks that the disorder required to enter in the red phase goes linearly with the distance i, implying that a larger chain needs a larger disorder strength.

Furthermore, it is interesting to note that decreasing the dipolar coupling Ω 0 can enhance the disorder-induced effects (not shown), allowing for a disorder-enhanced transport for smaller disorder strength W , and operating at shorter distances. In view of our model, we understand this by the fact that it would lead to a suppression of the first exponential decay, but would have less impact on the second one, the latter being mediated by cavity photons. The relative dominance of the polaritonic long-range transport over the dark state one would thus be larger, in a similar fashion as with the increase of dipole losses. At specific distances, a very small dipolar coupling can also lead to a revival of Anderson localization appearing after the disorder-enhanced and disorder-independent regimes, similarly to what was observed in Ref. [START_REF] Chávez | Disorder-enhanced and disorderindependent transport with long-range hopping: Application to molecular chains in optical cavities[END_REF]. However, decreasing the dipolar coupling, i.e., increasing the spacing d/a between the dipoles, would also go along with a drastic overall reduction of the propagation [START_REF] Brandstetter-Kunc | Nonradiative limitations to plasmon propagation in chains of metallic nanoparticles[END_REF].

Finally, to study the effect of driving frequencies ω d /ω 0 = 1.0, we show on Fig. 11 the steady-state amplitudes of the dipole moments in a density plot, as a function of both the driving frequency and the site i along a chain of N = 1000 dipoles. On panels (a), (b), and (c), we compare the propagation in a chain with disorder strengths W/ω 0 = 0, W/ω 0 = 0.1, and W/ω 0 = 1, respectively. We consider lossy dipoles, with γ/ω 0 = 0.02, for the purpose of maximizing the dominance of the long-range polaritonic transport over the dark state one. In Fig. 11(a), without disorder, a bright yellow spot around ω d /ω 0 = 1.0 demonstrates that the maximal amplitude is generated at short range by driving a dark state, supporting the fact that these states are more easily driven than polaritons. They present however inefficient long-range propagation (see the purple and black spots). At the bottom of the spectrum, a clear threshold can be seen around ω d /ω 0 = 0.85, corresponding to the minimum of the LP branch. The driven states below this threshold are not anymore eigenmodes of the system and hence present poor transport properties. Finally, the orange strips above and below the bright spot around ω d /ω 0 = 1.0 correspond to hybridized polaritons, with poor short-range propagation but very good long-range one, as already observed in Fig. 7.

By adding disorder, in Fig. 11(b), we clearly see the expansion of the effective dark state band around ω d /ω 0 = 1.0. This expansion goes along with a decrease of the transport efficiency, with the disappearance of the bright yellow spot at short range and the enlargement of the black spot at longer range. Polaritons outside the dark state band (orange strips), on the other hand, remain not affected by disorder, displaying the effect of cavity protection of their transport against disorder. When further increasing the disorder strength, in Fig. 11(c), all the frequencies up to ω d /ω 0 = 1.5 belong to the dark state band. All these now hybridized and semilocalized states have approximately the same propagation characteristics. This, crucially, leads to an improved long-range transport of the dark states around ω d /ω 0 = 1.0, which can be readily seen from the complete disappearance of the black spot. We note that it also leads to a slightly enhanced shortrange transport of the mainly photonic polaritons, the latter gaining a dipolar weight from the disorder-induced hybridization, allowing them to be more easily driven. Their longrange transport is however reduced by disorder. In agreement with what was observed in Sec. III through the fully dipolar and exponentially localized states visible as the black regions in Figs. 3(a) and 3(b), the eigenstates driven with a frequency essentially smaller than the lowest bare photonic state ω ph q=0 /ω 0 0.85 do not inherit polaritonic propagation features, the light-matter coupling being not strong enough to let them be coupled to photons by the disorder.

V. CONCLUSIONS

We provided a study of the interplay between light-matter coupling, dipolar coupling, and disorder, through the analysis of both the localization properties and the transport along a chain of disordered dipoles placed inside an optical cavity. We considered a multimode light-matter coupling Hamiltonian which goes beyond the typical single-mode coupling widely used in quantum optics, as well as an all-to-all 1/r 3 quasistatic dipole-dipole coupling. Crucially, losses in both the dipoles and the mirrors are taken into account. The disorder was considered on the individual resonance frequencies ω i of the dipoles.

In our multimode model, all the eigenstates are hybridized by the strong light-matter coupling. However, we can distinguish the dark states, with eigenfrequencies around the av- erage bare dipole resonance frequency ω 0 and almost zero photonic weight, and the polaritons, which have an eigenfrequency highly renormalized by the coupling and a significant photonic weight. We showed that at weak and medium disorder strengths, dark states suffer from Anderson localization, becoming exponentially localized on a few sites of the chain. On the other hand, the polaritons are more robust against frequency disorder, and show a cavity-protection effect [START_REF] Houdré | Vacuum-field Rabi splitting in the presence of inhomogeneous broadening: Resolution of a homogeneous linewidth in an inhomogeneously broadened system[END_REF][START_REF] Michetti | Polariton states in disordered organic microcavities[END_REF]. Notably, the only states impacted by disorder are the ones located in the same frequency window as the possible disordered individual dipole frequencies ω i , namely around ω 0 ± W/2. Increasing the disorder strength leads this frequency window to act as an effective dark state band, whose bandwidth increases at a rate W/2. At large disorder strength, we have shown that the states located both in this dark state band and with an eigenfrequency larger than the lowest bare photonic mode ω ph q=0 are subject to a disorder-induced hybridization with photons, allowing new coupling between matter and photonic degrees of freedom. Indeed, the photonic part of these states increases linearly with the disorder strength and such states become semilocalized [START_REF] Botzung | Dark state semilocalization of quantum emitters in a cavity[END_REF], i.e., localized in multiple nonadjacent sites. The eigenfrequencies smaller than ω ph q=0 , on the other hand, cannot be hybridized through disorder, and suffer from the usual Anderson exponential localization.

By studying the transport of an excitation along the chain in a driven-dissipative scenario, we have shown that the strong light-matter coupling regime leads to a second exponential decay in the propagation, following a first one which is solely due to nearest-neighbor quasistatic dipole-dipole interactions. The larger the photonic part of the driven state is, the poorer the short-range propagation is, but the more flat this second exponential regime is. At large frequency disorder, the dark states, which have become polaritons through the abovementioned mechanism of disorder-induced hybridization, inherit polaritonic transport properties, namely a poor shortrange transport but an efficient long-range one. Moreover, we have found that this disorder-enhanced transport is not efficient over too long distances, or at least requires very large disorder strength. Indeed, the amount of disorder required to enter into such a regime increases linearly with the system size.

Crucially, we have shown that increasing the dipole losses impacts essentially the short-range transport, mediated by the quasistatic dipole-dipole coupling, which consequently increases the dominance of the polaritons at longer range. In this way, we have found that by considering highly lossy dipoles, a large frequency disorder can increase the long-range transport up to a factor 10 4 as compared to the case without disorder, see, e.g., around the site i = 400 in Fig. 8(b). In terms of the power radiated by the dipoles, it translates into an increase of up to 8 orders of magnitude [see Eq. (D1)]. This drastic increase of long-range transport could be useful especially for setups using plasmonic nanoparticle chains, where experiments have witnessed short transport distances [36-38, 40, 41, 60, 61], notably due to large Ohmic losses. Indeed, when considering plasmonic nanoparticles, the dimensionless dipole strength k 0 a = 0.1 used throughout this work corresponds to nanoparticles with, e.g., radius a = 7.6 nm and average resonance frequencies ω 0 = 2.6 eV/ . Considering a driving field intensity I 0 = 1 MW/cm 2 , corresponding to, e.g., a laser of power P 0 = 7.1 mW focused on a spot of size πλ 2 0 0.71 µm 2 , the power radiated in the far field by a plasmonic nanoparticle with decay rate γ/ω 0 = 0.02 around the site i = 400 [see Fig. 8(b)] can then be increased from P 400 (W/ω 0 = 0) 10 -24 W to P 400 (W/ω 0 = 1) 10 -16 W. While the ordered case leads to a very small power which is clearly experimentally unreachable, a setup with subfemtowatt sensitivity should detect the radiated power at the end of a highly disordered plasmonic chain.

Our model, which allowed to study the particularly nontrivial interplay between disorder and light-matter coupling, opens the way to study more complex systems which could be of larger dimensionality, or feature topological properties [START_REF] Mann | Manipulating type-I and type-II Dirac polaritons in cavityembedded honeycomb metasurfaces[END_REF][START_REF] Mann | Tunable pseudomagnetic fields for polaritons in strained metasurfaces[END_REF][START_REF] Downing | Topological phases of polaritons in a cavity waveguide[END_REF][START_REF] Mann | Topological transitions in arrays of dipoles coupled to a cavity waveguide[END_REF], for which the interplay with disorder and a coupling beyond nearest neighbor has recently presented interesting phenomena [START_REF] Pérez-González | Interplay between long-range hopping and disorder in topological systems[END_REF]. Moreover, there is a growing interest in topological photonics [START_REF] Ozawa | Topological photonics[END_REF], which one of the key chal-lenges is to combine disorder-robust topological edge states with strong light-matter coupling. A detailed understanding of the interplay between disorder and strong light-matter coupling, allowed by our model, would then be of crucial importance to this field.

Note added. Recently, Ref. [START_REF] Engelhardt | Polarition localization and spectroscopic properties of disordered quantum emitters in spatiallyextended microcavities[END_REF], which considers a similar system as ours, in particular a multimode photonic cavity, appeared on the arXiv. 12. Polaritonic band structure of an ordered chain in the thermodynamic limit, in units of the bare frequency ω0, and as a function of the reduced wavenumber qd0 in the first Brillouin zone. In orange, we show the dispersion used in the main text, which is obtained using the RWA and given by Eq. (A6), and we compare it to the one obtained without using the approximation, given by Eq. (B3) and shown in green. The considered cavity height is Lx/a = 12.

where the dipolar Hamiltonian with counter-rotating terms reads [58]

H CR dp = N i=1 ω i b † i b i + N -1 i=1 N j=i+1 Ω ij b i + b † i b j + b † j .
(B2) In the limit of an ordered chain, similarly to what has been done in Appendix A, one can rewrite the Hamiltonian (B1) in Fourier space using periodic boundary conditions and diagonalize it using a Hopfield-Bogoliubov transformation. The diagonalization leads to the polaritonic dispersion

w pol qτ = Γ 2 q + τ Γ 4 q -(w dp q ) 2 (ω ph q ) 2 -4 ξ 2 q ω 0 ω ph q , (
B3) where we have introduced the quantity

Γ 2 q = 1 2 ω ph q 2 + w dp q 2 , (B4)
with the bare photonic dispersion given in Eq. (A3), and where

w dp q = ω 0 1 + 2f q Ω 0 ω 0 (B5)
is the dispersion of the dipolar Hamiltonian with counterrotating terms (B2). Details on this diagonalization procedure can be found in the supplementary material of Ref. [START_REF] Downing | Polaritonic Tamm states induced by cavity photons[END_REF], with the only difference that here we neglect the diamagnetic A 2term in the light-matter Hamiltonian [START_REF] Botzung | Dark state semilocalization of quantum emitters in a cavity[END_REF]. In Fig. 12, we show a comparison between the polaritonic dispersions obtained with (orange line) and without (green line) the RWA, given respectively by Eqs. (A6) and (B3). We observe that the two band structures overlap almost perfectly, with very slight differences noticeable only around the center Site i

|p i | / (Ω R /ω 0 )
With RWA Without RWA W/ω 0 = 0 FIG. 13. Steady-state amplitude of the dipole moment |pi| in units of the reduced Rabi frequency ΩR/ω0 along the sites i of an ordered chain. The results are shown within (orange) and without (green) the RWA. Two different frequencies are considered: ω d /ω0 = 1.0, which drives a dark state is shown as solid lines, and ω d /ω0 = 0.9, which drives a polariton with a large photonic part is shown as dotted lines. The cavity height Lx/a = 12, and the dipole and mirror losses are respectively γ/ω0 = 0.001 and κ/ω0 = 0.001.

of the BZ, where the light-matter coupling is maximal. This is in agreement with the fact that our model does not allow the ultra-strong coupling (USC) regime to be reached, since one enters into the latter regime when the effects of the counterrotating terms become sizable [START_REF] Forn-Díaz | Ultrastrong coupling regimes of light-matter interaction[END_REF], which is clearly not the case here. The fact that we remain in the strong-coupling regime only also justifies neglecting the diamagnetic A 2 -term in the light-matter coupling Hamiltonian [START_REF] Botzung | Dark state semilocalization of quantum emitters in a cavity[END_REF]. As a matter of fact, perceptible differences in the eigenspectrum caused by the absence of such diamagnetic term were found in the context of the ultra-strong or deep-strong coupling regime only [START_REF] Stefano | Resolution of gauge ambiguities in ultrastrong-coupling cavity quantum electrodynamics[END_REF][START_REF] Kockum | Ultrastrong coupling between light and matter[END_REF].

However, in the case of driven-dissipative transport computations, special attention has to be taken when using the RWA. Indeed, it has been shown [START_REF] Agarwal | Rotating-wave approximation and spontaneous emission[END_REF][START_REF] Agarwal | Rotating-wave approximation and spontaneous emission[END_REF][START_REF] West | On the rotating wave approximation[END_REF][START_REF] Ford | The rotating wave approximation (RWA) of quantum optics: serious defect[END_REF][START_REF] Fleming | The rotating-wave approximation: consistency and applicability from an open quantum system analysis[END_REF] that using the RWA along with a Lindblad master equation can lead to inaccurate dynamics, especially when dropping the counter-rotating terms directly in the Hamiltonian, as we have done in the main text. In Fig. 13, we thus compare our transport results obtained with (orange lines) and without (green lines) neglecting the counter-rotating terms in the Hamiltonian entering in the Lindblad master equation [START_REF] Zhong | Energy transfer between spatially separated entangled molecules[END_REF]. We show on a log-log scale the amplitude of the steady-state dipole moment |p i | in units of the reduced Rabi frequency Ω R /ω 0 , as a function of the dipole site i along the chain, when driving a dark state with ω d /ω 0 = 1.0 (solid lines), and a mainly photonic polariton with ω d /ω 0 = 0.9 (dotted lines). As can be seen from Fig. 13, the propagation is qualitatively the same within and without the RWA. The small differences arise essentially from the very slight frequency shift visible in Fig. 12, implying the fact that using a given driving frequency ω d does not drive an eigenstate with exactly the same photonic weight in the case of employing or not the RWA in the computation. 14)] as a function of the size of the chain N , for the eigenmodes associated to the eigenfrequencies ω pol n /ω0 = 0.9, 1.0, and 1.06, with increasing disorder strength from the left to the right panel. Panel (c) shows in addition the result for the eigenmodes associated to the eigenfrequencies ω pol n /ω0 = 0.7 and 1.55. Extended states show a participation ratio growing linearly with N , while localized states present a constant participation ratio. The black dashed line represents the maximum growth rate 2(N + 1)/3. The cavity height Lx/a = 12, and the data have been averaged over 100 disorder realizations.

Appendix C: On the nature of the semilocalized eigenstates

In the main text, we show that at large disorder strength the eigenstates of the Hamiltonian (10) can be of three very different types [cf. Figs. 3(a) and 3(b)]: (i) dipolar with a PR of order 1, (ii) hybridized with a PR of order 10, (iii) photonic with a PR of order N . By inspection of the probability density |Ψ i | 2 of these states along the sites i ∈ [1, N ] of the chain, we find that the dipolar ones are exponentially localized in the chain, with very small power-law tails, so that we coin them as localized. The photonic states are homogeneously extended along the chain, and we term them delocalized at the scale of the system size, since their localization length is larger than the size of the chain.

The hybridized states, however, are localized on multiple, noncontiguous sites [see panels (c), (f), and (i) of Fig. 4, where N = 1000], so that we coin them semilocalized, following Ref. [START_REF] Botzung | Dark state semilocalization of quantum emitters in a cavity[END_REF] were similar features were observed in a disordered Tavis-Cummings (TC) model. In the context of disordered systems, such neither localized nor delocalized states have been associated to multifractal, nonergodic extended states, which are notably a feature of the critical point of Anderson transitions [START_REF] Evers | Anderson transitions[END_REF]. Multifractal states have also been observed in disordered Floquet systems [START_REF] Roy | Multifractality without fine-tuning in a Floquet quasiperiodic chain[END_REF] as well as in long-ranged disordered systems [START_REF] Levitov | Critical Hamiltonians with long range hopping[END_REF][START_REF] Deng | Quantum Levy flights and multifractality of dipolar excitations in a random system[END_REF]. Recently, Dubail et al. [START_REF] Dubail | Large random arrowhead matrices: Multifractality, semilocalization, and protected transport in disordered quantum spins coupled to a cavity[END_REF] unveiled the multifractal nature of the semilocalized states present in a disordered TC model without short-range hopping, where the light-matter coupling acts as an effective long-range hopping term.

In this Appendix, we study the nature of the semilocalized states present in our multimode light-matter coupling model, to ascertain whether they have the same properties as in simplified single-mode models, or if their particular localization profile comes from finite-size effects. To this end, we first conduct a scaling analysis of the participation ratio with the system size N , then we perform a level statistics analysis, before we carry out a multifractal analysis through the scaling of the generalized participation ratio with the system size.

Scaling of the participation ratio

We show in Fig. 14 the scaling of the participation ratio PR(n) [defined in Eq. ( 14)] with the number of dipoles N ∈ [100, 5000], for increasing values of the disorder strength W/ω 0 = 0.01, 0.1, and 1 from the left to the right panel. We note that to minimize fluctuations, all the results have been averaged over 100 disorder realizations. Such scaling is of crucial importance to correctly conclude on the localization nature of the eigenstates. Indeed, the PR of a localized state must remain independent of the system size, while the one of an extended state should scale with the system size.

In Figs. 14(a) and 14(b), we display the results for the eigenstates corresponding to the eigenfrequencies ω pol n /ω 0 = 0.9, 1, and 1.06, i.e., the same as the ones studied in the main text, notably in Figs. 3 and4. On the one hand, at small disorder strength in panel (a), two polaritonic states (plotted by orange diamonds and blue squares) show a PR increasing linearly with the number of dipoles N , precisely following the growth rate 2(N + 1)/3 plotted as a black dashed line, which corresponds to maximally extended eigenstates in an ordered (W = 0) system. On the other hand, the PR of the dark state, plotted by red triangles, increases at small system sizes before it converges to a constant value for large enough number of dipoles N , demonstrating its localized nature despite its quite large localization length. At larger disorder strength in panel (b), the PR of both the mainly dipolar polariton (blue squares) and the dark state (red triangles) are drastically reduced and remain constant, indicating localized states, as already seen in Figs. 4(e) and 4(h). The mainly photonic polariton (orange diamonds) remains delocalized, although presenting a slightly reduced PR. Finally, in Fig. 14(c), a larger disorder strength W/ω 0 = 1 is considered and in addition to the three eigenstates discussed above which are now semilocalized, we also show the scaling of the PR of a mainly photonic polariton with eigenfrequency ω pol n /ω 0 = 1.55 (green circles) and of a fully dipolar dark state with eigenfrequency ω pol n /ω 0 = 0.7 (black crosses). In that context, one can see that all of the three previously discussed eigenstates show a PR that, after increasing at small sizes, converges to a constant value at larger sizes, confirming the fact that such semilocalized states are not extended. Moreover, the mainly photonic polariton (green circles) still has a growing PR, showing that such a state having an eigenfrequency ω pol n > ω 0 + W/2 remains extended. The eigenstates with frequencies ω pol n < ω ph q=0 newly allowed by the large value of the disorder strength are totally localized with a constant value of the PR around 1, as can be seen from the black crosses in Fig. 14(c). We note that we have also checked (not shown) that the photonic weight Ph(n) of the eigenstates is not affected when increasing the system size, and is independent of the number of dipoles N .

Level statistics analysis

We propose here an analysis of the nature of the semilocalized states observed in our model by means of level spacing statistics. To this end, we follow Ref. [START_REF] Oganesyan | Localization of interacting fermions at high temperature[END_REF] and compute the distribution of level spacing ratios

rn = min r n , 1 r n , (C1)
with

r n = s n s n-1 , (C2)
where s n = ω pol n+1ω pol n 0 is the level spacing, the eigenfrequencies ω pol n being sorted in ascending order. The average rn of such level spacing ratio over a given window of eigenstates n has been shown to display a universal value according to the level statistics [START_REF] Atas | Distribution of the ratio of consecutive level spacings in random matrix ensembles[END_REF]. Moreover, the average level spacing ratio is scale invariant at the critical point of an Anderson transition. This allows one to monitor such transition between an extended phase, with statistics related to the random matrix ensembles, and a localized phase, with eigenfunctions exponentially localized on random sites, not overlapping so that the eigenvalues follow Poisson statistics [START_REF] Oganesyan | Localization of interacting fermions at high temperature[END_REF][START_REF] Tarquini | Critical properties of the Anderson localization transition and the high-dimensional limit[END_REF][START_REF] Torres-Herrera | Level repulsion and dynamics in the finite one-dimensional Anderson model[END_REF][START_REF] Šuntajs | Spectral properties of three-dimensional Anderson model[END_REF]. In contrast to the standard level spacing distribution, the quantity (C1) has the advantage of not requiring any unfolding of the spectrum, since it is independent of the local density of states [START_REF] Oganesyan | Localization of interacting fermions at high temperature[END_REF][START_REF] Atas | Distribution of the ratio of consecutive level spacings in random matrix ensembles[END_REF].

In Fig. 15, we show such one-parameter scaling for increasing number of dipoles N , using the same parameters as in the previous subsection. The averaging of the ratio (C1) has been performed with 1/10 of the eigenstates of the spectrum centered around the middle of the dark state band in the case of an ordered chain, namely around the eigenfrequency ω pol n /ω 0 = 0.9981. The results were further averaged over 100 disorder realizations. As expected, there is no common intersection point between all the curves, demonstrating that there is no Anderson transition, as in usual 1d disordered systems. Instead, the average level spacing ratio converges to the exact value computed for the localized phase with Poisson statistics, rn PE = 2 ln 2 -1 0.39 [START_REF] Atas | Distribution of the ratio of consecutive level spacings in random matrix ensembles[END_REF], with diminishing disorder strength when the size of the system increases, as in the original 1d Anderson model [START_REF] Torres-Herrera | Level repulsion and dynamics in the finite one-dimensional Anderson model[END_REF]. At small disorder strength, the average ratio goes to 1, i.e., the value corresponding to an ordered system. Therefore, as expected in 1d [START_REF] Torres-Herrera | Level repulsion and dynamics in the finite one-dimensional Anderson model[END_REF], there is no extended phase with statistics related to random matrix ensembles.

Interestingly, when increasing further the disorder strength, we observe in Fig. 15 a rise of the average level spacing ratio, which corresponds to the semilocalization and hybridization of the dark states, as discussed in the main text. The ratio rn thus takes values between rn PE and rn = 1/2, which corresponds to semi-Poissonian statistics [START_REF] Bogomolny | Models of intermediate spectral statistics[END_REF][START_REF] Atas | Joint probability densities of level spacing ratios in random matrices[END_REF]. Such behavior is in line with the results obtained in a disordered TC model where multifractal semilocalized states were found to follow statistics that range from Poissonian to very close to semi-Poissonian, even at the thermodynamic limit [START_REF] Dubail | Large random arrowhead matrices: Multifractality, semilocalization, and protected transport in disordered quantum spins coupled to a cavity[END_REF].

Remarkably, however, in our multimode model we observe that increasing the size of the system leads the rise of the average level spacing ratio at large disorder to be flattened, so that the ratio converges to the Poisson value even at large disorder strength, as visible for the case of N = 5000 displayed by a black solid line in Fig. 15. Hence, for larger systems, a larger disorder strength is required to push the semilocalized states away from a Poisson statistics. This suggests that for a fixed, finite disorder strength W/ω 0 with the parameters chosen here, namely d/a = 4 and L x /a = 12, the semilocalized states in our model would just follow a Poisson statistics at the thermodynamic limit N → ∞, as usual localized states.

Generalized participation ratio and multifractal analysis

To ascertain whether or not the semilocalized states of our multimode model follow a multifractal structure as in singlemode models [START_REF] Dubail | Large random arrowhead matrices: Multifractality, semilocalization, and protected transport in disordered quantum spins coupled to a cavity[END_REF], here we conduct a multifractal analysis by analyzing the scaling with the system size N of the general- The exponents have been extracted according to the procedure explained in the text (see colored dots). The black dashed line represents the analytical result found for a simplified disordered Tavis-Cummings model [START_REF] Dubail | Large random arrowhead matrices: Multifractality, semilocalization, and protected transport in disordered quantum spins coupled to a cavity[END_REF]. The inset exemplifies the scaling procedure for the value q = 0.12, with the circles showing numerical data. The cavity height Lx/a = 12 and the disorder strength W/ω0 = 1.

ized participation ratio, which is defined as [START_REF] Evers | Anderson transitions[END_REF] PR

q (n) = N i=1 |Ψ i (n)| 2 q N i=1 |Ψ i (n)| 2q ∼ N →∞ N τq(n) . (C3)
While localized and delocalized eigenstates are characterized, respectively, by a multifractal exponent τ q = 0 (with q > 0) and τ q = q -1, any other behavior of τ q as a function of q implies multifractality [START_REF] Evers | Anderson transitions[END_REF]. We extract the multifractal exponent τ q from a logarithmic linear regression of the generalized participation ratio (C3) as a function of the system size N , that we average over a frequency window ω pol n /ω 0 ∈ [0.9 ; 1.4] corresponding to hybridized, semilocalized states. To minimize fluctuations, data obtained with system sizes N ∈ [10 2 ; 10 3 ] and N ∈ ]10 3 ; 3 × 10 3 ] have been further averaged over, respectively, 100 and 10 disorder realizations.

Figure 16 shows the resulting multifractal exponents extracted from different scaling procedures, considering only system sizes N ∈ [10 2 ; 10 3 ] in blue, N ∈ [10 3 ; 3 × 10 3 ] in orange, and N ∈ [3 × 10 3 ; 10 4 ] in green. In the inset we exemplify such scaling behavior for the value q = 0.12, with numerical data represented as grey circles, and the linear regressions with solid lines. We observe that when increasing the system size N , the multifractal exponents τ q are getting closer and closer to their value for localized states, namely τ q = 0. This drastically differs from the multifractal behavior found analytically at the thermodynamic limit in a simplified disordered TC model (see Eq. ( 8) of Ref. [START_REF] Dubail | Large random arrowhead matrices: Multifractality, semilocalization, and protected transport in disordered quantum spins coupled to a cavity[END_REF]) which we show as a black dashed line in Fig. 16. This suggests that the semilocalized states of our multimode model are not related to multifractal, nonergodic extended states at the thermody-namic limit. We note that in a single-mode version of our model, namely considering the Hamiltonian [START_REF] Engelhardt | Unusual dynamical properties of disordered polaritons in microcavities[END_REF] with only the photonic mode n z = 1, the exponent τ q shows the same multifractal behavior as in the simplified disordered TC model considered in Ref. [START_REF] Dubail | Large random arrowhead matrices: Multifractality, semilocalization, and protected transport in disordered quantum spins coupled to a cavity[END_REF].

To conclude, the analysis conducted in this Appendix suggests that the hybridized, semilocalized states present in our model through states localized in multiple, noncontiguous sites of the chain [see panels (c), (f), and (i) of Fig. 4 where N = 1000], behave in the thermodynamic limit as usual localized states following Poissonian statistics. This is in line with the results of our transport simulations discussed in the main text, where we show that disorder-enhanced transport, which is a consequence of the presence of hybridized, semilocalized states, vanishes for very large system sizes and is relevant at intermediate distances only (see Fig. 10). Therefore, in our multimode model, the semilocalization properties of the dark states that are hybridized through disorder fade out as the system size increases. In this Appendix, we detail our approach to find the steadystate solution of the Lindblad master equation ( 16) for the amplitudes of the dimensionless dipole moments p i = b i + b † i along the sites i of the chain. The latter quantity is directly proportional to the dipole moment at site i, p i = Q /2M ω 0 p i . The associated power radiated in the far field by a dipole at site i is then given by the classical Larmor formula [START_REF] Jackson | Classical Electrodynamics[END_REF]

P i = 2ω 4 0 3c 2 (p i ) 2 = 4πa 2 3 (k 0 a) 4 I 0 p i Ω R /ω 0 2 , ( D1
)
where I 0 = cE 0 2 /8π is the driving field intensity. In addition to p i , we introduce a corresponding momentum π i = i b ib † i , as well as the equivalent quantities for the photonic degrees of freedom, p ph nz = c nz + c † nz and π ph nz = i c nzc † nz . Using the identity Ȯ = Tr { ρO}, valid for any operator O, one obtains from the master equation ( 16) the system of coupled first order ordinary differential equations

ṗi = -ω i π i - γ i 2 p i + nz ξ inz p ph nz - j =i Ω ij π j (D2a) ṗph nz = -ω ph nz π ph nz - κ nz 2 p ph nz - i ξ inz p i (D2b) πi = ω i p i - γ i 2 π i + nz ξ inz π ph nz + j =i Ω ij p j + 2Ω R f (t) δ 1,n (D2c)
πph nz = ω ph nz p ph nz - κ nz 2 π ph nz - i ξ inz π i . (D2d)
Being interested in the stationary transport regime, we solve the above system for the steady-state solution using the complex representation p i = A i e iω d t , and we recall that we consider a drive f (t) = e iω d t . The latter continuous periodic drive leads to a vanishing time-averaged dipole moment p i t . However, its time-averaged amplitude |p i | t is nonzero. This procedure yields to steady-state solutions P

= (|p 1 |, |p 2 |, . . . , |p N |) as P = M -1 D, (D3)
where the N -dimensional driving vector is given by D = (2Ω R , 0, . . . , 0), and the N × N matrix M reads

M = ω d -i γ 2 ω d -i κ 2 1 N -M ξ (1 N -D 2 ) M t ξ × M -1 dp ω d -i γ 2 ω d -i κ 2 1 N + M ξ D 0 M t ξ -M dp . (D4)
In the above matrix equation, 1 N is the N × N identity matrix, the matrix M ξ is defined by its elements as

[M ξ ] ij = ξ ij , where i, j = 1, . . . , N , the diagonal matrices D β = Diag Ω β 1 , Ω β 2 , . . . , Ω β N , with the β-dependent function Ω β nz = ω ph nz β (ω ph nz ) 2 -ω d -i κ 2 2 , (D5)
and finally the matrix M dp is given by

M dp = M dp -M ξ D 1 M t ξ , (D6)
where M dp is defined by its elements as

[M dp ] ij = ω i δ ij + Ω ij (1 -δ ij ) . (D7)
We have checked (not shown) that such steady-state solutions |p i | are well recovered after a finite time when one directly solves the system (D2) using a numerical ODE solver.

Appendix E: Effects of cavity losses

In the main text, we fixed the cavity loss rate to κ/ω 0 = 10 -3 . Here, we motivate this choice by comparing in Fig. 17 the propagation characteristics of an ordered chain for different values of cavity loss, from κ/ω 0 = 10 -6 to κ/ω 0 = 10 -1 . The amplitude of the steady-state dipole moment |p i | in units of the reduced Rabi frequency Ω R /ω 0 is shown as a function of the dipole sites i along the chain. Two driving frequencies ω d /ω 0 = 1.0 and ω d /ω 0 = 0.9 are shown in Figs. 17(a) and 17(b), respectively. The dipole losses are fixed to γ/ω 0 = 0.001, while the cavity height is chosen as L x /a = 12.

In both panels, the qualitative behavior is substantially the same for cavity losses κ/ω 0 = 10 -6 (purple lines) and κ/ω 0 = 10 -3 (pink lines), justifying the choice made in the main text. We observe that considering larger cavity losses κ/ω 0 = 10 -2 (red lines) does not qualitatively change the results when driving a dark state in panel (a), showing the robustness of dark state propagation against cavity losses. As expected, however, the propagation of the mainly photonic polariton shown in Fig. 17(b) is more affected by the increase of cavity losses. Crucially, by increasing further the cavity loss rate to κ/ω 0 = 10 -1 (yellow lines), the second exponential regime, that is, the cavity-enhanced transport, completely disappears in Fig. 17(a), showing that the system is not anymore in the strong-coupling regime. The dark state propagation is then very well described by the exponential regime mediated solely by the nearest-neighbor quasistatic dipole-dipole coupling, given by the analytical result of Eq. ( 17) and shown as a dashed grey line. In the case of the mainly photonic polariton in Fig. 17(b), we observe similar results, namely that the long-range transport starts to become suppressed for cavity losses κ/ω 0 = 10 -1 . This transition from the strong-to the weak-coupling regime is expected, since the latter large cavity loss rate is of the order of the Rabi splitting frequency, Ω RS 0.12. We note that similar conclusions can be drawn when the dipole losses are fixed to a larger value.

Finally, we recall that the cavity quality factor Q is typically of the order of the inverse of the cavity loss rate in units of the average dipole frequency ω 0 /κ. This implies that the effects discussed in this work should be observed with cavity quality factors as low as Q ∼ 10 2 . Appendix F: Absorption spectra from transport computations

In this Appendix, we complement the discussion on transport of Sec. IV by presenting absorption spectra obtained from our driven-dissipative computations. To this end, we consider that the sum of all the steady-state dipole moment amplitudes, N i=1 |p i |, is at first order proportional to the electric field absorbed by the dipoles [START_REF] Sturges | Extreme renormalisations of dimer eigenmodes by strong light-matter coupling[END_REF].

Ordered chain

We begin with the case of an ordered chain, and show in Fig. 18 such absorption spectrum for cavity heights L x /a = 7 (blue lines) and L x /a = 12 (red lines). Two different dipole losses γ/ω 0 = 0.001 (solid lines) and γ/ω 0 = 0.02 (dotted lines) are considered.

While in weak coupling regime, L x /a = 7, the absorption spectrum is bell-shaped and centered around the dark state eigenfrequencies at ω d /ω 0 = 1.0, in the strong-coupling regime, L x /a = 12, a rich structure appears. In addition to the central peak, which is slightly higher, we observe high absorption for driving frequencies associated to hybridized polaritons (see the dispersion in Fig. 2), namely the bottom of the LP branch, starting around ω d /ω 0 = 0.85, and the UP branch, visible as a second peak starting around ω d /ω 0 = 1.05 and separated from the LP branch by a gap. The fast oscillations are peaks corresponding to the different eigenfrequencies of the system, visible one by one due to the steep slope in the dispersion in these frequency regions. We can hence precisely recover the eigenmodes from this driven-dissipative scenario. However, the absorption spectrum of Fig. 18 does not reflect the transport properties along the whole chain, the sum over the amplitudes |p i | being obviously dominated by the first few terms. It then provides information mostly about the shortrange propagation. This explains why the absorption is maximal for the dark eigenstates, and decreases when the driving frequency rises to more photonic eigenstates. Indeed, as discussed in Sec. IV, the dark states present the most efficient transport at short distances. By increasing the dipole losses, from the solid to the dotted lines, the central peak is drastically reduced and the full width at half maximum is increased, as expected. The polaritons at the bottom of the LP and UP branches are less distinct from the dark states of the central peak, the latter being the most affected by dipole losses. Crucially, however, the two bands are still distinguishable and separated by a gap despite the larger dipole losses, indicating that the system is still in the strong-coupling regime.

Disordered chain

The case of a disordered chain is shown in Fig. 19, where the sum of the amplitudes of the dipole moments along the chain for a given disorder strength, We choose a cavity height L x /a = 12 and we consider small dipole losses γ/ω 0 = 0.001. By increasing the disorder strength from W/ω 0 = 0 (blue line) to W/ω 0 = 0.01 (orange line), and to W/ω 0 = 0.1 (green line), the absorption is more and more suppressed for the driving frequencies corresponding to dark states, in a growing frequency window around ω d /ω 0 = 1.0. This is in agreement with what we observe in Sec. III, namely that the effective dark state band grows with the disorder strength W around ω 0 , while the polaritons outside the dark state band, with ω pol n ω 0 + W/2, are not affected by disorder. The polaritons hence profit from this cavity-protection effect due to their eigenfrequencies being far from ω 0 .

With a disorder strength W/ω 0 = 1 (black line in Fig. 19), the dark state band includes now frequencies up to 1.5 ω 0 . Due to the disorder-induced hybridization mechanism, eigenstates at frequencies around 1.5 ω 0 have gone from almost purely photonic polaritons to semilocalized polaritonic states with a significant dipolar part. This leads their absorption to become larger than the one without disorder, the dipolar part allowing them to be more easily excited and to acquire a better short-range propagation mediated by the quasistatic dipole-dipole coupling. Importantly, we thus observe here a disorder-enhanced short-range transport, from almost fully photonic modes acquiring a dipolar part from the increased disorder strength, and inheriting efficient short-range dipolar transport.

FIG. 2 .

 2 FIG. 2. Polaritonic dispersion ω polqτ [see Eq. (A6)] of an infinite, ordered dipole chain for cavity heights Lx/a = 7 (left panel) and Lx/a = 12 (right panel), in units of the bare frequency ω0 and as a function of the reduced wavenumber qd0 in the first Brillouin zone. The colormap represents the photonic weight Phqτ of the mode [see Eq. (A7b)], and the dotted grey line in panel (b) indicates the value of ω pol ±π/d 0 ,-. The blue and green dashed lines represent the bare dipolar dispersion ω dp q [see Eq. (A2)] and the bare photonic one ω ph q [see Eq. (A3)], respectively. The orange, red, and blue symbols point out eigenstates that will be studied in detail in the sequel. As in the remaining of this paper, the cavity aspect ratio is fixed to Ly = 3Lx, the interdipole spacing is d0 = 4a, and the dimensionless dipole strength k0a = 0.1.

FIG. 3 .

 3 FIG. 3. (a)-(b) Colormaps of the photonic weight Ph(n) and of the participation ratio PR(n) as a function of both the eigenfrequencies ω pol n and disorder strength W in units of the average dipole frequency ω0. (c)-(d) Horizontal cuts of panels (a)-(b) at the eigenfrequencies ω pol n /ω0 = 0.9, 1.0, and 1.06, pointed out in panels (a)-(b) by orange, red, and blue triangles, respectively. The cavity height is Lx/a = 12, the number of dipole is N = 1000, and the data have been averaged over 100 disorder realizations.

  3(a) and 3(b)]. In that respect and as detailed in the next subsection, for a given large enough disorder strength W 0.3 the model hosts three distinct phases as can be seen from Figs. 3(a) and 3(b): an exponentially localized, dark state phase composed of the eigenfrequencies ω pol n ω ph q=0 (black region), a phase of mostly photonic states, not yet affected by disorder due to their eigenfrequencies ω pol n ω 0 + W/2 (yellow-orange region), and a phase of mainly dipolar polaritons with an inter-mediate PR, made of states in the dark state band with eigenfrequencies ω ph q=0 ω pol n ω 0 + W/2 (purple region). The latter intermediate phase presents the most interesting properties. We highlight them in Figs. 3(c) and 3(d), by showing horizontal cuts of panels (a) and (b), for the eigenmodes corresponding to the eigenfrequencies ω pol

78 PRFIG. 4 .

 784 FIG.[START_REF] Gonzalez-Ballestero | Uncoupled dark states can inherit polaritonic properties[END_REF]. Probability density |Ψi(n)| 2 along the sites i of a chain of N = 1000 dipoles, for different eigenfrequencies ω pol n /ω0 = 0.9, 1.0, and 1.06 (rows), and different disorder strength W/ω0 = 0.01, 0.1, and 1 (columns). In the figure, data are not averaged over disorder realizations and the cavity height Lx/a = 12.

FIG. 5 .

 5 FIG. 5. (a),(c) Photonic weight Ph(n) and (b),(d) participation ratio PR(n) as a function of both the reduced eigenfrequencies ω pol n /ω0 and cavity height Lx/a, for disorder strengths W/ω0 = 0.01 [left panels (a)-(b)] and W/ω0 = 1 [right panels (c)-(d)]. The data have been averaged over 100 disorder realizations.

  LP branch, that we previously coined the dark state band, remains composed of almost fully dipolar states [see panel(a)]. Looking at the corresponding PR in panel (b), one sees that the modes near the band edges of this dark state band are more localized than the ones in the middle, confirming the behavior of an effective band.In the right panels, Figs. 5(c)-5(d), we consider now the case of a larger disorder strength W/ω 0 = 1. As we increase the cavity height L x /a, we observe an increase of the share of semilocalized states, visible in purple both in Figs.3(a)-3(b) and 5(c)-5(d). The frequency range occupied by these semilocalized states matches the part of the dark state band overlapping with the bare photonic band, thus being bounded from below by the lowest bare photonic state ω ph

L x /a = 7 L x /a = 12 W/ω 0 = 0 FIG. 6 .

 71206 FIG.6. Steady-state amplitude of the dipole moment |pi| on site i in units of the reduced Rabi frequency ΩR/ω0 in an ordered (W/ω0 = 0) chain of N = 2500 dipoles, resulting from a monochromatic drive on the first dipole at a frequency ω d /ω0 = 1.0, corresponding to a dark state. Results are shown for different cavity heights Lx/a = 7 (red lines, weak coupling) and Lx/a = 12 (blue lines, strong coupling), and different dipole losses γ/ω0 = 0.001 (solid lines) and γ/ω0 = 0.02 (dotted lines). The dashed (dotted) grey line represents the analytical estimate[START_REF] Lerario | High-speed flow of interacting organic polaritons[END_REF] with γ/ω0 = 0.001 (γ/ω0 = 0.02).

W/ω 0 = 0 FIG. 7 .

 07 FIG.7. Same quantity as in Fig.6, with dipole losses (a) γ/ω0 = 0.001 and (b) γ/ω0 = 0.02, but for a fixed cavity height Lx/a = 12 and different driving frequencies ω d /ω0 = 0.9, 1.0, and 1.06, exciting respectively a mainly photonic polariton, a dark state, and a mainly dipolar polariton. The grey dashed lines show the analytical formula for ω d /ω0 = 1.0 given in Eq.[START_REF] Lerario | High-speed flow of interacting organic polaritons[END_REF].

FIG. 8 .

 8 FIG.8. Same quantity as in Fig.6, with driving frequency ω d /ω0 = 1.0 corresponding to a dark state and dipole losses (a) γ/ω0 = 0.001 and (b) γ/ω0 = 0.02, but for a fixed cavity height Lx/a = 12 and increasing disorder strength W . The data have been averaged over 100 disorder realizations.

3 W

 3 FIG.9. Steady-state amplitude of the last dipole moment of the chain, |pN |, normalized by the one obtained without disorder, and as a function of the disorder strength W in units of ω0. A chain of N = 500 dipoles is considered, the driving frequency ω d /ω0 = 1.0, the cavity height Lx/a = 12, and the data have been averaged over 100 disorder realizations.

FIG. 10 .

 10 FIG. 10. Steady-state amplitude of the dipole moment |pi|, normalized by the one obtained without disorder, as a function of both the site i along the chain and the disorder strength W in units of ω0. The dipole losses γ/ω0 = 0.02, the driving frequency ω d /ω0 = 1.0, the cavity height Lx/a = 12, and the data have been averaged over 100 disorder realizations.

FIG. 11 .

 11 FIG. 11. Density plots of the steady-state amplitude of the dipole moment |pi| in units of the reduced Rabi frequency ΩR/ω0, as a function of both the site i and the driving frequency ω d /ω0. Three are presented: (a) Without disorder, (b) W/ω0 = 0.1, and (c) W/ω0 = 1. The dipole losses γ/ω0 = 0.02, the cavity height Lx/a = 12, and the data have been averaged over 100 disorder realizations.

  FIG.12. Polaritonic band structure of an ordered chain in the thermodynamic limit, in units of the bare frequency ω0, and as a function of the reduced wavenumber qd0 in the first Brillouin zone. In orange, we show the dispersion used in the main text, which is obtained using the RWA and given by Eq. (A6), and we compare it to the one obtained without using the approximation, given by Eq. (B3) and shown in green. The considered cavity height is Lx/a = 12.

1 FIG. 14 .

 114 FIG.[START_REF] Lerario | Room temperature Bloch surface wave polaritons[END_REF]. Participation ratio PR(n) [defined Eq. (14)] as a function of the size of the chain N , for the eigenmodes associated to the eigenfrequencies ω pol n /ω0 = 0.9, 1.0, and 1.06, with increasing disorder strength from the left to the right panel. Panel (c) shows in addition the result for the eigenmodes associated to the eigenfrequencies ω pol n /ω0 = 0.7 and 1.55. Extended states show a participation ratio growing linearly with N , while localized states present a constant participation ratio. The black dashed line represents the maximum growth rate 2(N + 1)/3. The cavity height Lx/a = 12, and the data have been averaged over 100 disorder realizations.

FIG. 15 .

 15 FIG.[START_REF] Orgiu | Conductivity in organic semiconductors hybridized with the vacuum field[END_REF]. Average of the level spacing ratio rn for eigenstates within the dark state band as a function of the disorder strength, for increasing sizes of the chain N . The grey dashed line shows the exact value for Poisson statistics rn PE = 2 ln 2 -1 0.39. The cavity height Lx/a = 12, and the data have been averaged over 100 disorder realizations.

12 FIG. 16 .

 1216 FIG.[START_REF] Zhong | Energy transfer between spatially separated entangled molecules[END_REF]. Multifractal exponent of the semilocalized states τq(n) [as defined in Eq. (C3)] as a function of q. The exponents have been extracted according to the procedure explained in the text (see colored dots). The black dashed line represents the analytical result found for a simplified disordered Tavis-Cummings model[START_REF] Dubail | Large random arrowhead matrices: Multifractality, semilocalization, and protected transport in disordered quantum spins coupled to a cavity[END_REF]. The inset exemplifies the scaling procedure for the value q = 0.12, with the circles showing numerical data. The cavity height Lx/a = 12 and the disorder strength W/ω0 = 1.

  Appendix D: Steady-state solution of the Lindblad master equation

9 W/ω 0 = 0 FIG. 17 .

 9017 FIG. 17. Steady-state amplitude of the dipole moment |pi| as a function of the dipole sites i for values of the cavity loss parameter κ, in the case of an ordered chain. Panels (a) and (b) correspond to the driving, respectively, of a dark state with ω d /ω0 = 1.0 and of a mainly photonic polariton with ω d /ω0 = 0.9. The grey dashed line in panel (a) represents the exponential decay of Eq.[START_REF] Lerario | High-speed flow of interacting organic polaritons[END_REF]. The dipole losses γ/ω0 = 0.001, and the cavity height Lx/a = 12.

L x /a = 7 L x /a = 12 W/ω 0 = 0 FIG. 18 .

 712018 FIG. 18. Absorption spectrum as the sum of N = 1000 dipole moment amplitudes |pi| in units of the reduced Rabi frequency ΩR/ω0, as a function of the reduced driving frequency ω d /ω0. Results are shown without disorder for two different cavity heights Lx/a = 7 (red lines, weak coupling) and Lx/a = 12 (blue lines, strong coupling), and two different dipole losses γ/ω0 = 0.001 (solid lines) and γ/ω0 = 0.02 (dotted lines).

FIG. 19 .

 19 FIG.[START_REF] Hou | Ultralong-range energy transport in a disordered organic semiconductor at room temperature via coherent exciton-polariton propagation[END_REF]. Absorption spectrum for different disorder strengths W , normalized by the maximum of the one obtained without disorder. The results are shown for a chain of N = 500 dipoles, with a cavity height Lx/a = 12, dipole losses γ/ω0 = 0.001, and have been averaged over 100 disorder realizations.

N

  i=1 |p i (W = 0)|, has been normalized by the maximum of the one found without disorder, max[ N i=1 |p i (W = 0)|].

As an example, in the case of spherical metallic nanoparticles hosting localized surface plasmons, Q = -eNe and M = meNe, where Ne is the number of valence electrons of charge -e < 0 and mass me. The resonance frequency of the plasmon is then the Mie frequency of the nanoparticle ω 0 =

 (Nee

[START_REF] Feist | Extraordinary exciton conductance induced by strong coupling[END_REF] /mea 3 ) 1/2 , with the dipole length scale a being the nanoparticle radius. Note that throughout this paper, we use cgs units.

In contrast to single cavity mode models such as the Tavis-Cummings one, in a realistic multimode model there is no unambiguous definition of dark and polaritonic states[START_REF] Ribeiro | Multimode polariton effects on molecular energy transport and spectral fluctuations[END_REF], all the eigenstates being hybridized.

This unusual transition, occurring here at a cavity height Lx/a 9.1,

Nevertheless, due to cavity confinement, the small heights Lx/a 9 that lead here to the weak-coupling regime are not comparable to the weak

ACKNOWLEDGMENTS

We would like to acknowledge Stéphane Berciaud, Charles A. Downing, David Hagenmüller, Rodolfo A. Jalabert, Gaëtan Percebois, and Dietmar Weinmann for insightful discussions. This work of the Interdisciplinary Thematic Institute QMat, as part of the ITI 2021-2028 program of the University of Strasbourg, CNRS, and Inserm, was supported by IdEx Unistra (ANR 10 IDEX 0002), and by SFRI STRAT'US Projects No. ANR-20-SFRI-0012 and No. ANR-17-EURE-0024 under the framework of the French Investments for the Future Program.

Appendix A: Analytical diagonalization in the case of an infinite ordered chain

In this Appendix, we reproduce for the sake of selfcontainedness the diagonalization procedure of Ref. [START_REF] Downing | Polaritonic Tamm states induced by cavity photons[END_REF] in the case of an infinite ordered dipole chain within a cuboidal photonic cavity, and give analytical expressions for the eigenfrequencies and eigenstates of the system. As in the main text, we here employ the rotating wave approximation. In such translationally invariant case, one can use periodic boundary conditions, and it is then convenient to use the Fourier transform b j = N -1/2 q e ijqd0 b q (j = 1, . . . , N ), where the dipolar wavenumber q = 2πp/N d 0 , with the integer p ∈ [-N /2, +N /2]. Importantly, the periodicity of the system in the z direction implies that the photonic wavenumber k z = 2πm z /L z , m z ∈ Z, where L z = (N + 1)d 0 , is conserved with the dipolar one. In the following, we then write k z = q.

By quantizing the vector potential (8) in a cavity with periodic boundary conditions in the z direction [START_REF] Kakazu | Quantization of electromagnetic fields in cavities and spontaneous emission[END_REF], and by rewriting the dipolar Hamiltonian (2) in Fourier space, the polaritonic Hamiltonian [START_REF] Engelhardt | Unusual dynamical properties of disordered polaritons in microcavities[END_REF] becomes

where H q is the Bloch Hamiltonian, while φ † q = (b † q , c † q ). In the above expression, ω dp q is the quasistatic dispersion of the collective dipolar excitation, i.e., the eigenspectrum of the dipolar Hamiltonian (2). It reads

with the dipolar coupling strength Ω 0 = (ω 0 /2)(a/d) 3 , and where the lattice sum f q can be expressed in closed form in terms of the polylogarithm function Li s (z) = ∞ n=1 z n /n s as f σ q = Li 3 (e iqd ) + Li 3 (e -iqd ). Within such periodic boundary conditions, the photon dispersion [START_REF] Patton | Allphotonic quantum simulators with spectrally disordered emitters[END_REF] is

and the light-matter coupling (12) now reads

Using a bosonic Bogoliubov transformation, one can readily diagonalize the Hamiltonian (A1) as

where the upper (τ = +) and lower (τ = -) polaritonic bands are given by

with ∆ q = (ω ph qω dp q )/2 being the light-matter detuning between the bare photonic and dipolar dispersions.

Finally, the Bogoliubov operators diagonalizing the Hamiltonian (A1) are a linear combination of the dipolar and photonic ladder operators, β qτ = u qτ b q + v qτ c q . The modulus squared of the two coefficients u qτ and v qτ , which are normalized as |u qτ | 2 + |v qτ | 2 = 1, thus represent, respectively, the dipolar part D qτ and the photonic part Ph qτ of the polaritonic eigenmodes. These latter quantities read

and

Appendix B: Beyond the rotating wave approximation

Here, we provide justifications for the rotating wave approximation (RWA) used in the main text. To this end, we consider the effects of counter-rotating terms in both the eigenvalues and transport properties of an ordered chain.

Considering the counter-rotating terms, the polaritonic Hamiltonian (10) that we use in the main text transforms into