Anderson localisation in steady states of microcavity polaritons

We present an experimental signature of the Anderson localisation of microcavity polaritons, and provide a systematic study of the dependence on disorder strength. We reveal a controllable degree of localisation, as characterised by the inverse-participation ratio, by tuning the positional disorder of arrays of interacting mesas. This constitutes the realisation of disorder-induced localisation in a driven-dissipative system. In addition to being an ideal candidate for investigating localisation in this regime, microcavity polaritons hold promise for low-power, ultra-small devices and their localisation could be used as a resource in quantum memory and quantum information processing.

The localisation of diffusive waves due to an underlying disorder is pervasive throughout nature. A Nobel-prize winning description of the localisation of single-particle electronic wavefunctions in periodic lattices was provided by Anderson in 1958 1 . The essence being that for sufficiently strong on-site disorder, the eigenstates transition from extended Bloch states to exponentially localised. Since then, localisation has been realised in several distinct physical systems such as photonic waveguides 2 , microwaves 3,4 and sound waves 5,6 . True Anderson localisation of matter waves remained elusive for a long time, and was only realised in an atomic Bose-Einstein condensate in 2008 7 . In addition, the question of whether single particle states remain localised in the presence of inter-particle interactions was eventually resolved positively 8,9 , and has led to the field of many-body localisation 10 .
In general, localisation prevents the system from reaching the equilibrium state of the corresponding non-disordered system. The most striking example can be found in closed systems of cold atoms, where a random disorder potential prevents the system from reaching thermal equilibrium 7,11,12 . In contrast to closed regimes, we are interested in localisation in the steady states of driven-dissipative systems, the solutions that arise from the balance of flows from driving and dissipation. It is not immediately obvious if the disorder-induced Anderson localisation has a role to play in this regime, and so it is an interesting question to explore experimentally.
Microcavity polaritons are an ideal test-bed for this. They are hybrid bosonic quasiparticles that arise due to the strong interaction between microcavity photons and quantum well excitons 13 . They have lifetimes on the order of picoseconds, and thus the polariton population must be continually replenished by a laser source in order to form a steady state. Due to advanced fabrication techniques, the potential landscape experienced by the polaritons can be exactly engineered, and their distribution can be directly mapped by measuring the photoluminescence. We note that microcavity polariton localisation has been reported in a few different contexts such as: nonlinear-induced localisation when resonantly driving one micropillar of an interacting dimer 14,15 ; gain-induced localisation in finite width excitation spots 16 , a consequence of a finite polariton lifetime; and localisation in flat bands where the polaritons have an infinite effective mass 17 . However, to date, no systematic study of the effect of disorder-induced Anderson localisation has been performed.
Here we experimentally demonstrate a signature of the Anderson localisation of exciton-polaritons in two-dimensions. To do this, we study the steady-state polariton distribution under non-resonant excitation in a set of eight hexagonal lattices with increasing levels of static disorder. The static off-diagonal disorder is introduced by adding a random displacement to each lattice site with a controllable maximum amplitude. The localisation is characterised by the inverse participation-ratio, and is shown to monotonically increase as a function of disorder strength. The experimental results are supported by numerical simulations of a Gross-Pitaevskii equation, with which we also explore the effect of the polariton nonlinearity.

