Role of Multiple Scattering in Single Particle Perturbations in Absorbing Random Media

Speckle patterns produced by disordered scattering systems exhibit a sensitivity to addition of individual particles which can be used for sensing applications. Using a coupled dipole model we investigate how multiple scattering can enhance field perturbations arising in such random scattering based sensors. Three distinct families of multiple scattering paths are shown to contribute and the corresponding complex enhancement factors derived. Probability distributions of individual enhancement factors over the complex plane are characterised numerically within the context of surface plasmon polariton scattering in which absorption is shown to play an important role. We show that enhancements become more strongly dependent on individual scatterer properties when absorption losses are larger, however, amplitude enhancements $\sim 10^2$, comparable to low loss surface plasmons, are achievable through sensor optimisation. Approximate analytic expressions for the complex mean enhancements are also found, which agree well with simulations when loop contributions are negligible.


I. INTRODUCTION
Use of optical scattering for detection and measurement is a powerful and widespread approach underpinning techniques such as interferometric scattering microscopy (iSCAT) [1], dynamic light scattering [2] and diffusing wave spectroscopy [3]. This phenomenon has seen extensive application in the biological sciences and environmental monitoring, in turn driving development of scattering based sensors. As sensitivity gains have been made, so sensing has moved from monitoring of bulk properties to detection of individual nanometer-sized analyte particles, such as virions and proteins [4]. Such small dimensions however mean particles only scatter weakly, hence presenting a major challenge. To mitigate this issue strongly confined fields, which can enhance light-matter interactions, have been employed, for example, using high Q optical resonators [5], photonic crystals [6] and nanoapertures [7]. Plasmonic systems, supporting localized or propagating surface plasmon-polaritons (SPPs), are also particularly attractive for sensing, since in addition to confining optical fields they can be easily implemented on chip scale devices, are biocompatible, allow operation in aqueous/microfluidic environments and can exploit the existing wealth of functionalisation protocols required to maintain specificity [8,9]. Accordingly, SPPs have found applications in numerous sensing and particle tracking setups [10][11][12][13][14][15].
Interferometric plasmonic systems provide one route to yet further sensitivity gains [16][17][18][19][20] by leveraging coherent detection. Alternatively, nanostructured sensors, such as nanoparticle arrays, aperiodic gratings and randomly disordered substrates [21][22][23][24][25] have also shown significant promise. Random scattering, in particular, affords numerous opportunities in sensing by virtue of the * matthew.foreman@imperial.ac.uk diverse range of phenomena that can occur. For example, depending on the degree of multiple scattering (as parametrized by the scattering mean free path), scattering can give rise to long and short range correlations, weak and strong (Anderson) localisation and fluctuations in the local density of states. A substantial amount of work has been dedicated to study such phenomena in both the optical [26][27][28][29][30][31] and plasmonic domains [32][33][34][35][36][37][38] over the years. Indeed, exploitation of random scattering has a rich track record in optics. For example, correlations present in the speckle patterns have been used for refractive index sensing, spectrometry and imaging [39][40][41][42]. Speckle patterns generated by disordered multiple scattering environments have also been shown to depend on the properties of individual scatterers [43,44], such as their position, effective charge or orientation [45][46][47], whilst also providing enhanced sensitivity as compared to single scattering environments [48]. Approaches to extract the position of a single scatterer accounting for multiple scattering effects have thus been developed, for example based on diffusive models of light propagation [49] or extension of single scattering holography localisation techniques [50,51]. Recent advances in machine learning moreover present further opportunities to extract information from randomly scattered light, since such approaches do not require a detailed physical model and are hence applicable across a broad range of scattering regimes [52][53][54].
In plasmonics, random scattering has also seen employ, for example, in light harvesting, super-resolution imaging and sensing [25,[55][56][57][58]. Absorption associated with Ohmic losses in metals is, however, intrinsic to plasmonic systems [59]. For resonance tracking based sensors, absorption broadens the resonance lineshape and thus limits sensitivity. Statistical properties of speckle patterns in absorbing multiple scattering environments can however also be affected, for example, absorption can give rise to non-Rayleigh intensity statistics, as well as generate non-universal and reduced correlations [37,[60][61][62].
In this article, we address the open question as to how absorption affects sensitivity in random scattering based particle sensing. Particular emphasis is placed on surface plasmon based setups due to their prevalence and capabilities in this domain. To address this question, in Sections II A and II B we first derive three enhancement factors, arising from three distinct categories of multiple scattering paths, which describe the effect of multiple scattering on the electric field perturbation caused by the presence of an additional analyte particle. We recently studied the amplitude of these enhancement factors in the context of multiple scattering of SPPs [48] by randomly distributed scatterers on a metal surface, however here we study the full probability distribution of the enhancements on the complex plane, including phase effects. Approximate analytic results for the mean enhancement factors are derived in Section II C, before numerical results are given in Section III. The dependence of the achievable enhancements and the associated probability distributions on absorption loss is detailed in Section III B. In particular, through consideration of the role of scattering phase, propagation phase and absorption we identify a non-trivial dependence of the mean enhancement on tunable properties of the scattering configuration. This dependence is explored as a route to sensor optimisation in Section III C. As such the insights gained in this work allow us not only to understand the interplay of absorption and multiple scattering upon addition of an analyte particle, but also to guide future development of optimal random SPP sensors.

A. Coupled Dipole Model
The disordered scattering system we shall study is a collection of N coupled point dipole scatterers [63][64][65], situated in an environment with background dielectric function ε(r). A Green's tensor G(r, r ) can be defined for this system as the solution to Maxwell's wave equation where k 0 = ω/c = 2π/λ 0 , ω is the angular frequency, c is the speed of light, λ 0 is the wavelength in vacuum and I is the 3 × 3 identity matrix. When the point scatterers are illuminated with a monochromatic incident electric field E 0 (r), the total electric field E(r) at position r is where r j , α j and p j = α j E exc (r j ) are the position, dressed polarizability and dipole moment of the jth scatterer respectively, and E exc (r j ) = E 0 (r j ) + i =j G(r j , r i )p i is the exciting field incident on the jth dipole, consisting of the incident field and the field from all other dipoles [63,66]. Notably, α j includes the effect of self-interactions (e.g. due to reflections from the background medium). From Eq. (2) we can construct the set of linear equations for i = 1, 2, . . . N , where p 0,i = α i E 0 (r i ) is the dipole moment induced by the incident field in the ith scatterer, the matrix elements M ij are defined by for i, j = 1, 2, . . . N , and G ij = G(r i , r j ). Once Eq. (3) is solved for the N dipole moments, the field at any point can be calculated using Eq. (2). Throughout this analysis, we consider scattering of a vector field with corresponding Green's tensor, such that M ij are the tensor elements of an N × N matrix of tensors (or equivalently they are the 3 × 3 blocks making up a 3N × 3N matrix), which we denote M . Our analysis, however, is equally valid for scattering of a scalar field, if G, α, E and p i are replaced with scalar equivalents, in which case M ij are the scalar elements of an N × N matrix. For random positions, r i , the matrix M is a Euclidean random matrix, the statistics of which have been studied, for example, in the context of optical scattering and vibrational modes of glasses [67][68][69][70]. Within the single scattering regime, the off-diagonal terms describing coupling between the dipoles are negligible such that M ij ≈ Iδ ij and p i ≈ p 0,i .

B. Adding a Scatterer
We now consider perturbing the scattering configuration by introducing an additional point scatterer with polarizability α N +1 at position r N +1 . The perturbed system can be described similarly to above yielding the set of coupled dipole equations N +1 j=1 M ij p j = p 0,i (i = 1, 2, . . . N + 1) in terms of the modified dipole moments p j . We note that the matrix elements M ij for the perturbed system are again given by Eq. (4) albeit the indices i and j run from 1 to N +1 (hence M ij = M ij , for i, j ≤ N ). The new set of N + 1 dipole moments results in the perturbed field E (cf. Eq. (2)) Accordingly, the perturbation to the field δE = E − E caused by the addition of the scatterer is hence where δp j = p j − p j is the perturbation to the jth dipole moment and, since there is no (N + 1)th scatterer in the unperturbed system, we have dropped the prime from p N +1 . The first term of Eq. (6) corresponds to the field scattered by the added dipole p N +1 , whereas the second term arises because multiple scattering introduces dipole coupling whereby the presence of the additional scatterer modifies the N initial dipole moments. In the single scattering regime, the coupling between dipoles is negligible so that δp i = 0 and the second term vanishes. Similarly, p N +1 = p 0,N +1 such that the single scattering perturbation δE ss (r) reduces to The coupled dipole equations for the perturbed N + 1 scatterer system can be expressed in the form Using Eq. (3), Eq. (8) can in turn be rearranged to yield where M −1 ij is here used to denote the the (i, j)th 3 × 3 block (corresponding to rows 3i − 2 to 3i and columns 3j − 2 to 3j) of the inverse of the entire 3N × 3N matrix M , as opposed to (M ij ) −1 , the inverse of the 3 × 3 submatrix M ij (similarly, in the scalar case, it corresponds to the (i, j)th element of the inverse of the N × N matrix M ). Substituting Eq. (10) into Eq. (6) then gives where we have also defined the enhancement factor Expressing δE as such allows comparison with the single scattering result in Eq. (7). Specifically, it is evident that the perturbation to the dipole moments of the N initial scatterers from introduction of an additional scatterer is described by the factor γ 1 . Equivalently, dipole coupling through multiple scattering acts to modify the effective dipole moment of the additional scatterer such that p N +1 → γ 1 (r)p N +1 . The tensor nature of γ 1 reflects the fact that the polarization of the field perturbation can be modified by multiple scattering. Similarly, γ 1 is a complex quantity, implying multiple scattering can affect both the phase and amplitude of δE.
In addition to the dipole coupling captured in γ 1 , there remain further multiple scattering effects which cause p N +1 = p 0,N +1 . Specifically, the local field experienced by the additional scatterer is not solely dictated by the incident field E 0 , but also contains a contribution from scattering of the illumination field by the N initial scatterers. To demonstrate this, we substitute Eq. (10) and Eq. (3) into Eq. (9), which results in Rearranging for p N +1 yields Defining two further enhancement factors allows δE to be expressed as where Expressed in this way, it can be seen that the effect of multiple scattering is equivalent to changing the dipole moment from p 0,N +1 to γ 1 γ 2 γ 3 p 0,N +1 . In general, as with γ 1 , the enhancement factors γ 2 and γ 3 are complex matrices, meaning multiple scattering can change the phase, amplitude and polarization of δE. Each enhancement factor can be associated with a class of multiple scattering paths involving the additional scatterer as shown in Fig. 1. Firstly, the effect of rescattering of the field as it propagates to the observation point r after being scattered by the additional scatterer is accounted for by γ 1 . The factor of α j G j,N +1 freely propagates the scattered field from r N +1 to a scattering event at the jth scatterer, while M −1 ij propagates the field from the jth scatterer to the ith scatterer via all possible scattering paths involving the initial N scatterers. Free propagation from the ith scatterer to the observation point r is then described by G(r, r i ). Secondly, FIG. 1. Example multiple scattering paths for each enhancement factor: (γ1, left) rescattering between scattering from the additional particle and propagation to the observation point, (γ2, center) loop trajectories starting and ending on the additional scatterer and (γ3, right) multiple scattering of the illumination field onto the analyte particle.
the γ 2 factor describes the effect of loop scattering paths in which waves, after being scattered by the additional dipole, return (possibly multiple times) to the additional dipole via multiple scattering from the N initial dipoles. As with γ 1 , a factor of M −1 ij α j G j,N +1 propagates the scattered field from the additional scatterer to the ith scatterer via all possible scattering paths not including the additional scatterer. The factor of α N +1 G N +1,i then propagates the field back to the additional scatterer, from which it is scattered again, completing the loop. Summing over the number of loops yields a geometric series in terms of the single loop factor, and hence γ 2 can be expressed as a matrix inverse. This loop contribution is a self-interaction effect analogous to the surface dressing of polarizability. Finally, γ 3 accounts for the effect of scattering of the incident field onto the additional scatterer and therefore γ 3 describes the hotspot effect [71,72]. The incident field at the jth scatterer is multiply scattered to the ith scatterer, described by M −1 ij , and then propagated to the additional scatterer at r N +1 , as is described by the final factor of α N +1 G N +1,i .
Optical sensing often aims to detect particles at a surface where functionalisation of the surface can allow for specificity. For this reason, we shall henceforth consider scattering configurations in which the N initial scatterers are distributed over a planar surface at positions z i > 0. Physically, the scatterers could hence correspond to surface roughness features, bound receptors or nanoparticles, or nanostructures fabricated on a substrate. We also assume that the observation point r = (ρ, z) is taken in the far field as is the case in most sensing setups, which allows the form of γ 1 (r) to be greatly simplified. Specifically, the far field Green's tensor G ∞ is proportional to the 2D Fourier transform of the near field Green's tensor G(q, z, r ) with respect to the transverse position ρ = (x, y), i.e. [73] where k = n b k 0 is the wavenumber in the far field, n b is the refractive index at z, k = (k x , k y ) = k sin θ(cos φ, sin φ) is the 2D transverse component of the wavevector in the direction of observation and (r, θ, φ) are the standard spherical coordinates of r. Under the assumption of translational invariance in the transverse (x, y) plane and that the scatterers all lie in the same bulk medium of dielectric constant ε d , the far field Green's tensors for two different source positions are related. If the observation point is on the same side (z > 0) as the scatterers then whereas for observation points in the lower half space, for example if a thin film substrate is used, then where k z = ±(ε d k 2 0 − k 2 ) 1/2 , the upper (lower) sign is for observation points above (below) the interface and G dir ∞ , G ref ∞ and G tr ∞ are the direct, reflected and transmitted components of the Green's tensor respectively [73]. Under these assumptions, the Green's tensors in Eq. (12) cancel with the inverse Green's tensor factor, resulting in The function R ± (z i , z N +1 ) is derived in the appendix. Importantly, R − (z i , z N +1 ) = I for observation points z < 0 and R ± (z i , z i ) = I regardless of observation point. The remaining two enhancement factors do not depend on the observation point and thus do not differ between the farfield or near-field. We have so far only considered the perturbation to the electric field, however most experimental setups measure the intensity of light, I = |E| 2 . The intensity perturbation δI = |E | 2 − |E| 2 is therefore the typical signal in scattering based optical sensing, and can be related to the field perturbation through where, in addition to the intensity of the perturbation δE, there is a term corresponding to the interference between the field perturbation and the initial field. If, as will typically be the case for a large number of scatterers, |E| |δE|, the interference term dominates and the intensity perturbation can be significantly larger than the dark field case where only the intensity scattered by the analyte particle is present. This principle is central to iSCAT and related techniques [1,11,74], but it is not a multiple scattering effect and can be achieved equally well within a single scattering regime (and indeed typically is in iSCAT experiments), whether the interference is with other scattered fields or an external reference field. As such, this form of signal enhancement is independent of the scattering regime and different to the enhancement mechanisms we are considering. The phase difference between E and δE is random in both the single and multiple scattering regimes, so that the phase statistics of the interference term are essentially identical. The primary difference in the statistics of the interference term between single and multiple scattering lies in the different amplitudes |δE|.

C. Mean Enhancement Factors
For any given scattering configuration the value of each enhancement factor can be determined, however, it is valuable to characterise the distribution and average properties of the enhancement factors over the ensemble of different random configurations. In the following, the transverse positions ρ i of the initial scatterers are assumed to be independently randomly distributed with uniform probability across a 2D planar region of area L 2 on the surface of a substrate, with the same height z i = z s for i ≤ N . Furthermore, the initial scatterers are assumed to be identical and to have the same orientation relative to the surface, whereby α i = α for i ≤ N . Even for isotropic scatterers α may still be an anisotropic tensor due to surface dressing effects, however since the N scatterers are at the same height, the surface dressing effect is identical for each scatterer. Note that the polarizability α N +1 is not restricted and may be different to the background scatterers. We do however limit attention to the case where R ± (z i , z N +1 ) = I since this matches our simulations below and embodies all relevant physics in spite of the reduced mathematical complexity. Under these assumptions, γ 1 , S 2 and γ 3 can be calculated analytically, where γ 2 = (I − S 2 ) −1 and angled brackets denote averaging over realisations of the N background scatterer positions ρ i . It should be noted, owing to the inverse relationship between γ 2 and S 2 , their statistics have a more complicated relationship than the relationship between γ 1,3 and the corresponding sum terms appearing in Eqs. (12) and (17). Using Fourier analysis, γ 1,3 and S 2 can be expressed as where f (q) denotes the 2D Fourier transform of f (ρ) such that f (ρ) = f (q)e iq·ρ d 2 q/(2π) 2 and the function A(q, q ) is defined by In this form, the dependence of γ 1,3 and S 2 on the background scatterer positions is entirely described by A, such that their statistics are determined solely by the statistics of A. Accordingly, the means of Eqs. (23)-(25) can be calculated from A(q, q ) . In order to calculate A , we use the Neumann series ( Physically, Eq. (27) shows how M −1 ij corresponds to a sum over all scattering paths starting at the jth scatterer and ending at the ith scatterer, with P k ij corresponding to the contribution from all paths visiting exactly k scatterers. Each factor of αG lili+1 propagates the field to the next scattering event. The l i+1 = l i exclusion arises because a scattering path does not visit the same scatterer consecutively (as the self interaction is accounted for in α). Inserting this expansion into A(q, q ), the pth order contribution, denoted A (p) (q, q ) such that A(q, q ) = ∞ p=0 A (p) (q, q ), is given by where henceforth the limits and exclusions from the sums will be left implicit. Replacing each Green's tensor with its Fourier decomposition allows the dependence on the scatterer positions to be included within an exponential factor as follows The only random component of Eq. (29) is the exponential factors. Regrouping the exponents so that each r li term is in one exponential factor allows the sum to be rewritten as i,j,l1,...,lp−1 In general, there are terms in the sum in Eq. (30) where l i = l j even when i = j, meaning each exponential factor is not necessarily independent of the others and hence cannot be averaged individually. Following a similar approach to that taken in Ref. [69], we first consider only terms with no shared indices, where each l i is distinct. For these terms, each exponential e i(qc+1−qc)·ρ lc can be averaged independently from the rest. Since there are p + 1 different scatterers in such terms, there exist N (N − 1)(N − 2) . . . (N − p)) ≈ N p+1 terms in the sum with no repeated scatterers. Averaging a general function f (ρ i ) over a scatterer position ρ i corresponds to the integral f (ρ i ) = f (ρ i )d 2 ρ i /L 2 . Therefore, averaging over the p+1 different scatterer positions gives a factor of (L 2 ) −(p+1) , so that the contribution of these distinct scatterer terms scales as n p+1 , where n = N/L 2 is the areal scatterer density. If we now consider the contribution of terms in the sum with 1 repeated scatterer (corresponding to scattering paths involving loops), meaning p distinct scatterers are visited, choosing p scatterers out of N options gives N (N −1) . . . (N −(p−1)) ∼ N p such terms. In this case, averaging over the p scatterer positions give (L 2 ) −p , so that the contribution of these single repeated scatterer terms to the total sum is ∼ n p . It can be seen that the contribution of terms with r repeated scatterers to the total sum in Eq. (30) scales as n p+1−r . While methods to calculate the contribution from these loop paths exist [69], here we only take the leading order terms in n, i.e. the no loop contributions where all the indices i, j, l 1 , . . . , l p−1 are distinct. Within this approximation, in the limit of large system size and scatterer number, L → ∞ and N → ∞, while keeping the scatterer density n constant, the identity N j=1 e iq ·ρj → n(2π) 2 δ(q ) can be applied for each summation index. After averaging, each exponential factor in Eq. (30) can therefore be replaced with a Dirac δ-function. Thus, the pth order contribution to A can be approximated by Summing over p hence gives (32) The means of Eqs. (23)- (25), to leading order in n, therefore follow and are given by For the simple case of an incident (lossless) plane wave E 0 = A 0ξ exp ik in · r and isotropic polarizabilities, γ 3 reduces to a much simpler form, specifically The similar forms of γ 1 and γ 3 reflect the reciprocal symmetry present between scattering of an incoming plane wave scattering into an outgoing plane wave [75]. A notable feature of Eq. (33) is the divergence when I − n(k 2 0 /ε 0 )α G(k ; z s , z s ) is singular, or in the scalar case, when n(k 2 0 /ε 0 )α G(k ; z s , z s ) = 1. When this condition is close to being satisfied (i.e. det I − n(k 2 0 /ε 0 )α G(k ; z s , z s ) is close to zero), the mean will become very large, suggesting the multiple scattering environment is significantly more sensitive to the addition of a scatterer than the single scattering environment. In the scalar case, | γ 1 | has a maximum value | γ 1 | max > 1 provided Re[α G(−k ; z s , z s ] > 0, occurring at a density n opt,1 , where Analogous expressions for n opt,3 and | γ 3 | max arise in the lossless case, replacing k with k in in the argument of the Green's function. Physically, we can understand these conditions by considering the phase shifts involved in scattering. The plane wave component of the field scattered from one scatterer at wavevector q is phase shifted by arg[α G(q)] relative to the incident field. For any multiple scattering path, this phase shift is acquired at each scattering event, in addition to a propagation phase from travelling between scatterers. On averaging over realisations, the propagation phases cancel out, while the phase shift imparted by scattering events remains constant. When Im[α G(q)] = 0 and Re[α G(q)] > 0, there is no phase shift upon scattering and the averaged multiple scattering paths add up in phase, giving a maximum amplitude which, since the N → ∞ limit has been taken, diverges as there are an infinite number of scattering paths in this case. In turn, a divergence of Eq. (38) results. Of course, any given realisation need not be close to the mean, and the random propagation phase can play a large role for any given realisation. As a result, it is important to study the statistics beyond simply the complex means, which we do numerically below.

A. Numerical Model
In order to further study the statistical properties of the enhancement factors, Monte Carlo simulations were performed for scattering of SPPs propagating at a metaldielectric interface (with dielectric constants ε m and ε d respectively) by nanoparticles in the dielectric near the surface (see inset of Fig. 2). As discussed above, this choice of system is motivated by the use of SPP scattering in biological sensors [11,20,74]. Specifically, realisations of randomly distributed scatterers were generated and their corresponding scattered fields calculated by solving Eq. (3) and using Eq. (2). The simulation was repeated with an additional particle (cf. Eq. (5)) from which the field perturbation and individual enhancement factors were determined. Notably, a scalar model can be used to describe SPP scattering [76,77], with the scalar field corresponding to the out-of-plane component E z of the SPP field. When both z and z are near the interface, the Green's function can be approximated as a cylindrical wave [65,77,78] given by k SPP is the complex SPP wavenumber with corresponding absorption length l abs = (2 Im[k SPP ]) −1 and H (1) 0 (x) is the zeroth order Hankel function of the first kind. Simulations were performed using this Green's function. The incident field was taken to be a decaying SPP plane wave of the form E 0,z (x) = Θ(x) exp(ik SPP x), where Θ(x) is the Heaviside step function and we assume z N +1 = z s . Evaluating Eqs. (33)-(35) with these assumptions gives where we have defined k(n) = (k 2 SPP + 4nµ) 1/2 , µ = α(k 2 0 /ε 0 )A 0 exp[−2ak SPP z s ] and µ N +1 is defined analogously with α N +1 and z N +1 replacing α and z s . In addition, the SPP elastic scattering cross section σ SPP = 4|µ| 2 / Re[k SPP ] and corresponding scattering mean free path l s = (nσ SPP ) −1 can be defined for this model [76,77]. Note that the complex incident wavevector k SPP (i.e. the presence of absorption) means that γ 3 does not take the form of Eq. (36). In order to study the role of absorption, simulations were performed at two different wavelengths. Firstly, the 'low loss' case was simulated at λ 0 = 650 nm, for which ε d = 1.77 (corresponding to water) and ε m = −13.68 + 1.04i (corresponding to gold [79]), meaning that k SPP = (1.42 + 0.008i)k 0 . The 'high loss' case corresponded to λ 0 = 600nm, for which ε d = 1.77 (water) and ε m = −8.0 + 2.1i (gold) were taken whereby k SPP = (1.49 + 0.05i)k 0 . The absorption lengths were 9.9λ 0 and 1.6λ 0 respectively. In each case, the number of scatterers N was fixed (700 for the 'low loss' case and 800 for the 'high loss' case), and they were randomly distributed in a square of sides L. To vary the scatterer density n, L was varied between L = 9.3λ 0 and L = 118λ 0 in the low loss case and between L = 8λ 0 and L = 30λ 0 for the high loss case. Different sets of parameters were chosen for the two different wavelengths in order to ensure the density ranges in each case included both the single scattering and strong multiple scattering (l s < λ 0 ) regimes. In all simulations performed, the additional scatterer was identical to the other scatterers (α N +1 = α) and added at the fixed position r N +1 = (0, 0, z s ). All data points shown were calculated using 50,000 realisations of different scatterer positions unless otherwise stated.

B. Sensitivity Enhancements: Absorption Dependence
Figs. 2(a)-(c) show the complex mean enhancements γ i observed in the far field at 70 • to the surface normal in the backward direction (k = −ε 1/2 d k 0 sin(70 • )x) for λ 0 = 600 nm and assuming a polarizability α g1 corresponding to a 21.5nm radius gold sphere sitting on the gold surface. The mean amplitudes |γ i | are also shown in Fig. 2(d). The theoretical expressions (Eqs. (40) and (41)) are seen to describe γ 1,3 well over the entire density range. Both γ 2 and γ 3 remain close to unity, as do the corresponding mean amplitudes, indicating that the effects of the associated multiple scattering paths are negligible. As a result, γ 1 is the dominant factor in the behaviour of the total mean amplitude enhancement |γ 1 γ 2 γ 3 | (Fig. 2(d)), which scales very similarly to |γ 1 | and | γ 1 | (Fig. 2(a)).
Equivalent plots for the low loss, λ 0 = 650nm, case with a polarizability α g2 equivalent to that of a 40 nm gold sphere sitting on the surface and the same observation position are shown in Fig. 3, from which a few significantly different features are evident. In the low loss case, the enhancement factors show greater deviation in the complex means from unity ( Fig. 3(a)-(c)), even at mean free paths of several wavelengths, which is unsurprising because the attenuation of propagating SPPs means the amplitude of multiple scattering paths are negligible when l s > l abs . The other significant difference between the low and high loss cases is in the mean of the absolute value of the enhancement factors ( Fig. 3(d)). The statistics of this quantity are explored in more detail in Ref. [48], but here we note that in the low loss case, |γ 1,3 | are very different from | γ 1,3 |, by up to two orders of magnitude, whereas in the high loss case, the quantities are similar in value. Importantly, the low loss case allows for mean total amplitude enhancements |γ 1 γ 2 γ 3 | > 1, implying that multiple scattering increases the sensitivity, quite significantly, for a wide range of densities, whereas in the high loss case, multiple scattering only acts to decrease sensitivity on average. For the low loss case, the analytic results (Eqs. (40)-(42)) still provide an accurate description at lower densities/longer mean free paths, however at higher densities, significant deviations are seen, particularly for γ 3 , indicating that the loop scattering paths ignored in the derivation of the average enhancements play a significant role. Such loop paths are associated with weak localisation effects such as coherent back-scattering [80], which become significant at higher densities when Re[k SPP ]l s ∼ 1. Furthermore, in the region where the mean amplitude grows large, the complex mean is slower to converge due to the larger variance in the underlying probability distribution (see Fig. 4) and hence larger statistical fluctuations are seen in the simulated data. Indeed, the results plotted between nλ 2 0 = 0.21 and 3.87 in Fig. 3 are averaged over 150,000 realisations in order to improve convergence. Stronger Anderson localisation begins to play a role at the highest densities. The localisation length ξ = l s exp (πRe[k SPP ]l s /2) [81] becomes comparable to the system size for l s ≈ 0.73λ 0 , at which point Anderson localisation means only scatterers within ∼ ξ couple strongly with each other. As a result, the effect of the added scatterer is reduced, explaining the decrease in mean amplitudes at the very highest densities.
To study the underlying probability distributions in more detail we have plotted histograms of the relative frequency of the enhancement factors in the complex plane in Fig. 4 for the low loss case at different densities. Supplementary Movies 1 and 2 show the complete density evolution of the distributions for both the high and low cases [82]. In general, γ 1 and γ 3 appear to be distributed with rotational symmetry about their centres. Specifically, the standard deviations of the real and imaginary parts were found to typically be within 10% of each other for both γ 1 and γ 3 , although in some cases large outliers can cause significant differences. Similarly, the correlation coefficient between the phase and amplitude of the centred distribution γ 1,3 − γ 1,3 was never more than ∼ 0.02 across the density range considered. In contrast, γ 2 , associated with loop scattering paths, has a more complicated locus on the complex plane, reminiscent of the previously studied eigenvalue distributions of Euclidean matrices arising in similar scattering studies [68]. In the low loss case the distributions of γ 1,3 , while being narrow at low and high density become very broad for a range of intermediate densities. Thus, although the centre of the distribution remains close to the origin, the mean amplitudes |γ 1,3 | become very large as seen in Fig. 3(d). In fact, the centre of the distributions, starting from 1 at the lowest densities, move towards the origin with increasing density. This movement of the centre of the γ 1 distribution towards the origin is also seen in the high loss case (Fig. 2(a)), however, the distribution remains tight around the centre over the full density range. Similarly, γ 3 retains the narrow width for the entire density range, although in this case the centre remains close to 1. The similarity between the mean absolute values and the absolute value of the complex mean arises from these tight distributions.
In order to understand the significant difference in the widths of the probability distributions for the high and low loss cases, we must consider the relative role of scattering and propagation phases along different multiple scattering trajectories. Each scattering path has an associated phase and amplitude which are determined by contributions from scattering events (A scat e iΦscat ) and from propagation between scattering events (A prop e iΦprop ), such that the enhancement factors are determined from the sum over all possible paths ∼ paths A scat e iΦscat A prop e iΦprop . Changing realisations changes the propagation factors while the scattering contribution for a given sequence of scatterers is unchanged, since the scatterer positions change but not their properties. When averaging over realisations, the random Φ prop leads to cancellation of the propagation component and thus the complex mean simplifies to the sum of the deterministic A scat e iΦscat factors arising from scattering events. Absorption means that scattering paths longer than l abs have a small amplitude A scat and hence contribute negligibly to the enhancement factors for that particular realisation. In the low loss case (l abs = 9.9λ 0 ), a large number of scattering paths several wavelengths long contribute. As the paths extend over multiple wavelengths, the phases Φ scat are essentially uniform and random and thus the sum over scattering paths can give a significantly different result to the complex mean. Conversely, in the high loss case, only a small number of scattering paths shorter than l abs = 1.6λ 0 contribute significantly to the enhancement factor. Furthermore, since the amplitude decay due to absorption occurs on the wavelength scale (the amplitude decays by ∼ 20% over one SPP wavelength in the high loss case compared to ∼ 2% in the low loss case), very short sub-wavelength scattering paths for which Φ prop is close to zero will have significantly higher amplitude and contribute more to the total enhancement factors. As a result, the high loss case is close to the complex mean since the propagation has little effect. The behaviour of γ 1,3 in the high loss case is therefore dominated by the scattering phase shift.

C. Optimising Enhancements: Scatterer Dependence
The conclusion that multiple scattering has a more pronounced effect in the low loss case in unsurprising, but the fact that the high loss case shows great sensitivity to the phase acquired in a scattering event, which is determined by the individual scatterer properties, is significant. To illustrate this, we consider a case close to the divergence condition of Eq. (38). of directions in which light radiated to the far field is strongly confined [83,84]) in the glass substrate satisfies k = Re[k SPP ]. An observation position in the leakage ring furthermore has the additional benefit, from a sensing perspective, that the confinement of light means detected signals are stronger. While such a thin film configuration alters the Green's function and surface dressing, the functional form of the SPP remains the same for points in the lower index dielectric near the surface of the gold film, with only the parameter values changed (i.e. A 0 , k SPP , a and α). We thus now consider such an observation position, keeping in mind that the parameters in the model will no longer correspond to the same physical properties. Fig. 5 shows the results from further simulations of the high loss case, analogous to those shown in Fig. 2, albeit assuming k = − Re[k SPP ]x and that the polarizability is phase shifted by π, i.e. α = α g1 e iπ . Note that since the amplitude of the polarizability is unchanged, the cross-section and mean free path are also unaltered. The phase shift to α alters the absorption loss from a single scatterer, and also the phase difference between the scattered and incident field. The chosen phase means the divergence condition of Eq. (38) is nearly satisfied, i.e. the phase difference between the SPP incident on a scatterer and the SPPs radiated by the scatterer is small. We see significantly different behaviour in Fig. 5 as compared to Fig. 2. In particular, an optimum density n opt = 3.18λ −2 0 , at which |γ 1 γ 2 γ 3 | is maximised, is evident with a corresponding total amplitude enhancement of |γ 1 γ 2 γ 3 | = 196. The optimum density predicted from Eq. (37) is n opt,1 = 2.22λ −2 0 . Critically, these large enhancements occur even with l s > l abs when one might expect absorption to quench the effect of multiple scattering as was observed in Fig. 2. Results for the 'low loss' case with α = α g2 e 3iπ/4 (tuned near the divergence condition for the 'low loss' parameters) were also obtained (not shown), however, in contrast to the high loss case, the behaviour of the means shows very little difference qualitatively from the results of Fig. 3 and with similar levels of enhancement observed.
For the case shown in Fig. 5 and its low loss counterpart, the probability distributions over the complex plane behave analogously to the behaviour shown in Fig. 4   l opt,tot ) at which it occurs were calculated numerically as the phase of µ was varied (|µ| was again held constant and we assumed k = − Re[k SPP ]x). While arg(µ) is not dynamically tunable in general, it can be modified by changing various properties of the scatterers, for example their composition or geometry, or tuning the wavelength through a localised plasmonic resonance. More complex engineered scatterer structures such as core-shell nanospheres or nanorods allow further degrees of freedom for tuning α. In addition, the phase of µ can be altered via its dependence on z s and use of index-matched spacer layers. Fig. 6 shows the dependence of |γ 1 γ 2 γ 3 | max and l opt,tot on arg(µ), for both the high loss and low loss case. In the low loss case, we see that |γ 1 γ 2 γ 3 | max is always achievable regardless of arg(µ), with the value varying slightly with arg(µ), albeit remaining ∼ 10 2 for a broad range of phases. The optimum phase predicted from Eq. (38) coincides with the region where the enhancement is largest, and is also achieved at larger mean free paths (i.e. lower densities). Conversely, the high loss case has a range of arg(µ) for which no enhancement is possible on average, since absorption quenches any mul-tiple scattering enhancements. Tuning of arg(µ) does nevertheless allow a similar level of enhancement to the low loss case to be achieved, with the divergence condition introduced by Eqs. (37) and (38) providing a good predictor of the optimum phase. For the low loss case, long range scattering paths play a significant role as is discussed further in Ref. [48].

IV. CONCLUSION
To conclude, we have presented a general formalism to describe multiple scattering based enhancements to the field perturbation caused by adding an analyte particle into a random distribution of background scatterers. The approach presented is general and applicable to any wave scattering scenario, both vector and scalar, through appropriate choice of Green's tensor, for example scattering of acoustic waves or electromagnetic waves in free space, waveguides or photonic crystals [85][86][87][88]. Three enhancement factors were derived, each arising from a different class of multiple scattering paths and their statistics were FIG. 6. The maximum mean amplitude enhancement |γ1γ2γ3| (blue ) and mean free path lopt,tot (red ) at which it occurs for both high (a) and low (b) cases, as a function of the phase of α (or equivalently µ) relative to that of a gold nanosphere on the surface. The observation point was taken in the leakage ring (k = − Re[kSPP]x). Median value of |γ1γ2γ3| at lopt,tot is also shown (orange ). Light blue shaded region indicates |γ1γ2γ3| ≤ 1, meaning that the single scattering case is optimum and multiple scattering always reduces sensitivity on average (these points are not plotted). The dashed black line denotes the optimum phase at which Im[α G] vanishes and Eq. (38) diverges. studied in the context of scattering of planar SPP waves. Through a series of Monte-Carlo simulations we demonstrated that absorption can play an important role in the statistics of the enhancement factors, as it can quench long distance scattering paths. Supporting analytic calculations for the complex means of the enhancement factors were found to agree well when loop contributions were negligible. Whilst absorptive quenching was often seen to lead to an absence of any multiple scattering enhancement for high loss systems, the small propagation phases of short distance scattering paths imbues the system with a greater sensitivity to the scattering phase shift, and hence the individual scatterers. Consequently, we demonstrated that, by tuning the polarizability of the background scatterers, a mean total enhancement of up to two orders of magnitude can be achieved. Analytic expressions, capable of predicting the optimum polarizability, were also derived. Low loss systems were shown to exhibit contrasting behaviour. Specifically, it was found to always be possible to achieve an enhancement through appropriate tuning of the density of scatterers, regardless of the individual scatterer properties. Our results therefore demonstrate that multiple scattering can significantly enhance single particle detection, even in the presence of high losses, whilst insights gained can aid design of random scattering based nanostructured sensors, potentially enabling detection of weakly scattering particles such as single proteins or virions.

ACKNOWLEDGMENTS
This work was funded by the Engineering and Physical Sciences Research Council (EPSRC) (1992728) and the Royal Society (UF150335).

APPENDIX
In this section we outline the derivation of the function R ± (z i , z N +1 ) appearing in Eq. (21) defined through the equation G ∞ (r, r N +1 ) −1 G ∞ (r, r i ) = R ± (z i , z N +1 )e −ik ·(ρi−ρ N +1 ) e −ikz·(zi−z N +1 ) . (43) We assume that the upper interface of a planar stratified medium (such as a thin film structure) is located at z = 0, whilst the lowest interface lies at a position z = −d. Recall that all scatterers are assumed to lie in the upper half-space z i > 0 for i = 1, 2, . . . N + 1. We first note that from the translational invariance of the Green's function in the transverse plane, i.e. G ∞ (r, r i ) = G ∞ (r, z i ) exp −ik · ρ i it follows immediately that R ± (z i , z i ) = I. Considering the more general case of observations positions lying in the lower half space, i.e. for z < −d, it also follows trivially that R − (z i , z N +1 ) = I since there is only a transmitted component of the Green's function whereby from Eq. (20) G ∞ (r, r N +1 ) −1 G ∞ (r, r i ) = G tr ∞ (r, 0) −1 G tr ∞ (r, 0)e −ik ·(ρi−ρ N +1 ) e −ikz·(zi−z N +1 ) = e −ik ·(ρi−ρ N +1 ) e −ikz·(zi−z N +1 ) .
In the reflection case, the Fourier space Green's tensor, for observation points above the source point z > z i > 0 is [73] G(k ; r i ) = i 2k z H(k , z i )e −i(k ·ρi+kzzi) where H(k , z i ) = (1 + r s (k )e 2ikzzi )Γ s (k ) + Γ p (k )D(k , z i ), (46) D(k , z i ) is a diagonal matrix given by and Γ s,p are matrices projecting the source onto s and p polarized vectors, and can be expressed Γ s,p (k ) =ê s,p (k )ê † s,p (k ).
The unit vectorsê s,p (k ) are the s and p polarized unit vectors for a plane wave of wavevector k + k zẑ , given bŷ e s (k ) = (−k y , k x , 0) T /k (49) e p (k ) = (−k x k z , −k y k z , k 2 ) T /( √ d k 0 k ).