Results
Trapping potentials for the polaritons are designed by patterning the microcavity with mesas 18,19 , circular elongations of the cavity spacer that locally alter the cavity detuning. See Fig. 1 for a schematic of the device. We note that mesas are distinct from micropillars within which the photons are fully confined. We arrange the mesas in a hexagonal lattice, and close enough that the quantised modes of the individual mesas 18 hybridise into extended modes 19 whose features depend on the lattice geometry. The polariton population is fed by a reservoir of excitons created by the non-resonant excitation of a continuous-wave laser. From the photoluminescence images of the perfect hexagonal arrays, Fig. 2(a), we see the confinement of polaritons predominantly within the mesas and an approximately homogeneous distribution among them. We note that the coupling between the mesas was  In the absence of disorder we see a more homogeneous distribution among the mesas, whereas disorder induces the onset of patches of localisation. Also shown are plots that reveal how the IPR increases with disorder for (i) experiment, and (j) theory. The polariton nonlinearity in (j) is g = g 0 = 2.4 × 10 −3 meV.μm 2 . The inset of (j) shows the result for g = −g 0 (blue dotted line), g = 0 (green dot-dashed line), g = g 0 (red solid line), and g = 10g 0 (black dashed line). The lines are just a guide for the eye. In the experimental figures the cavity-exciton detuning is 2 meV. The simulation parameters are m = 5 × 10 −5 m e , R = 0.4 meV.μm 2 , γ = 0.5 meV, g = 2.4 × 10 −3 meV.μm 2 , V 0 = 9 meV, γ R = 2 meV, and P 0 = 2γ R γ/R. Here, V 0 is the maxima of the trapping potential, which we model as a radially symmetric sigmoid function for each mesa.
evidenced in the spectral distribution of the polaritons (see Supplementary Fig. S3) where one can see the periodic Bloch bands of the hexagonal lattice.
To introduce an off-diagonal disorder we shift the cartesian coordinates of each mesa by a random displacement in the range d[−δ, δ], where 0 ≤ δ ≤ 1 parameterises the amount of disorder, and d = 0.25 μm is the maximum possible displacement for the maximum disorder δ = 1. This positional disorder modifies the eigenstates of the system from Bloch states towards spatially separated patches of localisation. The effect of the disorder on the localisation, or clustering, of the polariton population can be seen by eye in the photoluminescence images, Fig. 2(a-d). To obtain a quantitative measure of the amount of localisation we calculate the inverse-participation ratio (IPR), which is functionally similar to imbalance measures used in cold atom experiments 11,12 , although it is not a site specific measure. The IPR has previously been used to quantify Anderson localisation in photonic systems 20 and is also applicable here. In essence it is a measure of inhomogeneity, and for a homogeneous distribution it is equal to unity. Thus as the onset of Anderson localisation causes some mesas to contribute more significantly to the total photoluminescence, the IPR increases also. To exploit this measure, we first normalise the data to account for the Gaussian background (see Methods), and then calculate the average occupation  Figure 2(i) shows the percentage change in the IPR from no disorder (δ = 0). We clearly see the increase in IPR with increasing disorder δ, which signals the onset of localisation. The error bars correspond to the standard error after repeating the experiment on 12 different regions sampled from a larger lattice for each disorder strength. In addition, similar results have been reproduced for several different laser powers (see Supplementary Figure S1).
We successfully modelled the experimental results with a generalised Gross-Pitaevskii equation 21 describing the evolution of the polariton wavefunction ψ(r, t), where m and γ are the effective mass and decay rate of the polaritons, g is the strength of polariton-polariton interactions, R is the reservoir-polariton exchange rate, and V(r) is the potential landscape defined by the lattice. The reservoir n R (r, t) is described by the rate equation where γ R is the decay rate of the reservoir. The reservoir is populated by the continuous-wave pump P(r) which we model as a Gaussian with amplitude P 0 . We use the Runga-Kutta method of fourth order to evolve the dynamics until a stationary solution is achieved (approximately 50 ps). In Fig. 2(e-h) we show the results of these simulations for a periodic and increasingly disordered arrays. We can see the onset of patches of localisation when a disorder is introduced. We calculate the IPR using Eq. 1 in much the same way as is done for the experimental data, and plot the results in Fig. 2(j).
In the inset of Fig. 2(j) we show calculations performed for different polariton-polariton nonlinearities. The linear regime is shown by the green dot-dashed line. In the present experiment we are working in a weak nonlinearity regime, which we model with a small g in the simulations (red solid line). Nonetheless we observe that the positive interaction acts to suppress the localisation; see also the black-dashed line where we increase the non-linearity ten-fold. In addition we show that a negative interaction acts to enhance the localisation (blue www.nature.com/scientificreports www.nature.com/scientificreports/ dotted line). Although a negative g is not possible with our experimental setup, we include this simulations result as a point of interest for the reader. Such a regime could be accessed with spinor condensates tuned near the Feshbach resonance 22 .
In Fig. 3 we provide some further theoretical analysis, demonstrating the ability to tune the localisation through additional system parameters. For example, simply by varying the detuning of the cavity via the spacer width we can vary e.g. the polariton lifetime. As we see in Fig. 3(a), decreasing the lifetime of the polaritons results in an enhanced localisation. As another example we can modify the reservoir-polariton coupling, which could be accomplished e.g. by enhancing phonon mediated relaxation of excitons into polariton states. In this way we can suppress the localisation by effectively increasing the driving of the polaritons, as shown in Fig. 3(b).

Discussion
Our results demonstrate a signature of Anderson localisation in the steady states of microcavity polaritons. In general we observe a convincing qualitative agreement between the experimental and numerical results. In both cases the IPR increases monotonically with disorder, and by a similar magnitude, which signals the onset of localisation. The difference between the rates of increase of IPR is attributed to intrinsic disorder of the experimental samples, as well as slightly different data processing algorithms. In this work we were operating in a regime of small nonlinearity where interactions acted to weakly suppress the localisation.
Our results offer a signature of disorder-induced localisation in the steady states of driven-dissipative systems, a regime quite different to the prototypical case of closed systems. We believe that the phenomenology reported here should be generally observable in other driven-dissipative systems and may open a new chapter for basic science explorations. In addition, since both polaritons and localisation are foreseen to have a high potential for applications in optoelectronic devices and quantum information respectively, such a robust and controllable phenomenon could be of use in novel devices.
To advance our work, it would be interesting to investigate localisation of strongly interacting polaritons. For example, one could explore the interplay between the disorder-induced localisation seen herein and effects such as nonlinear localisation observed for strongly driven microcavities 14 . Moreover with larger polariton densities it may be possible to extract signatures of many-body localisation, which in a crude sense is the persistence of Anderson localisation in the presence of many-body interactions. One could also explicitly examine the role of localisation in driven-dissipative systems for preserving memory of initial conditions. For example, by preparing initially imbalanced population distributions, with a highly inhomogeneous pumping, and then switching on a homogeneous pumping to see if a signature of the initial state perseveres.

Methods experimental methods.
We study a GaAs λ-microcavity made of 24/20 pairs of GaAs/AlAs, with embedded 8 nm, In 0.04 Ga 0.96 As quantum well (with an exciton energy of 1.482 eV) which gives a Rabi splitting of 3.1 meV. We fabricate circular mesas with a radius of r = 1 μm by a 6 nm local elongation of the cavity spacer, which provides a trapping potential of 9 meV for the polaritons, leading to confined quantised modes in the individual mesas 18 . We arrange the mesas into a hexagonal array with lattice constant a = 2.5 μm which is sufficient for wavefunction overlap between neighbouring mesas, giving rise to new hybridised modes 19 . To introduce an off-diagonal disorder to the system the x and y coordinates of the mesas are offset by a random value in the range d [−δ,δ] where d = 0.25 μm. We consider eight different disorder levels from δ = 0 to δ = 1 in evenly spaced steps. We excite the system non-resonantly with a 660 nm continuous wave laser focused to a spot size of 25 μm. We measure the photoluminescence of the sample with a collection lens of numerical aperture 0.42NA, and image the real-space integrated energy emission in a CCD.
We are interested in the onset of localisation phenomena as evidenced by the increase in the IPR. Thus we renormalise the experimentally obtained signal by the pump beam profile to eliminate its effect on the IPR. Thus we fit a Gaussian to the image after filtering out higher frequencies, which leaves us with the overall shape of the pump profile, after relaxing from the higher nonresonant energy. We then crop the image to a region within the central pumping region where the signal-to-noise ratio is sufficient. Then we use a peak-finding algorithm to determine the location of the mesas. We then calculate the average intensity I n of each mesa (labelled by the integer n) in a region around these peaks of the same width as the mesas. The IPR can then be calculated from Eq. (1). numerical methods. We use the Runge-Kutta method of fourth order to solve Eq. (2). We denote the mesa positions as = + + δ n n d R a a n 1 1 2 2  , where a 1 = a(0, 1) and = a a ( /2)(1, 3) 2 are the lattice vectors, a is the lattice constant, and n labels the integers n 1 and n 2 . Here, δ  is a random vector whose coordinates are uniform random variables sampled in the range [−δ, δ], 0 ≤ δ ≤ 1 parametrises the disorder, and d = 0.25 μm is the maximum possible displacement in each direction. We take the zero-energy to be E 0 , the bottom of the polariton band. So the wavefunction ψ is technically the slowly oscillating envelope of the real wavefunction Ψ = ψexp(−iE 0 t/). This is done to minimise numerical errors that can accumulate with fast oscillations.
We seed the process with a weak initial state ψ(t = 0) = P 0 /10 and then allow the system to evolve until the stationary solution is achieved (approximately 50 ps). To obtain the IPR of the solution we process the data in an analogous way to that done experimentally. First we normalise the wavefunction by that solution ψ back which is obtained by simulating the system without any mesas, in other words the effective 'background' of the polaritons. Then we calculate the average density ρ n of polaritons in the vicinity of each mesa (labelled by n). We then discard those mesas located outside the pumping area, i.e. those mesas with positions |R n | > 2σ. Finally, we calculate the IPR as per Eq. (1) in the main text. Schematically the process looks like 1. Obtain density normalised to background: ρ(r) = |ψ(r)| 2 /|ψ back (r)| 2 2. Calculate average density at each mesa: ρ n = 〈ρ(r)〉 |r − Rn| < a 3. Discard mesas outside the pumping area: S = {n ∈ N, |R n | ≤ 2σ} 4. Calculate IPR for the array of mesas:  ρ ρ = ∑ ∑ ∈ ∈ N( )( ) n S n n S n 2 2 Here, N is the number of mesas used in the calculation (length of set S). The simulation parameters are m = 5 × 10 −5 m e , R = 0.4 meV.μm 2 , γ = 0.5 meV, g = 2.4 × 10 −3 meV.μm 2 , V 0 = 9 meV, γ R = 2 meV, and P 0 = 2γ R γ/R. Here, V 0 is the maxima of the trapping potential, which we model as a radially symmetric sigmoid function for each mesa.

Data availability
All raw data and source code is available from the corresponding author upon reasonable request.